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

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


Introduction
One of the striking characteristics of freely rising and settling particles at sufficiently large Reynolds numbers is the existence of non-vertical paths.This type of motion is a common characteristic for the dynamics of blunt bodies in general (Ern et al. 2012) and is related to the presence of a (periodically) oscillating wake with vortex shedding, akin to that forming behind a fixed geometry (Gerrard 1966;Perry et al. 1982;Williamson 1996).The phase regimes of freely rising or settling bodies are more complicated than their fixed counterparts, however, due to the inherent coupling between the motion of the body-fluid interface and the surrounding flow morphology.This results in a complex, often only quasi-periodic motion which generally is difficult to predict a priori.Therefore, a complete solution or model of these systems remains elusive in many cases such that new insights rely on empirical results and parameter studies based on experiments or numerical simulations.Understanding and modelling of particle dynamics is important in many fields, for instance to predict the the spread of (plastic) pollutants in the ocean (Sutherland et al. 2023).
Properties of the paths and dynamics observed for buoyancy/gravity driven motion are determined by the strength of the coupling between particle and fluid.When this coupling is weak (Horowitz & Williamson 2006) or when the degrees of freedom of motion are limited (Williamson & Govardhan 2004), then the fluid response will be similar to that of a fixed geometry.On the contrary, regimes exist where particle kinematics are strongly affected by the fluid motion, and vice versa leading to alterations in flow morphology and particle trajectory and dynamics, 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 et al. 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 (2D) cylinders by introducing a Centre-Of-Mass (COM) offset.This approach is motivated by recent work of Will & Krug (2021b), where the COM of freely rising and settling spheres was varied experimentally.This study found that the rotational dynamics of the spheres are strongly affected by the internal mass distribution, which in turn strongly affects all other aspects of particle motion.Stronger rotational dynamics are 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 timescale ratio between the 'pendulum' frequency, induced by the offset, and the driving frequency, which is determined by the vortex-shedding.This model captured several key features of the particle dynamics when COM offset was present.It further predicted additional dependencies on the particle-to-fluid density ratio Γ , the dimensionless Moment Of Inertia (MOI) of the particle I * , and implicitly on the Galileo number Ga, which is similar to the Reynolds number Re as it is a measure of the ratio between inertial to viscous forces but uses an a priori defined buoyancy velocity scale instead of a dynamical one, and governs the onset and transitions in between the various wake-structure topologies.However, due to experimental and physical constraints, the parameter space available in Will & Krug (2021b) is insufficient to test these conclusively.Therefore, we aim to systematically investigate these dependencies by means of numerical simulations of 2D cylinders with COM offset, to show that the underlying physics are similar between the 2D and the three-dimensional (3D) case, and that the results presented here can shed light on the remaining open questions.
The problem of freely rising or settling cylinders is an extension of the classical case of vortex-induced vibrations (Bearman 1984;Parkinson 1989;Govardhan & Williamson 2000;Williamson & Govardhan 2004), where a cylinder is placed in a free stream with only limited degrees of freedom.The applied restrictions indicated that the degrees of freedom and, therefore, the amount of particle motion (tuned by body inertia, spring stiffness, and damping of the supports) strongly influence the wake structure and coupled dynamics.Williamson & Roshko (1988) presented qualitatively similar results for a cylinder that was forced periodically in a free stream and classified the resulting wake patterns as a function of the driving amplitude and frequency.Building on this, the work by Jeon & Gharib (2004) showed that the type of vortex-shedding depends on transverse and streamwise oscillations as well as on their relative phase.Similarly, for elastically mounted cylinders at subcritical Reynolds numbers (Re ⩽ 30) the effects of forced rotations were examined in recent work by Bourguet (2023), uncovering significant alterations in flow structure and amplitude of oscillation depending on Re and rotational magnitude and frequency.For the case of freely rising and settling 2D cylinders, a critical mass density ratio (Γ ), the ratio of particle to fluid mass density, was encountered, marking the threshold between reduced coupling at high Γ where the particle dynamics and its wake barely influence each other.On the contrary, below this threshold, particles show large path amplitudes and substantial alterations in the wake vortex shedding frequency (Horowitz & Williamson 2006, 2010b).Similar density ratio related transitions in the regime of motion have also been observed for spheres (Horowitz & Williamson 2010a;Auguste & Magnaudet 2018;Will & Krug 2021a).Following the same train of thought, the rotational moment of inertia was also investigated as a relevant parameter, governing the dynamics of rising and settling 2D cylinders in the numerical work by (Mathai et al. 2017), as well as experimentally for spheres (Mathai et al. 2018;Will & Krug 2021a).Due to these previous observations, we investigate the effects of Γ and MOI separately and in combination with effects induced by a COM offset.
Before proceeding with the problem definition, caution is warranted when interpreting the results as the 2D assumption in this work effectively corresponds to the limiting case of particle motion for very long cylinders settling or rising in a three-dimensional (3D) environment.Beyond a certain Reynolds number, the flow will become inherently 3D, even for a cylinder of infinite length; for a fixed cylinder, this is found to occur around Re ≈ 190 (Henderson 1997;Williamson & Brown 1998;Aleksyuk & Heil 2023).Moreover, the cylinder length and associated end-effects play an important role (Inoue & Sakuragi 2008), such that the motion of these cylinders becomes inherently 3D, showing both horizontal (azimuthal) cylinder rotation (around the vertical axis) and fluttering motion (around a horizontal axis) (Toupoint et al. 2019).Therefore, we would like to preface this work by stating that the principle goal is not to perfectly predict dynamics of 3D 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.

Problem definition and equations of motion
In this work, we concern ourselves with the dynamics and kinematics of freely rising and settling 2D cylinders in an otherwise quiescent, infinite fluid.The fluid phase has a constant mass density ρ f and a kinematic viscosity ν.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 x p , v, and a, respectively.
The particle (see schematic in figure 1 (a)) has a circular diameter D, and an effective mass m p , mass density ρ p , and a volume V, per unit length.The COM of the cylinder (designated with point G) is displaced by a distance ℓ from the geometric centre C. Subscripts C and G throughout will refer to these points.A unit pointing vector p is defined between these points from G to C. From here, we define the angles θ between p and e y (the vertical unit vector), and θ v between p and 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.
The dynamics of freely rising and settling buoyancy-driven particles are governed by the linear and angular momentum balances.In the present work, all particle dynamics are for unconstrained motion, implying that the only two forces acting on the geometry are a body force due to gravity (F g = −ρ p Vge y ) and the force exerted by the fluid on the particle F F .For convenience we split this total fluid force into a contribution due to buoyancy and a time-varying part F F = F b + F f , where F b = ρ f Vge y .Therefore, the conservation of linear momentum for the cylinder is given by with Γ = ρ p /ρ 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 2D cylinder with COM offset reads 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 (ge y +a C )×p on the right-hand side results from the COM offset.In the current work we treat I C , the moment of inertia 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 ℓ 2 .Further, γ denotes the ratio γ = 2ℓ/D.The rotational dynamics described in (1.2) resemble those of a forced pendulum (Will & Krug 2021b), for which, when linearised, the natural frequency is given by with I * = I C /I Γ , and I Γ = m p D 2 /8 a reference value of a homogeneous cylinder with identical mass m p .Therefore, the physically relevant range is 0 ⩽ I * ⩽ 2, i.e. 0 ⩽ I C ⩽ 1/4m 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 timescale t = t/τ v , where τ v = D/V b is the vortex shedding timescale and therefore represents the typical timescale 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 = µDLV b T * f (Jordan & Fromm 1972;Bouchet et al. 2006).Here, L is the length of the cylinder which 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 a C ∼ F f /m p , and using inertial scaling Applying this to 1.2 we obtain: Here, the Galileo number is defined as Ga = |1 − Γ |gD 3 /ν.The dimensionless prefactor γ/(|1 − Γ |I * ) in front of the pendulum term is proportional to the square of the ratio T of the vortex shedding to the pendulum timescale as defined in Will & Krug (2021b).For a 2D cylinder, this timescale ratio is equal to (1.5) Note that the control parameter T is solely dependent on the prescribed particle and fluid properties.It was shown by Will & Krug (2021b) that this control parameter governs the rotation dynamics of spheres with a COM offset.We will validate and expand upon this finding in the present work.
In § 2 we will first outline the numerical approach used to obtain the results as well as the data processing applied to the data set.Then, in § 4 - § 6, our results are presented and discussed.This is split into four sections, where § 3 contains a general discussion on effects of COM offset, the resonance mechanism, and the effect of fluid inertia for freely rising cylinders.Next, in § 4 the role of the particle-to-fluid density ratio and how it affects COM-induced phenomena is discussed along with the differences between rising and settling, followed by a discussion on MOI effects in § 5 and Galileo number effects in § 6.Finally, in § 7 the primary results and findings are summarised.

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) with u the fluid velocity vector, p the kinematic pressure and f the boundary forcing from the immersed boundary method (IBM) that enforces the no-slip condition (described in detail § 2.2).Note that the hydrostatic component is absent from p as it is explicitly added to the forces acting on the cylinder.The velocity at the inflow is set to zero to simulate an asymptotically quiescent fluid.At the outflow, a convective boundary is imposed (see e.g.Kim & Choi 2006).The side walls of the two-dimensional domain are periodic.The domain size is set to 60D in the gravity direction and 20D in the transverse direction, which was found to be sufficiently large to avoid box size effects.The grid spacing was kept constant within a square of size 2D adjacent to the cylinder.Outside this domain, the mesh spacing expands linearly towards the edges of the domain.A heat equation is solved to smoothen the mesh 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.

Numerical method
The numerical scheme closely follows the immersed boundary projection method originally developed by Taira & Colonius (2007) and Lǎcis et al. (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 non-linear 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: where Du n+1 = 0 and Êu Here, N (u, v C ) denotes the non-linear operator, Ĝ the gradient operator, L the Laplace operator, D the divergence operator, Ĥ the spreading (regularisation) operator, Ê the interpolation operator, φ n+1/2 the discrete incremental pressure, p n the discrete pressure, L the Lagrangian marker coordinates with respect to geometric centre and f n+1/2 the discrete analogue of f in (2.1a).The interpolation and regularisation matrices Ê and Ĥ, respectively make use of a discrete three-point δ-function (Roma et al. 1999).
For the particle, the non-dimensional Newton-Euler equations are advanced in time via with (2.4) )/∆t and Q the matrix that interpolates the velocity inside the cylinder (cf.Kempe & Fröhlich 2012).The time step is limited with CFL = 0.4.Vector 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 v C = v G + γ/2 ω × 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 a C × p term.The terms from the latter contributions are collected in ζ n .Finally, we have for g n and ζ n : , discretised components relating to a C , respectively.
Equation (2.2a) together with constraints (2.2b), and (2.3) can be rewritten as with q n+1 = Ru n+1 , E = ÊR −1 , diagonal matrix R ≡ [∆y j , 0; 0, ∆x i ], A = M I/∆t − L/(2Ga) , and r n , r n B containing the explicit terms.We approximate the inverse of A as )/2, 0; 0, (∆y j + ∆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 A 0 0 ) Here, q * in (2.7a) is solved for via a well tested factorisation procedure (see e.g.Verzicco & Orlandi 1996).The solution of φ n+1/2 and f n+1/2 are obtained via the PETSc library (Balay et al. 1997(Balay et al. , 2019) ) using the algebraic multigrid method BoomerAMG as the preconditioner and the general minimum residual method (GMRES) to solve (2.7b).This combination of solvers was found to be robust and converge within 12 to 17 iterations depending on the selected grid size and time step.A relative tolerance was set to 10 −13 , to obtain solutions that satisfy the divergence free condition and no slip condition up to the limit of double-precision calculations.Once the solution vector is found, we update the pressure field (cf.Verzicco & Orlandi 1996) (2.8) The additional solving routines for equations (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 firstorder convergence rate in the L 2 norm for the velocity field and first-order in time (owing to the approximation of A −1 ≈ ∆tM −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 Γ = 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 v * C in r n B , since the corrector velocity v n+1 C is solved for simultaneously with q n+1 .

Data set and processing
In this work, a total of 938 cases have been simulated for different combinations of the four control parameters: T , Ga, Γ and I * .The main goal is to investigate the effect of the COM offset in combination with the other parameters.For this, we varied T between 0 and a maximum of 0.6, Ga in the range between 50 up to 2000, Γ between 0.001 and 5, and I * from 0.5 to 16. Compiled input and output parameters for all cases in our data set are included in the supplementary materials.
All results presented in this study are obtained after a statistically steady state has been reached.To ensure this, first, a moving average is computed of the time trace of the vertical velocity with an averaging window much larger than a single period of the typical fluctuations.This processed signal is compared to the terminal velocity of that case (determined by the average of the last 10% of the time signal).The initial transient is considered to have ended once the filtered time signal deviates less than 5% from this terminal velocity.For Ga = 200, this typically is the case after a transient time of 60D/V b , which is short compared to the total average simulation time of (1.9 × 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 x p (t) and θ.The reference θ = 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, Γ, T , I * ), or one full oscillation time.Thus, we obtain a 'centreline' (x cl (t), with horizontal mean drift velocity v d = ⟨|dx cl, x /dt|⟩, documented in the supplementary data) which is subtracted from the actual position and orientation time signal to remove any low frequency effects.The absolute value of the signal processed this way is used to determine a list of the individual peak amplitudes (A) for the path and (θ) for the rotational oscillations, the mean of which is denoted by Â and θ, respectively.Note that, as a consequence, this can mean that the particle does not exhibit rotational oscillations around θ = 0 (where p is pointing upwards).Instead, especially for small offsets we observed a behaviour where θ might drift away from θ = 0 followed by a large rotation back to the reference state when the rotational amplitude becomes large.
The phase lag ∆ϕ between the horizontal component of the Magnus lift force F m and the horizontal body acceleration a x is calculated via cross-correlation of these quantities.Similarly, ∆ψ denotes the phase lag between the angular acceleration α and the fluid torque T f .The lag obtained from the cross-correlation is divided by the length of an oscillation period 1/f and then expressed as a phase angle ranging from -180 • to 180 • .The respective components that define ∆ϕ are illustrated in figure 1(b).Figure 1(c

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 nondimensionalised fluid vorticity field (ω z = ∂ y u x − ∂ x u y ) along with particle tracks (black lines) for six different COM offsets in the range γ ∈ [0, 1.23] (increasing from left to right).All cases here are for Ga = 200, Γ = 0.5, and I * = 1.
The cylinder with zero COM offset in figure 2a is seen to rise almost straight, with regular vortex shedding occurring at 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 T = 0.174 (figure 2b), but beyond this value, e.g.T = 0.201 shown in panel figure 2c, remarkably different kinematics are encountered.Both amplitude and wavelength of the path oscillations are significantly larger in this case, and the wake now exhibits an irregular vortex shedding pattern as can be seen in supplementary video 1.For this case it is observed that the wake structure intermittently switches between several modes: (i) path oscillations with no significant shedding events (denoted by 0 in figure 2), (ii) one single vortex pair as in the 2S-regime per oscillation, 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 timescale pattern.This chaotic shedding pattern occurs for cases close to what we will call "resonance", where the rotational forcing induced by the path oscillations occurs at the same frequency as the inherent pendulum timescale.For even higher values of T beyond resonance, represented by T = 0.285 in panel figure 2d , we observe that the large amplitude path oscillations persist albeit with a reduced wavelength.Further, the vortex shedding returns again to an unperturbed 2S mode now with staggered vortex cores due to the strong path oscillations.Finally, with even larger offsets, figure 2(e, f ), the amplitude of the path oscillations begins to gradually reduce, returning to a state very much like that for the zero offset case (see supplementary video 1 for T > 0.3 cases).The results shown here are representative of the Γ -range where the resonance phenomenon is present.In the following, we will evaluate how this resonance behaviour depends on all of the other governing parameters.

On the importance of fluid inertia
In order to investigate the effect of fluid inertia, we first consider the mean rotational amplitude ( θ) for a constant density ratio (Γ = 0.6) as a function of the timescale ratio T as shown in figure 3(a).Focusing initially on the case Ga = 200 (corresponding to figure 2), we see that at T = 0 the rotational amplitude is small θ = 0.4 • .Introducing a small amount of offset results in a marginal increase in this amplitude up to θ = 1.4 • at T = 0.16.However, around T = 0.2, there is a strong increase reaching a maximum amplitude of more than θ = 35 • at 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 T .
When comparing results across different Ga numbers, the same characteristic behaviour is observed for all cases in figure 3(a).However, the value of T at which resonance appears (marked by the steep increase in θ) is consistently shifted towards higher values as Ga decreases.Such a variation with Ga is not surprising since the definition of T does not incorporate any viscous effects.However, for low values of Ga one would expect the Stokes layer surrounding the particle to contribute significantly to the total rotational inertia of the body and thereby also to modify the particle pendulum timescale.We can account for this effect by additionally including the rotational inertia I a resulting from the Stokes layer with thickness δ as illustrated in figure 3(b) in our analysis.As a result, the modified pendulum frequency and timescale ratio of the system become respectively.Here I * a is the dimensionless fluid inertia defined as I * a ≡ 8I a /(m f D 2 ), the ratio of the Stokes layer's rotational inertia to that of the displaced fluid.The total with c 1 as the only free parameter.We find that choosing c 1 = 2.3 results in a reasonable collapse of the resonance regime for different Ga when plotting θ against T as shown in figure 3(c).The corresponding thickness of the Stokes layer and magnitude of the added fluid inertia as a function of Ga are provided in figures 3 (d, e), respectively.For Ga = 200, the thickness of the fluid layer is approximately 0.33 particle radii and the rotational inertia amounts to about 2 times that of the displaced fluid.Beyond Ga = O(10 3 ), the value of I * a changes much more slowly, explaining the weak Ga dependence at higher Ga observed in figure 3(a) as well as in previous work on spheres (Will & Krug 2021b).Note, however, that I * a is still 0.72 of the displaced fluid mass at 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.

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 a C × p, the latter being a consequence of the COM offset.Here, we will investigate their respective role with respect to the resonance behaviour.We know from the analysis in §3.2 that the maximum in rotational amplitude is related to resonance between the vortex shedding timescale τ v and the pendulum timescale τ p .However, both the viscous and translational driving will occur at the vortex shedding frequency, such that a distinction of their effects is not possible on this basis only.Answering the question of 'who's driving' also provides insight in the effectiveness of COM offset in specific regimes of motion.
In order to untangle the effects of both contributions, simulations were performed where, after a statistically steady state had been reached, the (a C × 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 Ga = 200, Γ = 0.5, and T = 0 (grey line) and T = 0.285 (red line).At t = 0, the coupling term a C ×p is turned off for the case with offset.Figure 4 (b) displays the rotation rate ω for the same simulations.These results clearly indicate that as the coupling-term is turned off, the particle dynamics return to those of the case without offset.Note here that the pendulum term ( e y × p) is still present for t ⩾ 0, but evidently it has no effect without translational driving of the rotational dynamics.Therefore, we conclude that the rotational resonance phenomenon is linked to the translational coupling.As a consequence, we expect COM offset to have no impact on particle dynamics for cases where no horizontal path oscillations (i.e.no horizontal accelerations) are present, e.g. at low Ga.This also suggests that the resonance behaviour might also be triggered by outside periodic forcing, as would be present in a turbulent flow environment.It would be interesting to study how the settling/rising velocities of low Galileo number bodies with COM offset are affected in turbulence via this mechanism.
On the role of T f , it is further instructive to consider the phase lag ∆ψ between T f and the rotational acceleration α, which is shown in figure 4 (c) for the full range for Γ and T at Ga = 200.For zero or very small offsets, T f is driving the (weak) rotational dynamics as evidenced by α and T f being close to in phase.However, as the offset increases towards resonance and beyond, ∆ψ switches swiftly to values close to 180 • , such that the viscous torque predominantly acts as damping in these cases.In essence, these trends also hold for higher Ga.However, the dynamics become somewhat more chaotic at higher Ga, as will be shown in §6, resulting in slightly lower values of ∆ψ on the order of 120 to 150 • .

Frequency of oscillation
In the following sections, we will discuss how the effect of the COM offset varies with the density ratio.In doing so, we focus on the representative case of Ga = 200 and I * = 1.We first consider the frequency of the path oscillations (f ), as this parameter also corresponds to the frequency at which the rotational dynamics are forced.In figure 5(a), we plot the data in the form of the Strouhal number as a function of T and Γ .The marker colour in the figure indicates the exact Str obtained from the simulations as can be read from the legend, the iso-contours and background colours represent a linear interpolation of this data.The transparent white area bordered by the black dashed line indicates the region where γ > √ 0.5.This corresponds to the theoretical state where I G , the MOI of the particle around the centre-of-mass, is zero in accordance with the parallel axis theorem (I C = I G + m p ℓ 2 ) as a consequence of keeping I C ≡ m p D 2 /8 = const.(i.e.I * = 1).Furthermore, we also add a line where γ = 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.
Considering the zero offset (T = 0) cases in figure 5(a) first, we find that Str varies quite significantly from Str = 0.195 at high density ratios (Γ ⩾ 0.2) to Str = 0.127 for Γ ⩽ 0.1.This transition appears to be quite sudden, suggesting the existence of critical density ratio as previously observed for rising and settling blunt bodies (Namkoong et al. 2008;Horowitz & Williamson 2010b;Mathai et al. 2017;Auguste & Magnaudet 2018;Will & Krug 2021a).The change in the path frequency at the lowest Γ is also associated with an increase in the path amplitude as is evident from the trajectories at T = 0.15 (which resemble those at T = 0) in figure 5(b).Now let us consider the effects of COM offset for varying Γ (i.e.moving vertically in the figure).Depending on Γ , three distinct effects of increased offset can be observed.First of all, for Γ ⩾ 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 (a C × p) provides the driving.The restoring torque is proportional to |1 − Γ | −1 and therefore goes to infinity for Γ → 1.This is not the case for the driving term which scales according to Γ −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 Γ 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 T the dynamics switch between the base state (near identical to γ = 0) and the resonance state.This is best shown in the inset of figure 5(c) where we see that at low values of T for intermediate density ratios (0.2 ⩽ Γ ⩽ 0.8) Str stays constant at approximately 0.195 (upper branch).However, as the offset increases there is a sharp jump to the lower branch of Str .The upper branch corresponds to a system state with minimal body rotation and translation, and in the lower branch the vortex shedding latches on to body motion and is affected by the pendulum frequency.The cases Γ = 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 Γ ⩽ 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 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 Γ , there is no effect of increasing offset, at intermediate density ratios we see a large difference between different T , and at low Γ , we observe large path oscillations even at small/zero offset.
Beside the jump at the onset of resonance, Str also varies significantly beyond the resonance state.The isocontours of Str in this region approximately follow the lines of constant T (white dashed lines) and in particular the minimum of Str coincides roughly with T = 0.08.The correlation between Str and 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 Γ 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 T .As T becomes very large, Str appears to return to the trend at large Γ where higher density ratios have a slightly higher Str .It is further clear that while T is the relevant parameter to describe the behaviour after the transition from the low Str to the high Str state, the transition itself does not coincide with T = const., but occurs at lower values of T for lower Γ .The case with the lowest density of Γ = 0.001 spans only a tiny range in terms of 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/Γ , it will likely dominate the pendulum torque, which does not diverge for small Γ .
As mentioned at the beginning of this section, the frequency of the path-oscillations is important for the driving of the rotational dynamics through equation (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 / fp , which we show in the main panel of figure 5(c) as a function of the timescale ratio T .Curves of constant Str corresponding to the two different states are indicated by the black dashed lines.Importantly, we find that f / fp = 1 occurs around 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 Str = 0.127 to meet f / fp = 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 ⩽ T ⩽ 0.12 is indicated by a grey shaded area in the figure background throughout this work.Finally, we included the results for rising and settling spheres with COM offset from the work by Will & Krug (2021b) as black circles in figure 5(c).The good agreement with the present results suggests that the underlying physics of the resonance mechanism are indeed the same and that results and trends presented here are also relevant for spherical bodies in a 3D flow environment.

On the transition to resonance
In this section, we will discuss the transition to resonance in the intermediate Γ regime, i.e. for 0.2 ⩽ Γ ⩽ 0.8 in more detail.In the range 0.3 ⩽ Γ ⩽ 0.7, the transition from the high Str number mode to the low Str one occurs within a narrow band of T .This is also visible from figure 6(a), where the power spectra normalised by the maximum amplitude (F) of v x /V b are shown for T = 0.159, 0.195, and 0.225 in the vicinity of the transition point at Γ = 0.6 (yellow diamonds, pink pentagons and purple hexagons respectively).For both, T = 0.159 and T = 0.225, the spectra feature singular peaks only at Str ≈ 0.1 and Str ≈ 0.2, respectively.The former peak also dominates for the intermediate case T = 0.195, however, a weaker secondary peak at Str ≈ 0.21 is also seen to emerge at this offset value.Similar trends can be observed across the range 0.3 ⩽ Γ ⩽ 0.7 with varying ratios of relative peak height, suggesting that the transition between modes happens in a narrow band of T , but is not entirely discrete.
In §4.1, we mentioned that Γ = 0.2 and 0.8 were on the edges of the Γ -range for which a sharp transition to the resonance regime was encountered.We will investigate these cases in more detail here.For Γ = 0.2, the range of T where multiple modes are observed widens significantly as compare to 0.3 ⩽ Γ ⩽ 0.7.We observe multiple peaks in the spectra for 0.138 ⩽ T ⩽ 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 γ = 0 for Γ close to 0 and is similar to the cases of Γ ⩽ 0.1.However, looking at figure 5(c), we still observe that the cases at Γ = 0.2 and γ ≈ 0 follow the upper branch in terms of f /f p and Str , making the behaviour transitional between the Γ -regimes.
At Γ = 0.8, we find only a single case (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 ωD/V b in figure 6(b).For both quantities, a modulation of the amplitude but also of the frequency is evident at timescales much larger than that of the path oscillations.This behaviour stands in stark contrast to the rest of the cases which exhibit very regular periodic motion.Specifically, the parameter combinations that show this characteristic behaviour are: Γ = 0.8 with T = 0.327 and T = 0.365, and Γ = 0.7 with T = 0.411.A video showing this behaviour along with the vorticity in the wake can be found in supplementary video 2.
To quantify the frequency variations, we plot the instantaneous Str in figure 6(c) based on determining the distance between the maxima and minima in the signal v x /V b .Instances where the frequency is relatively low are indicated by the blue regions in figure 6(b,c).The corresponding particle kinematics, showing a marked reduction in the lateral amplitude, and the associated vortical structures for the low Str mode are illustrated in figure 6(d ).In this case, the cylinder rises almost vertically and the length of the attached wake is at its maximum extent.After this, in the transitional period marked in yellow, the path amplitude remains low but the frequency of the oscillations increases.Figure 6(e) shows that this is linked to the rapid shedding of the build-up attached wake, quickly reducing its length.This state is similar to the dynamics observed at higher Γ .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 Γ = 0.8 and T = 0.327 at this Galileo number and is indicative of the density regime transitions, where the dynamics exhibit signs of both regimes.These observations highlight that in the transition range, multiple states can coexist.

Drag coefficient and Magnus force
In this section, we will concern ourselves with the mean vertical velocity, i.e. the terminal rising/settling velocity, which is of particular practical relevance.We define the drag coefficient C d obtained from the time-averaged vertical force balance between the buoyancy and the drag force, given that the particle has reached terminal velocity.For a two-dimensional cylinder, this results in Note that in this definition C d solely reflects variations in the vertical velocity of the particle.When plotted as a function of Γ and T (figure 7(a)), the drag coefficient exhibits considerable variations almost up to a factor of 4 across the parameter space that are predominantly induced by COM effects.At T = 0, C d is found to be lowest for Γ = 0.2.Moving to the right in the figure, i.e. towards increasing γ, the previously ( §4.1) defined Γ regimes can again be noted.For Γ ⩾ 0.8, no increase in C d is observed as a consequence of the COM offset.However, for Γ ⩽ 0.7 the resonance behaviour manifests itself in a strong increase in C d .These trends are explicitly plotted in terms of T in figure 7 (b, c).
For Γ ⩽ 0.7, the resonance behaviour reaches a maximum for T ≈ 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 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 Γ .Beyond 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.Γ ⩾ 0.8), the drag decreases slightly.It appears that the magnitude of the decrease is inversely correlated with the mass density, resulting in larger reduction for lighter particles (figure 7(c)).This phenomenon does not appear to occur at fixed values of T or T .It is noteworthy, though, that a similar drag reduction was encountered around similar values of T for settling and rising spheres with COM offset (Will & Krug 2021b).
In previous work by Will & Krug (2021b), the connection was made between the drag increase and the maximum in the enhancement of horizontal particle acceleration (a x ) through the rotation induced Magnus lift force (F m, x ∼ −ω 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 (∆ϕ) between F m, x and a x .Examples of time series with different phase lags for three values of T are shown in figure 1(c).When F m, x and a x are in phase (∆ϕ = 0), the enhancement of the horizontal particle motion by the Magnus force is maximum.The connection between ∆ϕ and C d is established for the current data set in the inset of figure 7(d ).In the main panel of figure 7(d ), the phase lag is plotted explicitly vs. 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 ∆ϕ = 0 • around T = 0.11 for all cases where resonance is present, whereas acceleration and Magnus force are significantly out of phase (∆ϕ ≈ −60 • ) 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 ∆ϕ = 0.The good agreement with the experimental sphere data of Will & Krug (2021b), included in figure 7(d ), again underlining the fact that the resonance phenomenon in 2D is indeed comparable to that in 3D.

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 three-dimensional (3D) chaotic regime when varying the MOI of spheres (Will & Krug 2021a) and both of these quantities did not result in an adequate prediction of drag for the spiralling regime.
To add to this, we investigate how variations in C d relate to the presence and strength of rotational and translational fluctuations in the present data.To this end, we define the dimensionless fluctuating velocity 8(a), and dimensionless root-mean-squared rotation rate ω * = ⟨ω⟩ rms D/V b , shown in figure 8(b), for all rising cases at Ga = 200.Since particle velocity fluctuations are dominated by their horizontal component, they are qualitatively similar to the amplitude of the path oscillations Â/D.A striking difference between the distributions of v * and ω * concerns the region at low density ratios (Γ < 0.2) and small offsets (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 ω * , respectively.In all panels of figure 8, the marker edge colour is white if ∆ϕ < 0 and black if ∆ϕ ⩾ 0. For a subset of the markers labelled with a white border, an approximately linear scaling of C d with v * 2 can be observed in figure 8(c)).This linear range corresponds to the region of small offsets and low density ratios featuring path oscillations but importantly no rotations.These cases adhere to the scaling C d (v * ) = 2.0036v * 2 +1.1 as indicate by the dashed black line in figure 8(c).However, once rotation begins to become significant, we find that it dominates the drag behaviour.This is demonstrated in figure 8(d ), from which it is clear that for the markers with black borders C d is approximately proportional to ω * 2 .For cases beyond resonance, results are reasonably well represented by the fit C d (ω * 2 ) = 0.0024ω * 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.

Settling particles
Up until this point, the focus was exclusively on light (rising) 2D cylinders.For heavy particles, it was already demonstrated in Will & Krug (2021b) that the feedback between the Magnus lift force and particle acceleration becomes negative, effectively suppressing the resonance mode.This implies that the strong coupling between rotation and translation and the associated drag increase are absent, but not necessarily that the COM offset has no effect at Γ > 1.To explore this, the density ratio range from Γ = 1.1 up to 5 was studied with T ranging from 0 to the contour I G = 0 at Ga = 200, 500 and 700 and I * = 1.We present only the results for Ga = 200 in figure 9(a-e) since the trends for the higher Galileo numbers are similar.
Figure 9(a) confirms that the drag coefficient does vary as a function of T also for settling particles.Yet, the magnitude of the increase in C d (from around 1.2 to 1.8) is much smaller compared to that observed for rising bodies (from around 1.1 up to 4).Furthermore, the drag increase is more pronounced at larger Γ and contrary to rising cylinders, the contours of constant C d do not align well with isocontours of either T or T (dashed white lines), suggesting that the mechanism of drag increase here is not resonance related.This is further evidenced by figure 9(b), where for the the phase lag ∆ϕ is shown for the same cases as in figure 7(d ).Unlike for rising particles, there is no monotonic increasing trend between offset and ∆ϕ and no regime where ∆ϕ = 0 can be identified.In fact, the phase lag is strongly positive (between 90 • and 135 • ) 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 T /T | I G =0 .Doing so reveals that the drag is maximum for low, but non-zero, rotational inertia (around T /T | I G =0 = 0.9).Similarly to rising cylinders, the increase in drag coincides with a decrease in Str (see figure 9(c)) although this effect is much smaller here as compared to the resonance mode encountered for Γ < 1.Consistent with the trends established for light cylinders at zero offset, the lateral (figure 9(d )) and rotational (figure 9(e)) amplitudes are also elevated in this parameter region.This behaviour can also be observed in supplementary video 3 showing six values of T for Γ = 2.5.While the magnitude of all the path-oscillations remains significantly lower as those encountered in the resonance regime for Γ < 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 ω * 2 is not exactly identical.Furthermore, due to the absence of rotational-translational coupling, the observed increase in body rotation is not reflected in a similar increase in the path-oscillations.We suspect that for settling the increase in drag is primarily resulting from the rotational motion given the minute increase in the translational dynamics in this case.Finally, we observe that the magnitudes of C d , Â/D, and θ all decrease towards Γ = 1.This behaviour is consistent with the previous results for Γ < 1 and the explanation provided in §4.1, namely the divergence of the pendulum term.

Effects of varying moment of inertia
In this section, we will explore the effects variations in the particle MOI around G, which thus far has been kept constant.Since particle rotation proved to be a critical aspect in the preceding analysis, it is anticipated that the variations of the MOI will also affect the overall dynamics.We investigate this by varying the dimensionless MOI around the geometric centre (I * ) in the range I * ∈ [0.5, 16] for the cases of Γ = 0.4, Ga = 200 and T ∈ [0, 0.5].The corresponding results for C d , ∆ϕ, Str , Â/D, and θ are presented in figure 10.In these figures, physically feasible boundaries are represented by the solid black line, which indicates the line where γ = 1 and the dashed black line indicates where I G = 0.The grey shaded region marks parameter combinations beyond both these two criteria.This region is probed less extensively and therefore no linear interpolation of the data is provided there.
The results for C d as a function of I * and 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 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 T resolves this dependence.The latter can be seen from the white dashed lines in figure 10(a), but is even more evident from the inset, where the same data is plotted directly vs. T ; the maxima in C d collapse onto T = 0.11.Besides the dependence on T , our data further show that higher values of I * lead to an increased peak drag coefficient at resonance.
In figure 10(b), ∆ϕ is shown for the same data set.Identically to the results described in §4.3, we find that peak drag occurs for ∆ϕ = 0 when the system is in resonance.The phase data is the best way to asses the validity of the inclusion of fluid inertia as shown in figure 10(b), the isocontours of T almost exactly match the interpolated ∆ϕ 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 ∆ϕ, isocontours of all these quantities also line up with lines for which T = const.Also in line with the C d results, the resonance induced changes become stronger with increasing I * in all quantities considered.This behaviour can be understood by considering the systems as a driven and damped harmonic oscillator noting that when the rotational inertia of the system increases, the damping ratio decreases resulting in larger rotational amplitudes θ.This, in turn, then affects the other parameters (C d , Str and Â/D), for which the known effects of rotation become enhanced.
As an aside, we would like to remark on the zero offset case, γ = 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 matching Galileo number (Ga = 500) and overlapping range of I * and Γ .Finally, we would like to point out that the conclusion based on the present data, i.e. that rotation plays a very marginal role in affecting regime transition in absence of a COM offset, is also consistent with the experimental study on the effect of varying MOI for rising spheres by Will & Krug (2021a).

Galileo number effects
In this section we revisit the Galileo number, previously discussed in §3.2, however here we will take a broader scope and will look beyond the effects of fluid inertia.In this work seven Galileo numbers, ranging from 50 up to 2000, were examined for varying COM offset at fixed Γ = 0.6, and I * = 1.These results are presented in figure 11.Furthermore, for Ga = 500 and 700 we also varied Γ from 0.001 up to 5 identical to what was previously presented for Ga = 200.The only difference in the results for higher Galileo number settling particles is that for 0.16 < T < 0.19 a significant horizontal drift is encountered v d /V b > 0.1, resulting in oblique trajectories.Results for the drift velocity for all cases can be found in the supplementary data.
In figure 11 (a) the drag coefficient is shown as a function of Ga and 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 γ = 0, a minimum drag (C d = 1.2) occurs for Ga = 200.The increase in C d towards lower Ga is related to the viscous dominance in this regime and for higher Ga the increase in C d is associated with increasing path and rotational oscillations.Furthermore, the isocontours of 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 T .This also reveals that the onset of resonance occurs near at constant T .Additionally, the range of offsets where C d is affected extends to larger T for increasing Ga in a similar way as decreasing Γ or I * would.
Figure 11 (b) shows ∆ϕ for the same parameter range.Here, the T = const.contours match with constant ∆ϕ predominantly over the range −45 • < ∆ϕ < −15 • (0.08 < 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 ⟨ θ⟩, figure 3 (c).An alternate choice would have been to fit c 1 to align ∆ϕ = 0 across all Ga.Doing so results in a value c 1 of approximately 0.5 (note that this yields a more than 80% reduction in I * a for Ga < 2000), resulting in T contours overlapping the ∆ϕ = 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 T , Γ , I * , and Ga and is, in fact, not constant in time altogether.
For all Galileo numbers a reduction in the Strouhal is encountered around 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 • at high Galileo numbers are sufficient to cause a drastic drop in Str .In figure 11(e) the path-amplitude response is depicted.The behaviour here is similar to that of C d and θ except for the large increase in path-amplitude at high Ga and small values of T .We can see that as Ga increases, the path amplitude at γ = 0 also grows for Ga ⩾ 500, whereas the rotational amplitude remains low in the same parameter region (see figure 11 (c).This behaviour is akin to that observed for spheres as reported in Auguste & Magnaudet (2018); Will & Krug (2021a), where a Γ threshold is encountered demarcating the transition between a vertical and chaotic rise mode, notably in the absence of strong rotation.The Γ value of this threshold was shown to increase with Ga which is what we find here as well since for Ga = 200 this transition was encountered for 0.1 < Γ < 0.2, see §4.1.And indeed, for Ga ⩾ 500 the behaviour of the 2D cylinders is chaotic with large fluctuations in both path and rotational amplitudes without periodto-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 ( Â/D) for Ga = 700 in figure 12 (a) as well as the standard deviation of Â′ /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 Ga = 200, for which Â′ /D ⩽ 0.1 even for the transitional cases discussed in §4.2.This stresses that behaviour at higher Ga is in fact chaotic and it is therefore quite remarkable that the mean quantities are reasonably well behaved and in-line with results for lower Galileo numbers.This can in part be explained by a second observation pertaining to figure 12 (b), namely that when the resonance threshold is exceeded, i.e.T > 0.11, the value of Â′ /D drops drastically.The periodicity imposed by the pendulum frequency becomes dominant and stabilises the chaotic motion resulting from the body-fluid interaction.Qualitatively, this chaotic behaviour at low offset and the stabilising effect of large offsets is visible in the supplementary video 4 and supplementary video 5 showing results for Ga = 500 and 700, respectively.
Finally, we focus on the lowest Galileo number (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 γ = 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 ( Â/D = 0.06 at γ = 0) combined with COM offset are enough to trigger resonance behaviour similar to that observed for the 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 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 γ = 0 exhibits path-oscillations.

Summary and conclusions
In this work, we have systematically studied Centre-Of-Mass (COM) offset effects for a freely rising or settling cylinder in a quiescent fluid via direct numerical simulations employing the Immersed Boundary Projection Method.The non-dimensional parameter characterising the COM offset is given by the timescale ratio (1.5) T ≡ τ v /τ p , with τ v defining the vortex shedding frequency timescale (set by the buoyancy velocity and particle diameter), and τ p a timescale originating from the "pendulum"-like restoring torque resulting from the offset between the centres of mass and buoyancy (cf.Will & Krug 2021b).The main goal of this work has been to confirm that the behaviour of COM offset can be predicted in term of this timescale ratio, which depends on both Γ and I * .Additionally, a dependence on the Galileo number is expected but this is not reflected in the definition of T by Will & Krug (2021b).These dependencies could not be adequately confirmed experimentally in previous work due to physical constraints, therefore, a numerical study was desirable since then one can exactly set all of the control parameters.Simplification to 2D, i.e. cylinders, allowed us to examine the 4-dimensional control parameter space which is not feasible for 3D.The thus studied parameter space ranges 0 ⩽ T ⩽ 0.6 for COM offset for Galileo numbers ranging from 50 up to 2000, 0.001 ⩽ Γ ⩽ 5, and 0.5 ⩽ I * ⩽ 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 Γ , constant (or large) Ga, and constant I * appeared to be well described by T (as was the case for Will & Krug (2021b)).However, for lower Γ or Ga we found that an additional effect was playing a role.This was identified as the contribution of rotational fluid inertia (I a ).We modelled this contribution as an annulus surrounding the cylinder, the thickness of which scales according to a boundary layer (1/ √ Ga, which is identical to 1/ √ Ga when C d is constant).This hence introduces a Galileo (or Reynolds) number dependence in the definition of the timescale ratio resulting in (3.1b), that was previously not explored.For rising cylinders, for which the resonance phenomenon is present, this modified timescale ratio ( T ) was seen to capture the trends in the observed resonance behaviour with resonance occurring at 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 T = 0.11.But more importantly, this also affects the frequency of oscillation, path-amplitude, and drag coefficient (terminal velocity).Altering the COM a couple of percent can induce a more than three-fold increase in C d .This increase in drag can be attributed to an increase in both rotational and path oscillations.However, the effect of rotation typically is way more significant and instances exhibiting translational oscillations without rotation featured only relatively small drag increases.While 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 Γ 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 a C × 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 Γ > 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 Γ > 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 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 Γ is around unity.This can be explained by the fact that the ratio of the pendulum torque over the driving torque (Γ/|1 − Γ |), in (3.1b) goes to infinity for Γ → 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 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.verifies that the velocity field is divergence-free.Next, let e = Eq n+1 − (v n+1 C − ω × L) define the error vector of the no-slip condition.
The maximum error of the velocity boundary condition max |e| is presented for various times steps in figure 13.To put the results in perspective, we compare them to the findings of Breugem (2012) who studied the accuracy of the multi-direct forcing for a fixed sphere.Overall, results for L ∞ are found to be of O(10 −14 ) or lower, which confirms that the pressure solver accurately obtains the solution φ n+1/2 such that q n+1 satisfies the noslip condition.Achieving a much more accurate no-slip condition up to machine precision is a clear advantage over the multi-direct forcing scheme, but is computationally more demanding.

A.2. The order of convergence
In the following we assess the convergence properties of the fluid solver coupled to the particle equations of motion.Here, we select a two-dimensional problem of a cylinder under the influence of gravity in a viscous fluid, which is initially at rest.The domain is square and has a width of 12.288D.The particle-to-fluid density ratio is set to Γ = 1.1, and the Galileo number to Ga = 100.We position the cylinder at the centre of the domain.For each case, we use constant grid spacing and time steps.
First, we investigate the spatial and temporal convergence of the vertical velocity field.We start with the spatial convergence by fixing the time step to ∆t = 1.0 × 10 −4 D/V b , and vary the grid spacing ∆x ∈ [0.016, 0.256].After a time period of 6D/V b we compare the vertical velocity field u y for each case to that of the reference case (∆x = 0.008).We present a snapshot of the latter in figure 14(a).Note that at this point in time, the horizontal particle velocity is still negligible (v x /v y = O(1 × 10 −12 )).Let e i define the error between the vertical velocity field of a simulated case and the reference case, i.e. e i = [u ref y (i) − u y (i)]/V b , with u ref y (i) the reference velocity interpolated onto the coarser grid using a third order interpolation method.Note that u ref y (i) and u y (i) are rewritten from a two-dimensional field to the vector form.In this study, we shall make use of the L 2 and L ∞ norm to report the convergence rate of the implemented IBPM.For this, we define a normalised non-negative vector norm as with p = 1, 2, . . .an integer, and Ω 0 the corresponding domain surface.For the L 2 -norm and L ∞ -norm it is readily obtained that when ∆x i = ∆y i = ∆x, and N the number of entries in the field.Figure 14 presents the spatial convergence results.We observe that the L 2 error converges at a rate of O(∆x), whereas the L ∞ norm becomes of O(∆x) for small enough grid spacings only.These results for the convergence are consistent with those reported in Stein et al. (2017) for a two-dimensional test problem.Note that Lǎcis et al. (2016) employed ∥x∥ 2 = √ x T x/N instead of the definition for ∥x∥ 2 provided in (A 2), where the present definition is preferred because it is consistent between continuous and discrete forms.This difference explains the somewhat faster convergence rate reported in their case.
Next, we investigate the temporal convergence for u y , by selecting a grid spacing of ∆x = 0.016D for each case, and vary the time step ∆t ∈ [1.0 × 10 −3 , 1.0 × 10 −2 ]D/V b .A reference case is chosen with the same grid spacing and a smaller time step of ∆t = 1.0 × 10 −4 D/V b .At t = 0.9D/V b the convergence results are assessed.Note that the reference case has reached a velocity of v y ≈ 0.32V b at this point in time.The results of L 2 and L ∞ for the temporal convergence test are depicted in figure 14(c).Here, we observe first-order convergence for both norms.This result is expected due to the approximation of A −1 ≈ B 1 , where the leading term of B is of ∆t and is in agreement with that of Lǎcis et al. (2016) who observed the same temporal convergence for a freely falling cylinder.Now, we turn our attention to the spatial and temporal convergence of the pressure field.The datasets for this analysis are the same as those used for the vertical velocity field.We start with the analysis of spatial convergence.Figure 15(a) presents the pressure field of the reference case at t = 6D/V b for the spatial convergence analysis.Let e i = (p ref (i) − p(i))/V b , define the error with p ref the interpolated reference data and p the pressure of a coarser grid.Here, we note that the pressure nodes that reside in the particle's interior are excluded from e i to measure only the convergence from the pressure field surrounding the particle.By doing so, we find the spatial convergence to be of O( √ x) for the L 2 norm and zeroth-order for the L ∞ norm.The absence of convergence for L ∞ was also observed by Stein et al. (2017) for a cylinder in a prescribed flow.They the lack of smoothness (in the vicinity of the particle) as the main cause for the reduced convergence and proposed a forcing scheme that makes the solution globally smooth.The latter approach is promising, as the grid only needs to be uniform in a small neighbourhood of the boundary, with the same width as the regularised δ-function, but is not implemented here.
Figure 15 presents the temporal convergence, where we observe L 2 and L ∞ to be of first order, owing to the approximation of A −1 = O(∆t), similarly as observed for the velocity field.

A.3. Validations for a freely moving cylinder
In this analysis, we validate the implemented IBPM against the case of a freely rising or settling cylinder.The reference data is taken from Namkoong et al. (2008).In that study, an implicit coupling approach within a finite element framework was used, with a body fitted mesh and adaptive refinement to resolve the cylinder wake.The settling particle has a particle-to-density ratio of Γ = 1.01 and the rising of Γ = 0.99.Our computations are performed on a domain size of 20D ×60D.The grid spacing is constant in the vicinity of the cylinder and complies with D/∆x = 50.
In figure 16(a), we present our results of Ga versus Re t , where Re t ≡ V t D/ν is obtained from the terminal vertical velocity V t .Here, we observe a good agreement with Namkoong et al. (2008) for the falling and rising cylinders.Next, we compare the corresponding Strouhal number based on the mean vertical velocity of the particle V t for the same data set.We find the agreement of Str to be very good as well.

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

B.1. History torque of a sphere
Here, we will present the main results for the history torque of a sphere, which is the dominating contribution of the Stokes layer (see e.g.Auguste & Magnaudet (2018)), and then find an approximate analogue history torque for cylinder.For small time t the leading contribution of the history torque for a sphere is (see e.g.Feuillebois & Lasek 1978) with the rotational velocity ω assumed to be continuously differentiable.Discretising the integral in (B 1) for small t one finds Our goal is to find the history torque T h acting on a cylinder.Here, we assume that T h approximately takes a similar form as that of a sphere described in § B.1.To this point we consider a rotating cylinder in a viscous fluid that is at rest initially.The flow is assumed to be axisymmetric and squares of the velocity field are neglected, yielding the expression for the azimuthal velocity component ûθ (non-dimensionalised with length scale D and velocity scale ωD) with Re θ ≡ 0.25ωD 2 /ν.Here, we solve (B 5) numerically, by using the fourth order Runge-Kutta scheme for the time-discretisation and a second order central difference scheme for spatial gradients.In addition, we point out that ω was given a fixed value when solving (B 5).We selected multiple Re θ ≈ O(1) and integrated up to times such that t ≪ D 2 /ν (the used time interval ranged from 10 −4 D 2 /ν down to 10 −8 D 2 ν).In this time-interval we assume the history torque to be dominating the torque on the cylinder.The torque on the cylinder in (B 6) has the fitting parameter c 1 , which we find by calculating the actual torque on the cylinder from our numerical experiment via T = 0.25D 2 µω 2π 0 [∂ r ûθ − ûθ /r] r=D/2 dθ.Alternatively, the analytical expression of the viscous torque derived by Mallick (1957) may be used: = 4πµ (ω 0 − ω)a 2 + ωa 3 π ∞ 0 y 1 (ax)J 0 (ax) − y 0 (ax)J 1 (ax) J 2 1 (ax) + y 2 1 (ax) e −x 2 νt dx , (B 7) with ω 0 the initial rotational velocity of the cylinder, a the radius of the cylinder, and J n , y n Bessel functions of the first and second kind, respectively.Given that the initial field is quiescent, we have the initial cylinder rotational velocity ω 0 = 0.
In the limit of time t → 0, we find for the analytical expression and numerical experiment that T h = T when the fitting parameter takes the value of c 1 = 1.
(B 8) Now that we have an approximate form of the torque on the cylinder, we can find the analogue form of (B 4) for a cylinder.We then plug the obtained expression for the torque in the angular momentum balance and find We fitted c 1 such that the rotational data depicted in figure 3(a) collapses with respect to T * (presented in 3b).For this fit we found c 1 ≈ 2.3 to yield good results for cylinders.
To compare with the history torque value from the analysis in § B.2, we follow Mathai et al. (2018) and set the timescale as half the oscillation time yielding ∆t = 0.5Str D/V b .From this it follows that δ = (2πStr Ga) −1/2 .By plugging this relationship into (B 11) and assuming Str ≈ 0.11 (the value for resonance where rotation is dominant), we find a value c 1 = 2.4, which does match the leading term of our fit (c 1 = 2.3) in (B 13) surprisingly well.In previous work by Mathai et al. (2018) only the leading order term was taken into account, however the inclusion of the additional higher order terms in (B 13) results in a large mismatch in the actual value of I * a ; more than a factor 2 at a Galileo number of 50.However, we found that inclusion of the higher order terms resulted in a better collapse of the data presented in the current work, their inclusion is therefore recommended.

Appendix C. Comparison to previous work
For the cases at Ga = 500 with zero offset, we examined a parameter space spanning Γ and I * that was already explored extensively in the work by Mathai et al. (2017).Here, we take a closer look into the differences between that work and the present one.
We compare our result for Ga = 500 with γ = 0 and Γ ranging from 0.001 to 0.99 to results of Mathai et al. (2017) at identical parameters.A note on the difference in convention: in the work by Mathai et al. (2017) the parameter I * is equal to I * Γ in the present work.We extracted the results from Mathai et al. (2017, figures 2 a, b), where, due to the difference in the definitions of the nondimensional rotational inertia, our results lie on the diagonal m * = I * , i.e. the line from (A) to (D).One of the main findings in Mathai et al. (2017) was change in the vortex shedding mode from a 2S mode at high Γ and I * to a 2P mode at low values of these two parameters.We found no such transition as all cases in the present work exhibited a 2S mode.Furthermore, Mathai et al. (2017) reported that this transition was accompanied by large increases in the path and rotational amplitudes, Â/D and θ, respectively.A direct comparison between the results for these two parameters is presented in figure 17 (a, b).While the general trends of a gradual increase for decreasing Γ (i.e.low I * Γ ) for both amplitudes is consistent between the works, there are large deviations in the magnitudes of both translational and rotational amplitudes for identical cases, especially towards lower density ratios.
These differences for identical parameters combinations raise the question if their employed virtual mass approach (Schwarz et al. 2015) could explain these variations.The use of the virtual mass approach (VMA) was required to stabilise the explicit scheme in Mathai et al. (2017).We tested this hypothesis by modifying equations (1.1), (1.2) to include a virtual mass contribution on both sides scaled by coefficient C v (for which we discretised the added time derivatives on the right hand side with a forward Euler scheme) A value of C v = 0 corresponds to the present approach, while typical values to stabilise explicit schemes are on the order of the added mass, i.e.C v = 1 for a cylinder (Schwarz et al. 2015).We did not observe appreciable changes in the particle dynamics when varying C v in the range C v ∈ [0, 5] and it therefore appears that the discrepancies between our work and Mathai et al. (2017) are not related to the VMA and may be caused by other unknown factors.

Figure 1 .
Figure 1.(a) Schematic of the cylinder with the centre-of-mass (G) displaced by a distance ℓ from the volumetric centre (C).The pointing vector p is a unit vector in the direction from G to C and θ is the angle between p and the vertical (y-direction).The forces acting on the body are buoyancy (F b ) and the remaining fluid forces F f (in C) and gravity 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 ∆ϕ.(c) Time signals of the horizontal component of the Magnus force (F m ∼ −ωzvy, blue) and horizontal particle acceleration (ax/g, red) for three different offset cases, showing different phase lags.Note that, since ⟨vy⟩ = O(1), the value on the left y-axis is indicative of body rotation.
) provides three examples of signals with varying T showing a negative, zero and positive value of ∆ϕ.

Figure 2 .
Figure 2. Snapshots of particle trajectories and wake structures of rising cylinders with Galileo number Ga = 200 and density ratio Γ = 0.5 for six different centre-of-mass offsets γ (and T ) (a-f ).The offset increases from left to right as indicated by the listed parameters in the top left of each subfigure.Particle trajectories are indicated by the black lines, the grid spacing has dimensions of the particle diameter D. Coloured contours represent the normalised vorticity field (ωzD/V b ).

Figure 3 .
Figure 3. (a) Mean rotational amplitude θ as a function of Galileo number versus the timescale ratio T , here Γ = 0.6 and I * = 1.(b) Schematic showing the parameters of the fluid inertia model.(c) θ for the same cases as in (a) plotted against the modified timescale ratio T , which includes the effects of a Galileo number dependent added fluid inertia as per (3.1b).(d ) Thickness of the fluid inertia layer δ and (e) added inertia as a function of Ga, based on empirical collapse of the data.

Figure 4 .
Figure 4. (a, b) Results for Ga = 200 and Γ = 0.5 for two cases; one without offset (grey line) and one with offset (red line).(a) Dimensionless horizontal velocity of the cylinder (vx) and (b) dimensionless rotation rate (rad.)versus dimensionless time.During these runs at t = 0 the aC × p term for (1.2) is turned off, showing that in absence of this coupling term the dynamics of particles with offset almost completely revert back to those of particles without offset.(c) Phase lag ∆ψ between the rotational particle acceleration (α) and the viscous torque (T f ).

Figure 5 .
Figure 5. Results for the path oscillation frequency (f ) of rising particles at Ga = 200 as a function of T and Γ .In (a), the marker colour indicates the dimensionless frequency (Str = f D/V b ) according to the colour bar provided below.The marker type indicates the different regimes in terms of the resonance behaviour discussed in the following.The isocontours are based on a linear interpolation of the data.Dashed white lines represent isocontours of T , the timescale ratio including effects of fluid inertia.(b) Horizontal particle position over the cylinder diameter (x/D) as a function of dimensionless time grouped in three values of T for three values of the Γ 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 fp vs. the timescale ratio T .Here the marker colour indicates Γ 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 Γ 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. T .

Figure 6 .
Figure 6.(a) Single sided amplitude spectrum (F) based on the particle horizontal particle velocity (vx/V b ) normalised by the maximum amplitude.(b) Dimensionless horizontal velocity and rotation rate (rad./s)versus time for Γ = 0.8 and T = 0.327.We see that the dynamics exhibit a cyclical behaviour on a timescale much greater than that of the vortex-shedding dynamics.This behaviour is split into three parts as indicated by the colours in the background of the figure.(c) Instantaneous Strouhal number as a function of time for Γ = 0.8 and T = 0.327, calculated based on the peak-to-peak times of |vx/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.

Figure 7 .
Figure 7. (a) Particle drag coefficient as a function of the particle-to-fluid density ratio Γ and the timescale ratio T for rising particles at Ga = 200 and I * = 1.(b) Drag coefficient plotted explicitly versus T , the timescale 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 ∆ϕ between the horizontal Magnus force and the horizontal component of instantaneous acceleration versus T alongside the experimental results for spheres (Will & Krug 2021b).The inset shows the correlation between C d and ∆ϕ.

Figure 9 .
Figure 9. Results for settling (Γ > 1) 2D cylinders at 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 T .

Figure 10 .
Figure 10.Investigation on the effect of the dimensionless moment of inertia I * in combination with the timescale ratio (T ) for Ga = 200 and Γ = 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 γ = 1 and IG = 0.In (a) the inset shows the same data as the main panel, however is plotted against the modified timescale ratio T to include effects of the rotational added mass due to the Stokes layer.

Figure 11 .
Figure 11.Exploration of the combined effects of COM offset and Galileo number for fixed density ratio Γ = 0.6 and dimensionless moment of inertia I * = 1.With (a) showing the drag coefficient C d and the inset highlighting the scaling in terms of T , (b) the phase lag ∆ϕ, (c) the mean rotational amplitude θ, (d ) the Strouhal number Str , and (e) the mean path-amplitude Â/D .The white dashed lines indicate isocontours in T , the black dashed line indicates γ = √ 0.5.

Figure 12 .
Figure 12.(a) Amplitude of the path oscillations for Ga = 700 and I * = 1 as a function T and Γ .(b) For the same parameter space shows the standard deviation of the peak path-amplitude.(d-e) Trajectories and wake structure for Ga 50, Γ = 0.6, I * = 1 and T = 0, 0.3, and 0.6, respectively.The colour gradient in the wake indicates non-dimensional fluid vorticity (ω f D/V b ) as indicated by the colour bar.

Figure 14 .
Figure 14.Convergence analysis for the vertical velocity field around a settling cylinder.(a) Snapshot of vy/V b at time instance t = 6D/V b for the reference case (∆x = 8 × 10 −3 ) of the spatial convergence analysis.(b) Spatial convergence.(c) Temporal convergence with constant ∆x = 1.6 × 10 −2 for each case, and a reference temporal spacing of ∆t = 1 × 10 −4 .

Figure 15 .
Figure 15.Convergence analysis for the pressure field around a settling cylinder for the same problem as in figure 14.(a) Snapshot of vy/V b at time instance t = 6D/V b for the reference case (∆x = 8 × 10 −3 ) of the spatial convergence analysis.(b) Spatial convergence.(c) Temporal convergence with constant ∆x = 1.6 × 10 −2 for each case, and a reference temporal spacing of ∆t = 1 × 10 −4 .

Figure 16 .
Figure 16.Comparison of the present results for a freely falling or rising cylinder with Namkoong et al. (2008).
2) with ∆ω ≡ ω n+1 − ω n , and the difference between time levels n and n + 1 being equal to ∆t.If one then introduces δ = ν∆t/(πD 2 ) (B 3) one can write for the rotational equation (cf.Auguste & Magnaudet 2018) History torque of a cylinder expression for continuously differentiable ω we apply Duhamel's principle and find (by approximating ∆ω = ω n+1 − ω n ) (B 10) teaches us that the added inertia I * a takes the form ofI * a = 16δ.(B 11)B.3.Dataset fit and comparisonThe moment of inertia for an annulus is obtained via r 1 = 0.5D and r 2 = (0.5 + c 1 / √ Ga)D, and c 1 a constant.Plugging in the latter radii and calculating I * a ≡ I a /I f yields

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

Table 1 .
Overview of the grids.The first column denotes the Galileo number Ga.The second column represents the number of grid points per diameter of the cylinder.The third column is the grid resolution for the fluid phase.