# Fluttering motion of a confined cylinder falling freely in fluid at rest

Dylan Letessier¹, Véronique Roig¹ and Patricia Ern¹,†

¹Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse and CNRS, 31400 Toulouse, France

(Received 2 May 2024; revised 9 October 2024; accepted 22 November 2024)

We investigate experimentally the planar paths displayed by cylinders falling freely in a thin-gap cell containing liquid at rest, by varying the elongation ratio and the Archimedes number of the cylinders, and the solid-to-fluid density ratio. In the investigated conditions, the oscillatory falling motion features two main characteristics: the mean fall velocity $\overline{u}_v$ does not scale with the gravitational velocity, which overestimates $\overline{u}_v$ and is unable to capture the influence of the density ratio on it; and high-amplitude oscillations of the order of $\overline{u}_v$ are observed for both translational and rotational velocities. To model the body behaviour, we propose a force balance, including proper and added inertia terms, the buoyancy force and vortical contributions accounting for the production of vorticity at the body surface and its interaction with the cell walls. Averaging the equations over a temporal period provides a mean force balance that governs the mean fall velocity of the cylinder, revealing that the coupling between the translational and rotational velocity components induces a mean upward inertial force responsible for the decrease of $\overline{u}_v$. This mean force balance also provides a normalization for the frequency of oscillation of the cylinder in agreement with experimental measurements. We then consider the instantaneous force balance experienced by the body, and propose three contributions for the modelling of the vortical force. These can be interpreted as drag, lift and history forces, and their dependence on the control parameters is adjusted on the basis of the experimental measurements.

**Key words:** flow-structure interactions

## 1. Introduction

Isolated bodies with widely different sizes, weights and shapes may describe complex paths, in many instances oscillatory paths, when rising or falling under the effect of

†Email address for correspondence: ern@imft.fr

© The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
buoyancy in fluid otherwise at rest. Familiar examples of non-straight paths are provided by paper cards or autumn leaves falling in quiescent air, or by bubbles rising in still water. In this problem, the geometrical anisotropy of the body strongly influences the coupling between the motion of the body and the hydrodynamical loads that it experiences due to the motion generated in the surrounding fluid (Ern et al. 2012). A better knowledge of the kinematics and dynamics of a freely moving body is, in particular, desirable for fibres, rods or finite-length cylinders (having diameter *d* and length *L* > *d*), as this geometry is present in industrial and environmental situations spanning a large range of solid-to-fluid density ratios, elongation ratios *L*/*d* and Archimedes numbers.

The variety of different modes of path instability of fluid-solid systems involving non-isotropic bodies comes to the fore in theoretical investigations using global stability analysis in the case of short-length cylinders (Tchoufag, Fabre & Magnaudet 2014) and numerical simulations regarding oblate spheroids (Zhou, Chrust & Dušek 2017; Moriche, Uhlmann & Dušek 2021). The consideration of faceted objects (Gai & Wachs 2024) and of two-dimensional non-homogeneous cylinders involving centre-of-mass offsets (Assen et al. 2024) also recently enlarged the understanding of the role of rotational effects on the dynamical behaviour of complex bodies. Concerning elongated bodies freely falling or rising at moderate Reynolds numbers in three- or two-dimensional configurations, experimental investigations also revealed the diversity of paths displayed by these bodies when the control parameters are varied (see D'Angelo et al. 2017; Toupoint, Ern & Roig 2019; Will et al. 2021, and references therein). Different types of periodic motion emerge in association with specific unsteady wakes, sharing specific symmetry properties. The elongation ratio of finite-length cylinders has in particular been shown to have a noticeable influence on the wake structure, its stability and symmetry properties, when the body is free to move (Lorite-Diez et al. 2022) as well as when it is fixed in an incident flow perpendicular to its revolution axis (Inoue & Sakuragi 2008; Yang, Feng & Zhang 2022).

In experiments, attention is focused principally on the characterization of paths presenting displacement deviations from the rectilinear motion that are comparable with the largest dimension of the body, typically *L* for the finite-length cylinder investigated here. For densities close to that of the fluid, and close enough to the transition from rectilinear to periodic motions, scaling laws characterizing the cylinder kinematics have been found (Toupoint et al. 2019). For a larger range of density ratios, in particular heavier bodies, the description of the body kinematics resorts to correlations presenting complex dependences on the control parameters, associated with the drag force experienced by the finite-length cylinder (see e.g. Clift, Grace & Weber 1978; Chow & Adams 2011; D'Angelo et al. 2017, and references therein). As discussed in the present paper, the difficulty arises from the significant oscillation amplitudes reached by the translational and rotational velocities of the body, which make the simple balance between drag and buoyancy forces along the vertical incomplete. Consequently, the gravitational velocity scale, commonly relevant for oscillating paths displaying velocity fluctuations one order of magnitude smaller than the mean vertical velocity (yet with centre-of-gravity displacement oscillations of the order of the body dimension *L*), no longer provides the order of magnitude of the mean fall velocity. Here, we investigate a representative case of such paths, featuring translational and rotational velocity fluctuations comparable in magnitude to the mean fall velocity experienced by the body. As will be shown, nonlinear forces coupling the fluctuating degrees of freedom of the body then have to be considered in the force balance governing the body behaviour.

Determining the hydrodynamical loads acting on finite-length cylinders (diameter *d* and length *L*) is a demanding task, even for isolated bodies fixed in an incident flow. As the problem depends on three dimensionless parameters – the Reynolds number *Re*,
Fluttering motion of a confined freely falling cylinder

the cylinder’s elongation ratio *L*/*d*, and the orientation of the incident flow relative to the cylinder’s revolution axis – the exploration domain is vast, and only a few inertia flow configurations have been presently investigated for fixed finite-length cylinders (Pierson *et al. 2019*; Kharrouba, Pierson & Magnaudet 2021; Pierson, Kharrouba & Magnaudet 2021; Cabrera *et al. 2022*). For freely moving cylinders, a dynamical model embodying their complex paths is still to be devised and adjusted, but can be searched on the basis of the investigations performed for the related case of the paths displayed by freely falling plates. Empirical models for these paths have been proposed in specific cases based on aerodynamic force models. Most often, the forces considered take the form of drag and lift contributions related to the instantaneous degrees of freedom of the body, for instance, as performed by Belmonte, Eisenberg & Moses (1998) for confined falling plates. Andersen, Pesavento & Wang (2005) also elaborated a dynamical model involving a quasi-stationary drag force and a circulation-related lift force to reproduce the fluttering and tumbling motions of thin two-dimensional plates falling freely in a three-dimensional (3-D) tank. The model successfully describes the different trajectory types observed experimentally, but focuses on one density ratio, and proposes separate coefficient adjustments for the four different elongation ratios considered. Extensions for this model have further been sought for heaving and pitching aerofoils following a prescribed motion (Moriche, Flores & García-Villalba 2017) and also, recently, for freely falling plates with appendages allowing for variable centre-of-mass offsets (Li *et al. 2022*). Along the same lines, the purpose of the present paper is to investigate experimentally the oscillatory motions of confined freely falling cylinders, and to propose a model for these paths that is able to take into account their dependence on the three control parameters, the elongation ratio, the density ratio and the Reynolds number.

## 2. Experimental results

The experimental set-up consists of a thin vertical cell, having height $h = 80$ cm, width $l = 40$ cm, and thickness $w = 1$ mm. We investigate the free fall of finite-length cylinders (diameter $d$, length $L$) in liquid at rest. The cylinders are released from the upper part of the cell, initially positioned with their axis of revolution horizontal. The diameter of the cylinders is constant and equal to $d = 0.8$ mm, comparable to the gap size $w$. Conversely, the dimensions of the cell, $l$ and $h$, are large relative to the length $L$ of the cylinders. We considered cylinders with elongation ratios $\xi = L/d \in \{3, 5, 7, 10, 12, 20\}$ made of different materials (density $\rho_c$) and released in water (density $\rho_f$, kinematic viscosity $\nu_f$), corresponding to body-to-fluid density ratios $\rho_c/\rho_f \in \{1.12, 1.4, 2.7, 4.5\}$. These density ratios were obtained with the following materials: 3-D printing plastic from the Visijet brand (M2R-CL resin), polyacetal plastic POMC, aluminium, and titanium, respectively. In the present configuration, each value of the density ratio corresponds to a constant value of the Archimedes number $Ar = \sqrt{(|\rho_c - \rho_f|/\rho_f)gd}$ ($d/\nu_f \in \{25, 60, 100, 155\}$). Note that it might also be interesting to introduce the Archimedes number $Ar_{3-D} = \sqrt{(|\rho_c - \rho_f|/\rho_f)gd_{3-D}}$ ($d_{3-D}/\nu_f$), based on the diameter of the sphere having a volume identical to that of the cylinder, $d_{3-D} = (\frac{3}{2}d^2L)^{1/3}$. Cylinders with different $\xi$ and $\rho_c/\rho_f$ but equal $Ar_{3-D}$ falling in the same liquid experience the same buoyancy force, which is the driving force for their motion. The explored range of $Ar_{3-D}$ is large, varying from 75 to 720. The cylinder motion was tracked using a high-speed camera (Photron RS-3000 APX) with acquisition frequency 250 Hz. The dimensions of the field of view were 10.8 cm × 10.8 cm, the size of one pixel corresponding to 105 µm.
<figure><img src="image_1.png" /><figcaption>Figure 1. Paths of aluminium cylinders falling in liquid at rest. The elongation ratios of the cylinders are, respectively from (a) to (f), $\xi = 3, 5, 7, 10, 12$ and 20. The centre of gravity evolutions are shown in blue. The cylinders are displayed in black at constant time intervals of 40 ms.</figcaption></figure>

## 2.1. *The different types of paths*

The strong confinement of the cylinders in the gap of the cell, $d/w = 0.8$, was chosen to reduce the degrees of freedom of the cylinder to two translations in the plane of the cell and a rotation about the perpendicular direction. In fact, Gianorio *et al.* (2014) showed that this confinement ratio prevents the rotation of the cylinder about its revolution axis, while a displacement of the cylinder within the gap in association with this rotation can occur for confinement values $d/w \le 0.6$, as also observed for fixed confined cylinders allowed to move only within the thickness of the cell (Semin *et al.* 2012). Our study is therefore carried out in a configuration where such specific unsteady motions related to confinement have not been reported, so that, in the following, the cylinder motion is assumed to reduce to three degrees of freedom.

In the studied range of parameters, we observe two main types of planar trajectories, i.e. the rectilinear path and the oscillatory path, also referred to as fluttering. The rectilinear path corresponds to a constant fall velocity of the cylinder with its revolution axis horizontal, perpendicular to the velocity direction. The fluttering motion is illustrated in figure 1, showing, for various elongation ratios, in blue the trajectories of the centres of gravity of aluminium cylinders ($\rho_c/\rho_f = 2.7$) and in black the cylinders at several instants. This motion is also observed for several elongation ratios of titanium cylinders ($\rho_c/\rho_f = 4.5$), as can be seen in figure 2. It features a displacement oscillation of the centre of gravity of the cylinder in the vertical plane of the cell, coupled with a rotation of the cylinder axis about a horizontal direction perpendicular to the cell plane. These trajectories were observed previously in thin-gap cells by Gianorio *et al.* (2014) and D’Angelo *et al.* (2017) for stronger lateral confinement of the cylinders ($0.055 \le L/l \le 0.94$), and by Toupoint (2018) for experimental conditions comparable to ours. In our case, the cell width is considerably wider than the length of a cylinder, and $0.009 \le L/l \le 0.057$.

A third type of trajectory was identified and will be termed here ‘apparently chaotic’, as it is characterized by a non-reproducible and irregular behaviour. We observed this phenomenon in the case of titanium cylinders having intermediate elongation ratios $\xi = 5$.
<figure><img src="image_2.png" /><figcaption>Fluttering motion of a confined freely falling cylinder</figcaption></figure>

Figure 2. Paths of titanium cylinders falling in liquid at rest. The elongation ratios of the cylinders are, respectively from (a) to (f), $\xi = 3, 5, 7, 10, 12$ and $20$. The centre of gravity evolutions are shown in orange. The cylinders are displayed in black at constant time intervals of 28 ms. Chaotic paths for $\xi = 5$ and 7 are extracted from the data of Toupoint (2018) and displayed with a different scale.

and 7 (figure 2). These findings are in line with those reported by Toupoint (2018). In this regime, the inclination of the body axis becomes from time to time so significant that the cylinder flips over (it tumbles). However, the ‘apparently chaotic’ trajectory is different from the tumbling motion observed by Belmonte *et al.* (1998) for thin plates confined similarly, as the tumbling motion is a predictable path associated with a periodic flipping over of the body, which drifts laterally in association with the constant-sign rotation. In addition to the fluttering and tumbling behaviours, similar ‘apparently chaotic’ motions were observed by Andersen *et al.* (2005) for two-dimensional plates falling in a 3-D tank (the length of the plates being close to the dimensions of the tank, with the other dimensions of the plates, width and thickness, being much smaller).

In figures 1 and 2, the representation at constant time intervals of the position and orientation of the body along its path indicates that for given $ρ_c/ρ_f$ or Ar, shorter cylinders fall more rapidly than longer ones. We can also observe that for given $ρ_c/ρ_f$ or Ar, the elongation ratio of the cylinders significantly impacts their dynamics and kinematics. For the fluttering motions displayed by aluminium cylinders, as L increases (d being kept constant here), the displacement of the centre of gravity of the cylinder presents an increasing oscillation amplitude, of the order of L. For titanium cylinders, the amplitudes of oscillation in translation and rotation depend on $\xi$ in a complex way, as the regime of ‘apparently chaotic’ motions emerges (for intermediate elongation ratios $\xi = 5$ and 7). A strong link between the instantaneous orientation of the body and its vertical velocity is also clearly visible, for instance for titanium cylinders with $\xi = 10$, where faster vertical displacements are related to higher inclination angles of the body revolution axis relative to the horizontal direction. After reaching a maximum inclination angle, the cylinder plunges more rapidly into the cell than during the uprighting phase where its revolution axis recovers a horizontal alignment.
D. Letessier, V. Roig and P. Ern

<figure><img src="image_3.png" /><figcaption>Figure 3. Types of paths observed in the map of parameters Ar₃-D and ξ. Colours correspond to different values of the density ratio: red for ρ_c/ρ_f = 4.5, blue for ρ_c/ρ_f = 2.7, green for ρ_c/ρ_f = 1.4, and violet for ρ_c/ρ_f = 1.12 and 1.16. The plastic cylinders used by Toupoint (2018) are in fact marginally heavier than those of the present study, but are presented with the same colour code for simplicity. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

The different types of paths observed are summarized in figure 3 as a function of the Archimedes number Ar₃-D and the elongation ratio of the cylinders, ξ. Symbols in figure 3 indicate the type of path, and colours indicate the density ratio. The same colour code will be used throughout the paper: red for titanium cylinders (ρ_c/ρ_f = 4.5), blue for aluminium ones (ρ_c/ρ_f = 2.7), green for POMC ones (ρ_c/ρ_f = 1.4), and violet for plastic cylinders obtained by 3-D printing (ρ_c/ρ_f = 1.12 and 1.16).

## 2.2. The mean fall velocity

In this subsection, we investigate the dependence on the governing parameters of the mean fall velocity of the cylinders, denoted $\overline{u_v}e_v$, where $e_v$ is the vertical ascending unit vector (see also figure 9). Figure 4 presents the evolution of $\overline{u_v}$ as a function of the elongation ratio of the cylinder, for the different density ratios investigated (materialized with different colours as specified previously). The present set of experiments is presented in comparison to that of Toupoint (2018), which covers a complementary range of parameters. The slightly different gap sizes, fluid viscosities, cylinder geometries (which are more precise in the present study due to a different manufacturing process) and densities may explain the slightly different falling velocities recorded in the two studies. Despite the variability observed from one run to another, the strong effects on the mean fall velocity of both density ratio and elongation ratio are conspicuous, in line with the observations by D'Angelo et al. (2017) for different ranges of parameters.

In order to obtain a predictive model for the mean fall velocity, the classical approach consists in considering the mean balance in the vertical direction between buoyancy and drag forces. This balance can be written using the inertial drag expression that introduces a drag coefficient C_d and a reference frontal area S = Ld, as

$$ (\rho_c - \rho_f)g\pi L \frac{d^2}{4} + \frac{1}{2} \rho_f C_d L d |\overline{u_v}| \overline{u_v} = 0. \quad (2.1) $$
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_4.png" /><figcaption>Figure 4. Mean vertical velocities of the cylinders as functions of their elongation ratio $\xi$ for various density ratios. Dashed lines correspond to (3.15) and (3.14). Filled symbols indicate present study; square empty symbols indicate data from Toupoint (2018).</figcaption></figure>

This equation provides an expression for the drag coefficient

$$C_d = \frac{\pi |\rho_c - \rho_f|}{2 \rho_f} \frac{gd}{\bar{u}_v^2} = \frac{\pi}{2} \left( \frac{u_g}{\bar{u}_v} \right)^2, \qquad (2.2)$$

and brings forth the gravitational velocity

$$u_g = -\operatorname{sgn}(\rho_c - \rho_f) \left( \frac{|\rho_c - \rho_f|}{\rho_f} g d \right)^{1/2}. \qquad (2.3)$$

For a constant $C_d$, as observed, for instance, in the case of cylinders falling freely in a large tank at comparable Reynolds numbers, the gravitational velocity scale yields a good estimate of the mean fall velocity (Toupoint et al. 2019). The values of the drag coefficient obtained in the present configuration are shown in figure 5(a) as a function of the Reynolds number $Re_v = |\bar{u}_v| d/v_f$, for various density ratios and elongation ratios. The strong dependence of $C_d$ on the different parameters $\rho_c/\rho_f$, $\xi$ and $Re_v$ is clearly visible. As a consequence, the gravitational velocity $u_g$ fails to scale the mean fall velocity $\bar{u}_v$ in the present configuration. This can also be seen in figure 5(b), where the normalization by $u_g$ fails to capture the dependences of $\bar{u}_v$ on both $\xi$ and $\rho_c/\rho_f$. The dependence on the elongation ratio $\xi = L/d$ is not surprising, as $C_d$ can be expected to depend on the body shape, this dependence being also related to the choice for the reference frontal area $S$. It is worth pointing out here that we tried different choices for the reference area, such as the average of the area projected by the cylinder on a horizontal plane during a period (which is a function of the maximal inclination angle attained by the cylinder), suggested by Chow & Adams (2011), without succeeding in removing the dependence on $\xi$. The dependence on the density ratio $\rho_c/\rho_f$, especially for cylinders exhibiting the same type of path, also deserves further attention, as it is not captured by the dependence of $C_d$ on the Reynolds number (figure 5a) or on the Archimedes number.

When using the classical balance between buoyancy and drag forces in the present configuration, we are therefore left with the result $C_d = f(Re_v, \rho_c/\rho_f, \xi)$ or
<figure><img src="image_5.png" /><figcaption>Figure 5. (a) Drag coefficient $C_d$ as a function of the Reynolds number $Re_v$ of the cylinder, for various elongation ratios $\xi$; same colour code for the density ratios as in (b). (b) Mean vertical velocity normalized with the gravitational velocity scale $u_g$. Filled symbols indicate present study; empty symbols indicate data from Toupoint (2018).</figcaption></figure>

$C_d = f(Ar_{3-D}, \rho_c/\rho_f, \xi)$, which is poorly satisfying. At this stage, a possibility is to propose correlations for $C_d$, as those proposed for cylinders in free fall in 3-D environments by Clift et al. (1978), Chow & Adams (2011), and references therein. These correlations are useful to predict the mean fall velocity of the cylinders in various practical situations, but lack of physical support regarding, in particular, the dependence of $C_d$ on the density ratio. To tackle this question, an alternative is to improve the modelling of the loads governing the cylinder motion, poorly represented by (2.1). This is the approach followed in this paper, by considering force contributions nonlinearly coupling the degrees of freedom of the body, in particular fluctuating ones. Therefore, we now turn our attention to the quantities characterizing the oscillatory motion of the body.

## 2.3. Characteristics of the oscillatory motion

We focus here on the oscillatory motion of the cylinders, called fluttering, occurring over a large range of parameters for both aluminium and titanium cylinders, as shown in figure 3. In this regime, the three degrees of freedom of the body – the vertical and horizontal velocities, denoted respectively $u_v$ and $u_h$, as well as the inclination angle of the axis of revolution of the cylinder relative to the horizontal, $\theta$ (see figure 9) – evolve in space and time in a periodic manner. Typical temporal evolutions of these quantities are displayed in figure 6 for two cylinders having the same density ratio $\rho_c/\rho_f = 2.7$, but different elongation ratios $\xi = 3$ (figure 6a) and $\xi = 10$ (figure 6b). As expected, the vertical velocity (green line) presents a mean component, and evolves at twice the frequency $f_\theta$ of oscillation characterizing both the horizontal velocity component (yellow line) and the inclination angle (dark red line). The comparison between the two cylinders also reveals a ratio close to a factor 2 between their oscillation frequencies $f_\theta$, the quantities for the cylinder with the smallest elongation ratio presenting oscillations with the shortest time period. The amplitudes displayed by the orientations of both cylinders are, however, similar, with values of $\theta$ exceeding 30°. In both cases, it is also conspicuous that the horizontal velocity oscillates with an amplitude larger than that of the vertical velocity, which is nonetheless large. The most important observation, however, is that the amplitude of oscillation of the horizontal velocity even exceeds the mean vertical velocity for the
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_6.png" /><figcaption>Figure 6. Temporal evolutions of the vertical velocity $u_v$ (green), the horizontal velocity $u_h$ (yellow), and the inclination angle $\theta$ (dark red dashed line) for two aluminium cylinders falling freely in water at rest: (a) $\xi = 3$, $Ar_{3-D} = 193$; (b) $\xi = 10$, $Ar_{3-D} = 352$.</figcaption></figure>

highest elongation ratio, despite the body falling with a significant velocity. This property is at variance with most cases investigated in the literature, where the mean vertical velocity is usually an order of magnitude larger than the fluctuation amplitudes, even for oscillatory displacements of the order of the body size (Ern et al. 2012).

The strong amplitudes involved in the fluttering motion observed in the present configuration can be expected to induce nonlinear effects, which will be investigated in the next section. The nonlinearity in the body path is already visible in the time signals of figure 6, where the signals depart from purely sinusoidal functions. This is confirmed by the spectral analysis of the signals, which provides a strong predominant frequency $f_\theta$ for $u_h$ and $\theta$, and $2f_\theta$ for $u_v$, along with higher-order frequencies of weaker amplitude, in particular harmonics induced by the nonlinear coupling between translational and rotational motions. To proceed with the analysis, we focus our attention on the main-frequency component of the fluttering motion defined as

$$ \left. \begin{aligned} u_h(t) &= \tilde{u}_h \sin(\omega t + \phi_h), \\ u_v(t) &= \overline{u_v} + \tilde{u}_v \sin(2\omega t + \phi_v), \\ \theta(t) &= \tilde{\theta} \sin(\omega t + \phi_\theta), \end{aligned} \right\} \qquad (2.4) $$

where $\omega = 2\pi f_\theta$. Quantities $\tilde{u}_h$, $\tilde{u}_v$ and $\tilde{\theta}$ denote, respectively, the amplitudes of the horizontal and vertical components of the velocity and of the inclination, while $\phi_h$, $\phi_v$ and $\phi_\theta$ are their phases. Figure 7(a) shows the dependence of the fluttering frequency $f_\theta$ on the elongation ratio $\xi$. Shorter cylinders oscillate at higher frequencies, as also observed for cylinders in 3-D environments (Marchildon, Clamen & Gauvin 1964; Chow & Adams 2011; Toupoint et al. 2019, to cite just a few). Important frequency differences
<figure><img src="image_7.png" /><figcaption>Figure 7. (a) Fluttering frequency $f_θ$ and (b) amplitude of inclination $\tilde{\theta}$ (in degrees) as a function of the elongation ratio of the cylinder, for two density ratios.</figcaption></figure>

are observed: cylinders with $\xi = 20$ oscillate at a frequency four times lower than that displayed by cylinders with $\xi = 3$. Conversely, for a given elongation ratio, the frequencies displayed by cylinders made of different materials are very close.

Figure 7(b) shows the dependence on the elongation ratio of the oscillation amplitude of the inclination, $\tilde{\theta}$ (in degrees). A non-monotonic evolution is observed as $\xi$ increases. A strong effect of the density ratio is also visible for given $\xi$, with the largest amplitude for the largest density ratio. In the range of parameters investigated here, the amplitudes of the inclination angle $\tilde{\theta}$ evolve in the fluttering regime over a large range of values, spanning from 10° to 70°. They exceed 90° for the cases marked at this value, during the tumbling phases observed for titanium cylinders having $\xi = 5$ or 7.

The amplitudes of the velocity components along the vertical and horizontal directions are presented in figure 8(a), normalized with the mean fall velocity $|\bar{u}_v|$, as a function of $Ar_{3-D}$ for various $\xi$. As noticed in the two examples of figure 6, the amplitudes of the vertical fluctuations are smaller than those recorded for the horizontal component over the whole parameter range. Furthermore, the amplitudes of the horizontal velocity compare to the mean vertical velocity values, or even exceed them for elongation ratios $\xi > 7$. A similar result is shown in figure 8(b) for the rotational velocity of the cylinder, with amplitude values $L\tilde{\theta}$ of the order of magnitude of $|\bar{u}_v|$. It can also be seen that the Archimedes number $Ar_{3-D}$ is not sufficient to organize the data consistently for the two density ratios investigated. For a given elongation ratio, both figures 8(a) and 8(b) show a dependence on the density ratio as well as on $Ar_{3-D}$. However, this dependence is more pronounced for the translational velocity fluctuations (figure 8a), where higher amplitudes are observed at the highest density ratio corresponding also to higher $Ar_{3-D}$ values. In contrast, the amplitude values for the angular velocity fluctuations remain relatively equivalent. In complement to the amplitudes, the phase difference between the body oscillations in horizontal velocity and inclination, $\phi_h - \phi_\theta$, is displayed in figure 8(c). For the two values of $\rho_c/\rho_f$ investigated, a similar evolution is observed as a function of $\xi$, the phase lag increasing slightly as the body elongates.

As a conclusion, the detailed analysis of the characteristics of the oscillatory motion reveals strong amplitudes of fluctuation for both translational and rotational velocity components, comparable in magnitude to the mean fall velocity. This indicates that a good prediction of the mean fall velocity of the cylinder can be achieved in the present case only
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_8.png" /><figcaption>Figure 8. Oscillation amplitudes of (a) the velocity components along the vertical (filled symbols) and horizontal (open symbols) directions, and (b) the angular velocity, normalized with the mean vertical velocity. (c) Phase difference $\phi_h - \phi_\theta$. Aluminium cylinders in blue; titanium cylinders in red.</figcaption></figure>

by considering altogether the different degrees of freedom, and by taking into account the coupling between translational and rotational velocity components. The next section is therefore devoted to build a force balance able to model the behaviour of the cylinders over the whole range of parameters investigated.

### 3. A model for the paths of cylinders falling in fluid at rest

We have seen that the behaviour of freely falling cylinders depends strongly, and in a complex manner, on the governing parameters, and that the different degrees of freedom of the cylinder, in both translation and rotation, present comparable levels of energy. On the basis of physical arguments, here we build up a model that, as large as possible in the ranges of $Ar_{3-D}$, $\rho_c/\rho_f$ and $\xi$ investigated experimentally, (i) proposes a mean force balance on the cylinder that governs its mean fall velocity, and (ii) is at instantaneous equilibrium when introducing the experimental temporal evolutions of the body's degrees of freedom. To do so, we will proceed by steps, and first determine the different forces that are expected to enter in the force balance.

A difficulty here is that the generalized Kelvin–Kirchhoff equations (Lamb 1932; Howe 1995; Mougin & Magnaudet 2002) that express the conservations of momentum and of moment of momentum, cannot be applied, as they are derived for a body free to move in an unbounded incompressible fluid at rest. However, on the basis of these equations, the following force balance can be considered to govern the motion of the body centre of mass:

$$F_I = F_B + F_\omega. \qquad (3.1)$$
D. Letessier, V. Roig and P. Ern

<figure><img src="image_9.png" /><figcaption>Figure 9. Scheme of the freely moving cylinder, indicating the different notations used in the model.</figcaption></figure>

On the left-hand side, $F_I$ stands for the inertial forces, involving both proper and added inertia contributions. On the right-hand side, $F_B$ is the force due to buoyancy, and $F_\omega$ corresponds to the force associated with the production of vorticity at the body surface, also accounting for its interaction with the confining walls. The former is known, and reads $F_B = (\rho_c - \rho_f)Vg$, where $V$ is the body volume, and $g$ is the gravitational acceleration, while $F_\omega$ is unknown.

In the formalism of the Kelvin–Kirchhoff equations, the inertial forces on a body of mass $m = \rho_c V$ read

$$F_I = (m\mathbb{I} + \mathbb{A}) \frac{du_c}{dt} + \mathbf{\Omega} \times ((m\mathbb{I} + \mathbb{A})u_c), \quad (3.2)$$

where $u_c$ is the velocity of the centre of mass, and $\mathbf{\Omega}$ is the angular velocity of the cylinder. The notation $\mathbb{I}$ stands for the identity tensor, and $\mathbb{A}$ for the added-mass tensor. For convenience, we also introduce the notations $F_{IT}$ and $F_{IR}$, standing for the inertia components associated with, respectively, the translation and the rotation of the body, so that the right-hand side of (3.2) reads $F_{IT} + F_{IR}$. The equations are expressed in a coordinate system having a fixed origin $O$ in the laboratory and axes rotating with the body. As shown in figure 9, the longitudinal direction is parallel to the revolution axis of the cylinder (unit vector $e_l$), and the transverse direction is perpendicular to it in the plane of the cell (unit vector $e_t$). The system of axes is completed with the direction perpendicular to the plane of cell, corresponding to unit vector $e_z$. In this system of axes, the added-mass tensor is diagonal, with added-mass coefficient $A$ along $e_l$, and $B$ along $e_t$. These coefficients depend on the body geometry, and are known for only several simple geometries. They are unknown in the present situation involving confining walls. The laboratory frame is given by $(O, e_h, e_v, e_z)$, where $e_h$ denotes the horizontal unit vector, while $e_v$ is the ascending vertical unit vector. The angle $\theta$ characterizes the inclination of the cylinder relative to the horizontal, as seen before, and the drift angle $\alpha$ measures the inclination of the velocity vector $u_c$ relative to the body revolution axis $e_l$. Note that in the present configuration, the velocity component in the direction $e_z$ is $u_z = 0$. Moreover, the angular velocity is given by $\mathbf{\Omega} = (0, 0, \dot{\theta})$. Finally, we define the velocity components as $u_c = u_h e_h + u_v e_v = u_l e_l + u_t e_t$. Equations (3.1) and (3.2) can therefore be recast in the form

$$(\rho_c V + A)\dot{u}_l - (\rho_c V + B)\dot{\theta}u_t = F_\omega^l - \Delta\rho Vg \sin\theta, \quad (3.3)$$

$$(\rho_c V + B)\dot{u}_t + (\rho_c V + A)\dot{\theta}u_l = F_\omega^t - \Delta\rho Vg \cos\theta, \quad (3.4)$$

where $F_\omega = F_\omega^l e_l + F_\omega^t e_t$ and $\Delta\rho = \rho_c - \rho_f$. For cylinders of finite length in an unconfined fluid, Loewenberg (1993) determined numerically the values of the added-mass coefficients in the longitudinal and transverse directions, $A$ and $B$,
## Fluttering motion of a confined freely falling cylinder

respectively, for several values of the elongation ratio. Pierson (2023) proposed the following fits for these data:

$$A_u = \rho_f \pi \frac{d^3}{4} \left( 0.655 - \frac{0.141}{1 + \xi^{1.17}} \right) \quad \text{and} \quad B_u = \rho_f \pi \frac{d^3}{4} \xi \left( 1 - \frac{0.828}{1 + \xi^{1.12}} \right). \quad (3.5a,b)$$

In the present configuration, as no data are available for the added-mass coefficients, we will assume an identical evolution with $\xi$. For infinite cylinders, Blevins & Plunkett (1980) showed numerically that the confinement modifies the values of the added-mass coefficients relative to an unconfined situation by a factor depending on the degree of confinement. They showed that the confinement increases the values of B (transverse direction), by a factor of 1.83 for a confinement equal to that of our configuration. No estimation was provided for the longitudinal component, but the multiplying factor can be expected to be weaker. We will therefore use here $B = 1.83B_u$ and $A = A_u$. Nonetheless, for the fluttering cases observed here for aluminium and titanium cylinders, it turns out that the term associated with the proper mass of the body $\rho_c V$ predominates over the added-mass contributions obtained in this manner. Knowing the exact values of the coefficients A and B is therefore not needed for the analysis and conclusions that are developed in the paper. We can now see that entering the experimental measurements of the body kinematics into (3.3)–(3.4) provides an instantaneous measure of the longitudinal and transverse components $F_\omega^l$ and $F_\omega^t$ of the vortical force.

### 3.1. Preliminary considerations

#### 3.1.1. Mean force balance along the transverse direction

Before proceeding with the modelling of $F_\omega$, we present in this section the information for this force that is provided by the measurements of the body kinematics, as these allow us to determine both $F_I$ and $F_B$. As the first goal of the model is to achieve a better understanding of the forces controlling the mean fall velocity of the cylinder, we will start by investigating the mean force balance along the transverse direction. Averaged over a time period of the fluttering motion, the conservation of momentum in the direction $e_t$ given by (3.4) reads

$$\langle F_\omega^l \rangle = (\rho_c V + A) \langle \dot{\theta} u_l \rangle + \Delta \rho V g \langle \cos \theta \rangle, \quad (3.6)$$

where, for a quantity Q, the notation $\langle Q \rangle$ corresponds to its average over a time period. Note that here, the inertia term in translation has a zero mean contribution, and is therefore omitted. We now investigate the different terms composing equation (3.6) for the fluttering cases provided by aluminium ($\rho_c/\rho_f = 2.7$) and titanium ($\rho_c/\rho_f = 4.5$) cylinders. For the two density ratios, figure 10 presents the mean contributions measured for the inertial force, the buoyancy force and the vortical force, when the elongation ratio of the cylinder varies. The plotted values correspond to the mean contributions of the considered forces normalized with the magnitude of the buoyancy force $\Delta\rho Vg$. The values for the vortical force are obtained by summing the two other contributions (inertia and buoyancy), according to (3.6). Several remarkable results are conspicuous. First, for all the parameters ($\xi, \rho_c/\rho_f$) investigated, we can see that the mean value of the measured vortical force in the transverse direction is nearly constant and very close to $\Delta\rho Vg$. Interestingly, this extends to unsteady paths a result holding for rectilinear motions, for which the transverse direction corresponds to the vertical direction, and the vortical force reduces to the drag force (see (2.1) or (3.11)). A second observation is that this
<figure><img src="image_10.png" /><figcaption>Figure 10. Measured mean contributions to the force balance in the transverse direction as a function of the elongation ratio of the cylinders, for two density ratios. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

result is achieved because the mean contribution of the inertial force associated with the body rotation has a contribution of opposite sign and complementary in magnitude to that of the buoyancy force. In particular, when the mean buoyancy contribution along the transverse direction decreases due to higher fluttering inclinations $\theta$, the coupling of the rotation rate $\dot{\theta}$ with the longitudinal velocity component $u_l$ provides a mean contribution equivalent to this decrease. The changes in the control parameters ($\xi$, $\rho_c/\rho_f$) also show that this compensation happens for a large range of magnitude of these forces (decrease in buoyancy reaching 50% for titanium cylinders). Finally, these results reveal the salient role and the important magnitudes that can be achieved by the inertia nonlinear term associated with the body rotation. When the amplitudes of oscillation of the body velocity in translation and rotation are weak relative to the mean fall velocity, this term is usually weak. However, when amplitudes of oscillation comparable to the mean fall velocity are at play and the fluctuating degrees of freedom are out of phase, a mean component arises at leading order in the mean force balance. As shown for the fluttering cases investigated here, the oscillatory motion then impacts the mean force balance, and classical estimations using the gravitational velocity scale fail to predict the mean fall velocity. Furthermore, this force balance enlightens how the density ratio enters additionally into play through this inertia effect. As will be shown in the following, the consideration of this inertia effect avoids assuming a dependence of $C_d$ on $\rho_c/\rho_f$ in the modelling of the drag term, and more generally, of the vortical force $F_\omega$. Finally, it is worth pointing out that the inertial force coupling the translational and rotational degrees of freedom of the body has two components: one associated with the proper inertia of the cylinder, and one with the added inertia of the surrounding fluid. Depending on the governing control parameters, the relative weights of these two terms can be different. Here, the proper inertia term is predominant.

### 3.1.2. Scaling law for the fluttering frequency

The analysis of the mean force balance along the transverse direction of the cylinder has shown that the translational and rotational velocity components couple, over a time period, in a way such that

$$ \langle F_{\omega}^{t} \rangle = (\rho_c V + A) \langle \dot{\theta} u_l \rangle + \Delta\rho V g \langle \cos \theta \rangle \approx \Delta\rho V g, \quad (3.7) $$
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_11.png" /><figcaption>Figure 11. Fluttering frequency normalized with $\sqrt{(\Delta\rho/\rho_c)(g/L)}$ for various elongation ratios and density ratios. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

in the range of values of $\xi$ and $\rho_c/\rho_f$ investigated. We have also seen that added mass can be neglected here relative to the proper mass of the cylinder. Equation (3.7) can therefore be written in dimensionless form, using $L$ as the characteristic length scale associated with the longitudinal velocity and $1/f_\theta$ as the time scale, as

$$ \frac{\rho_c L f_{\theta}^2}{\Delta\rho g} \langle \dot{u}^* u_l^* \rangle + \langle \cos \theta^* \rangle \approx 1, \quad (3.8) $$

with dimensionless quantities indicated with an asterisk. From this balance, the dimensionless quantity $\rho_c L f_{\theta}^2 / (\Delta\rho g)$ emerges, and allows us to propose the following scale for the fluttering frequency:

$$ f_{\theta} \sim \sqrt{\frac{\Delta\rho}{\rho_c} \frac{g}{L}}. \quad (3.9) $$

Figure 11 shows the frequency of oscillation measured experimentally, normalized with the prediction. We can observe that the data points gather consistently for the two density ratios about a mean value approximately 0.13. This result confirms that for all the cases investigated, the degrees of freedom, as well as their rate of change, couple in a very consistent manner over a time period to maintain $\langle F_{\omega}^t \rangle$ close to $\Delta\rho V_g$. Interestingly, the time scale of (3.9) was proposed and observed experimentally by Marchildon et al. (1964) and Chow & Adams (2011) for cylinders freely falling in an unbounded medium for Reynolds numbers in between 200 and 6000, and by D’Angelo et al. (2017) for confined cylinders in a complementary parameter range. Marchildon et al. (1964) and Chow & Adams (2011) introduce this frequency scaling by considering a torque equation on the cylinder that balances the time variation of the moment of momentum with a hydrodynamical torque related to a difference between the centre of gravity of the body and the centre of pressure, that depends linearly on the angle $\theta$ (denoted $\alpha$ in their works). We claim, however, that this torque balance is not satisfactory (except for high density ratios) due to the omission of the added-mass torque, which has a strong contribution to the torque balance (Ern et al. 2012). Their experimental results correspond
to a pre-factor 0.126, extremely close to the value 0.13 obtained in the present work. This leads us to consider that the results and conclusions obtained here from the mean transverse force balance (3.7) are generic and may also apply to non-confined cylinders exhibiting high amplitudes of fluttering motion. Depending on the density ratio and on the body geometry, including added mass could be necessary, through the more general normalization $f_\theta \sim \sqrt{(\Delta\rho/(\rho_c + A/V))(g/L)}$.

## 3.2. The different contributions to the vortical force

We now aim to obtain a model for the vortical force $F_\omega$ that will then be compared to the experimental determination for this force. The vortical force acting on the body can be expected to be associated with classical drag and lift forces, respectively, $F_D$ and $F_\Gamma$. Usually, drag and lift definitions are related to instantaneous quantities. In the present situation, history effects can also be expected to play a significant role. They will be taken into account by a force denoted $F_H$. To model $F_\omega$, we therefore consider three contributions:

$$F_\omega = F_D + F_\Gamma + F_H. \qquad (3.10)$$

Each term on the right-hand side will now be elaborated separately. To determine the drag force, rectilinear motions will first be analysed. The expression will then be extended to fluttering motions, using a quasi-steady approximation, by taking into account the instantaneous orientation of the velocity vector relative to the body revolution axis. We will then consider the lift contribution, and following the literature, we will suggest an expression in the form of a translational and a rotational lift, based on the circulation $\Gamma$. We will finally propose an integro-differential expression for the history force that enables (3.10) to match the experimental determination of $F_\omega$.

### 3.2.1. Drag force $F_D$ for cylinders following a rectilinear path

For bodies displaying a rectilinear motion, the drag force $F_D = F_D^l e_l + F_D^t e_t$ reduces to the transverse component, which balances the buoyancy force, so that

$$F_D^t = (\rho_c - \rho_f)g\pi L \frac{d^2}{4}. \qquad (3.11)$$

As the cylinders fall with their revolution axis aligned with the horizontal direction, we have $u_l = 0$, $\theta = 0^\circ$, $\alpha = 90^\circ$ and $u_v = u_t$, which will be denoted $u_0$ in the following model. The drag coefficient associated with the rectilinear motion is denoted $C_d^{\alpha=90^\circ}$, and the drag force then reads

$$F_D^t = -\frac{1}{2}\rho_f S C_d^{\alpha=90^\circ} |u_0| u_0, \qquad (3.12)$$

with $S = dL$. For consistency with what follows, we introduce the Reynolds number

$$Re = \frac{u_c d}{\nu_f}, \qquad (3.13)$$

based on the norm of $u_c$, $u_c = \sqrt{u_l^2 + u_t^2}$. For a rectilinear motion, we have $Re = Re_v$. For cylinders with $\rho_c/\rho_f = 1.12$ and 1.16 displaying a rectilinear path (figure 3), we seek a
## Fluttering motion of a confined freely falling cylinder

dependence of $C_d^{\alpha=90^{\circ}}$ on $\xi$ and $Re$ (violet data points of figure 5a) and obtain

$$C_d^{\alpha=90^{\circ}} = \frac{\kappa_1^{\alpha=90^{\circ}}(\xi)}{Re}, \quad \text{with } \kappa_1^{\alpha=90^{\circ}}(\xi) = 73 + 6.3\xi. \qquad (3.14)$$

Note that this expression quantifies the drag experienced by cylinders in rectilinear motion in this specific configuration, accounting both for their wake and for the presence of the confinement. Combining (3.11) and (3.12) then gives

$$u_0 = \operatorname{sgn}(u_g) \frac{\pi}{2} u_g^2 \frac{d}{\nu_f} \frac{1}{\kappa_1^{\alpha=90^{\circ}}(\xi)}, \qquad (3.15)$$

where $u_g$ is the gravitational velocity scale given by (2.3). At this stage, it is interesting to point out that the velocity scale $u_0$ presents a dependence on $\rho_c/\rho_f$ different from that of $u_g$. Note also that the dependence on $Re$ has a form similar to the expression $C_d^{\alpha=90^{\circ}} \propto Re^{-0.78}$ proposed by Clift et al. (1978) based on the data from Pruppacher, Clair & Hamielec (1970) for an unconfined long fixed cylinder held in an incident flow up to Reynolds number $Re = 400$. In the thin-gap cell configuration investigated here, an additional dependence on $d/w$ in (3.14) is expected to exist, but is out of the scope of the present paper (this parameter is not varied). The comparison of $u_0$ with the experimental results is provided in figure 4, where the evolution for $u_0$ is drawn with dashed lines for each density ratio, including those for which fluttering motions are observed. For a more detailed comparison, figure 12 also shows for all the cases the mean transverse velocity measured experimentally, normalized with $u_0$. As the model was adjusted for the rectilinear motions, good agreement is obtained for these cases (violet data points corresponding to the density ratios $\rho_c/\rho_f = 1.12$ and $1.16$), the model accurately predicting the falling speed within 6%. When considering also cylinders with $\rho_c/\rho_f = 1.4$ in rectilinear motion (green data points), the error increases to 9%. However, microscope examination of these cylinders revealed shape imperfections, notably the presence of filaments at their extremities, having lengths comparable to the diameter of the cylinder. Such defects may be expected to slightly increase the drag experienced by the cylinders. We therefore do not consider their deviation from the model to consist in a significant drawback. For the fluttering cases, as seen from the analysis of § 3.1.1, the balance provided by (3.11) is incomplete, and the velocity scale $u_0$, as well as $u_g$, strongly overestimates the mean transverse component of the velocity. However, on the basis of the modelling for rectilinear motions, we can now proceed with trying to improve the drag modelling for fluttering motions.

### 3.2.2. Drag force $F_D$ for both rectilinear and periodic motions

The simplest way to extend the previous drag model to fluttering motions is to adopt a quasi-steady approach that takes into account at each time the instantaneous orientation of the velocity vector relative to $e_l$, characterized by the angle $\alpha$ (see figure 9). This approach was, in particular, used by Andersen et al. (2005) for the modelling of the behaviour of freely falling two-dimensional plates, exhibiting paths comparable with the ones investigated here. The quasi-steady approximation makes possible to use the results for a body fixed in a uniform flow of velocity $u_c$ and oriented at angle $\alpha$ relative to the flow direction, when these results are known. However, for finite-length cylinders, even in this simpler situation, there is no general expression for the drag force. A possibility is then to adopt a result valid in the Stokes regime, based on the linearity of the equations,
D. Letessier, V. Roig and P. Ern

<figure><img src="image_12.png" /><figcaption>Figure 12. Mean transverse velocity component normalized with $u_0$. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

that considers that the cylinder receives in the proportion of cos α the force experienced when placed parallel to the incident flow (α = 0°) and of sin α the force that it experiences when perpendicular to the flow (α = 90°). We assume that this property still holds for the inertial regime of fluttering motions considered here. As, by definition, the drag force is aligned with the velocity vector $u_c$, we then get

$$C_d = C_d^{\alpha=0^{\circ}} \cos^2 \alpha + C_d^{\alpha=90^{\circ}} \sin^2 \alpha. \qquad (3.16)$$

At this stage, the choice of Andersen et al. (2005) was to adjust the values of $C_d^{\alpha=0^{\circ}}$ and $C_d^{\alpha=90^{\circ}}$ in order to match their experimental results for four different cases of freely falling plates. As we have a larger data set, we seek here an explicit dependence on $\xi$ and $Re$ for the drag coefficients $C_d^{\alpha=0^{\circ}}$ and $C_d^{\alpha=90^{\circ}}$. The latter has already been obtained in (3.14). As $C_d^{\alpha=0^{\circ}}$ is unknown, we assume $C_d^{\alpha=0^{\circ}} = \frac{1}{2}C_d^{\alpha=90^{\circ}}$, the theoretical relation derived in the slender-body approximation (for $d \ll L$) in the viscous regime and for an unconfined fluid (Batchelor 1970). Note that this relation was shown to still provide satisfactory results for a fixed cylinder with $\xi = 3$ and Reynolds numbers between 25 and 100 in the numerical investigation by Pierson et al. (2019). This series of considerations leads us to propose the following expression for the quasi-steady drag force experienced by the freely falling cylinder:

$$F_D = -\frac{1}{2} \rho_f C_d(\xi, Re, \alpha) S u_c(u_l e_l + u_t e_t), \qquad (3.17)$$

with

$$C_d(\xi, Re, \alpha) = C_d^{\alpha=90^{\circ}} \left(\frac{1}{2} \cos^2 \alpha + \sin^2 \alpha\right), \qquad (3.18)$$

where $S = dL$, $u_c = \sqrt{u_l^2 + u_t^2}$ and $C_d^{\alpha=90^{\circ}} = \kappa_1^{\alpha=90^{\circ}}/Re$ as determined previously for rectilinear motions in the present configuration (3.14). Though strong approximations were used, the strength of this simple expression is to take into account the instantaneous drift angle $\alpha$ in a way that is consistent with the rectilinear motion (when $\alpha = 0^{\circ}$) and with cylinders fixed in an incoming flow (when $\dot{\alpha} = 0^{\circ} \text{ s}^{-1}$).
### 3.2.3. Lift force $F_Γ$ due to the circulation $Γ$

To complement the quasi-steady drag force, we follow Pesavento & Wang (2004) and Andersen et al. (2005), and define a force $F_Γ$ proportional to the circulation $Γ$, as

$$F_Γ = \rho_f Γ(-u_t e_l + u_l e_t). \qquad (3.19)$$

As in these studies, we also assume two contributions to the circulation,

$$Γ = -C_T L d \frac{u_t u_t}{u_c} + \frac{1}{2} C_R d L^2 \dot{θ}, \qquad (3.20)$$

so that the force $F_Γ$ has two components, $F_Γ = F_{Γ_T} + F_{Γ_R}$. The former component, proportional to $C_T$ and $u_t u_t = u_c \sin(2\alpha)$, can be assimilated to a translation-induced lift. The latter component, proportional to $C_R$ and the angular velocity $\dot{θ}$, can be understood as a rotation-induced lift. Andersen et al. (2005) showed that the rotation-induced lift has an essential contribution to the force balance, and found values for $C_T$ and $C_R$ that allowed them to close the force balance for each path considered. However, they were not able to quantitatively evaluate the dependence of the coefficients $C_T$ and $C_R$ on the control parameters due to the limited number of experimental cases available (four). As done previously for the drag force, we here seek evolutions of the coefficients $C_T(Re)$ and $C_R(Re)$ with the Reynolds number $Re$ of the cylinders. As will be discussed later in detail (§ 3.3.1), such expressions can be found here satisfactorily to close the mean force balance along the transverse direction, obtained by averaging (3.4) over a period (there is no mean component for the longitudinal force balance). However, we were not able to achieve the instantaneous equilibrium of the forces along the longitudinal and transverse directions using the precedent expressions for the drag and lift terms. While, of course, the models for these terms could be revisited and improved, we attribute their failure to their dependence on the instantaneous kinematics of the cylinder. We therefore seek a correction in the form of a history force.

### 3.2.4. History force $F_H$

Here, we seek to incorporate in the model a force that will be able to guarantee the instantaneous force balance in both the longitudinal and transverse directions. With this force, we try to correct the quasi-steady character of the models proposed for $F_D$ and $F_Γ$. As the flow at a given time is clearly dependent on the flow at preceding times, this is also the case for the hydrodynamical loads experienced by the body. This dependence of the loads on the preceding times cannot be accounted for by the instantaneous values of the body kinematics. The simplest idea is therefore to take into account a rate of change for these quantities. However, the choice of the relevant quantities and the formulation of this history effect is not trivial. An alternative explored here is to consider the history of the temporal rate of change of the loads themselves. A careful analysis of the different force components, in particular their phase, led us to the following choice for $F_H$:

$$F_H = K_l St_\theta \int_{t-\tau_l}^{t} e^{-(t-t^*)/\tau_l} \frac{dF_\Gamma^t}{dt^*} dt^* e_l + K_l \int_{t-\tau_l}^{t} e^{-(t-t^*)/\tau_l} \frac{dF_D^t}{dt^*} dt^* e_l, \quad \text{with } St_\theta = \frac{df_\theta}{u_c}. \qquad (3.21)$$

The coefficients $K_l$ and $K_l$ allow us to adjust separately the amplitudes of the two components. While the expressions are different for the longitudinal and transverse directions, they have the same shape and the same meaning. Using an integro-differential
formulation inspired by the Basset-Boussinesq history force (Basset 1888) and by the force proposed by Lawrence & Weinbaum (1986) at low Reynolds numbers, the contributions of $F_H$ are composed with the product of a kernel with the temporal derivatives of the drag (for the transverse direction) and lift (for the longitudinal direction) forces. Different time scales are considered for the longitudinal ($\tau_l$) and transverse ($\tau_t$) components, though close values will be obtained when tuning the parameters of the model. The integration is performed over a characteristic time interval $\tau_l$ (and $\tau_t$, respectively) preceding the instantaneous time $t$, expected to be smaller than a period of oscillation ($\tau_l f_\theta < 1$ and $\tau_t f_\theta < 1$). The kernel is chosen, for simplicity, as an exponential decaying function of same characteristic time scale ($\tau_l$ or $\tau_t$). A similar kernel was used by Ern et al. (2009) to model the vortical torque experienced by freely falling disks. Note also that we introduce the Strouhal number $St_\theta$ to modulate the amplitude of the longitudinal component of $F_H$. The amplitude of this term is in fact strongly dependent on the elongation ratio of the body, shorter cylinders corresponding to larger Strouhal numbers. As the Strouhal number compares the rate of change of the flow surrounding the body ($f_\theta$) to the inertial time ($d/|\bar{u}_c|$), the quasi-steady approximation is more justified for lower values of $St_\theta$ (occurring for longer cylinders), and history effects are then also less important. This can be understood clearly by comparing the behaviour of long and short cylinders in figure 1 (or in figure 2). The vertical displacement of the body, plotted at equal time steps, is seen to scan very differently a period of oscillation in each case, due to the strong differences in velocity and oscillation frequency.

### 3.2.5. Values of the coefficients used in the model

The calibration of the eight coefficients ($C_T(Re)$, $C_R(Re)$, $K_l$, $K_t$, $\tau_l$, $\tau_t$ and two coefficients for $\kappa_1^{\alpha=90^{\circ}}(\xi)$) already introduced in (3.14)) involved in the model for the vortical force $F_\omega$ was performed by steps. Two of these parameters were first set from the investigation of the rectilinear motions, as seen in § 3.2.1. The consideration of the mean force balance in the transverse direction then allowed us to obtain suitable values for the coefficients $C_T(Re)$ and $C_R(Re)$. The optimization of these parameters was performed using a cost function that evaluates the average value over the entire data set of the sum of the different contributions entering in the force balance for each cylinder ($\rho_c/\rho_f$, $\xi$, $Ar_{3-D}$), expected to be zero for each case.

Assuming the laws for $C_T(Re)$ and $C_R(Re)$ obtained in this manner, we then look for the remaining parameters, $K_l$, $K_t$, $\tau_l$ and $\tau_t$, that close the instantaneous force balances in both the transverse and longitudinal directions. These were obtained by minimizing the sum of the contributions at each time for each experimental run, focusing on aluminium cylinders. The dimensionless coefficients $K_l$ and $K_t$ were initialized with values close to unity, and the initial values for the times $\tau_l$ and $\tau_t$ with values close to the period of oscillation $T$. After this first determination in consecutive steps, the obtained values were used as initialization values for the optimization process of the six coefficients let free. Values close, yet different, to those from the iterative procedure were found. These provide the following closure relations:

$$ \kappa_1^{\alpha=90^{\circ}}(\xi) = 73 + 6.3\xi, \qquad (3.22a) $$

$$ C_T(Re) = \frac{145}{Re} \quad \text{and} \quad C_R(Re) = \frac{66}{Re}, \qquad (3.23a,b) $$

$$ K_l = -14.43, \quad \tau_l = 0.53T, \qquad (3.24a,b) $$

$$ K_t = 2.63, \quad \tau_t = 0.63T. \qquad (3.25a,b) $$
## Fluttering motion of a confined freely falling cylinder

Similarly to the drag coefficient $C_d$, the coefficients involved in the lift force, $C_T$ and $C_R$, are inversely proportional to the Reynolds number $Re$ in the range of parameters investigated here. As concerns the history force, the time scales $\tau_l$ and $\tau_t$ correspond to approximately half a period of oscillation $T$, and have comparable magnitudes in both directions. Interestingly, the coefficient $K_l$ is negative, contributing along the longitudinal direction in a way opposite to the lift, the lift force providing a centripetal force contribution on the cylinder, and the history force a centrifugal one (see figure 16).

### 3.3. Discussion of the results

The comparison of the model with the experimental data is carried out in this subsection. We then discuss the relative weights of the different contributions involved in the vortical force $F_\omega$, and their evolution along the path of the moving body. For each force $F$, we define its components $F^l$ and $F^t$ along the longitudinal and transverse directions of the cylinders as $F = F^l e_l + F^t e_t$.

#### 3.3.1. Closure of the mean force balance along the transverse direction

The first goal of the force model elaborated for $F_\omega$ is to close the mean force balance in the transverse direction, investigated in § 3.1.1. Averaged over a time period, the conservation of momentum in the direction $e_t$ given by (3.4), along with the force decomposition of (3.10), reads

$$ \sum \langle F^t \rangle = -(\rho_c V + A) \langle \dot{\theta} u_l \rangle + \langle F_\omega^t \rangle - \Delta \rho V g \langle \cos \theta \rangle = 0, \quad (3.26) $$

with

$$ \langle F_{\omega}^{t} \rangle = \langle F_{D}^{t} \rangle + \langle F_{\Gamma_{T}}^{t} \rangle + \langle F_{\Gamma_{R}}^{t} \rangle + \langle F_{H}^{t} \rangle, \quad (3.27) $$

where, as already defined, the notation $\langle Q \rangle$ corresponds to the average over a time period of the quantity $Q$. Figure 13 presents the results obtained for the four density ratios and the seven elongation ratios, forming a total of 22 cylinders of parameters $(\xi, \rho_c/\rho_f, Ar_{3-D})$. The force model is also applied to the data from Toupoint (2018) to extend the parameter range. The sum $\sum \langle F^t \rangle$ is normalized for each cylinder with the magnitude of the buoyancy force, $\Delta\rho Vg$, where $\Delta\rho = (\rho_c - \rho_f)$. A satisfactory equilibrium is reached, resulting in average error 7%. Its main strength is to cover the whole range of parameters investigated, the departure from 0 being mostly due to the variability from one experiment to the other for given $(\xi, \rho_c/\rho_f, Ar_{3-D})$. The experimental cases for plastic cylinders made of POMC ($\rho_c/\rho_f = 1.4$) provide in average a less accurate balance, probably due to shape imperfections, as already discussed at the end of § 3.2.1.

We can now look at the mean contributions of the forces proposed to model the vortical force $F_\omega^t$. For the two density ratios (aluminium and titanium cylinders), figure 14 presents the mean contributions along the transverse direction for the proposed drag force, the lift forces in translation and rotation, and the history force, when the elongation ratio of the cylinder varies. As in the previous figure, the values of the forces are normalized with the magnitude of the buoyancy force $\Delta\rho Vg$. The inertial force and buoyancy force contributions, already determined and presented in figure 10, are also plotted for comparison. We can observe that the history force has no contribution to the mean force balance along the transverse direction, as expected. Moreover, we observe that the mean contributions of both the drag and circulation-induced forces are necessary to close the mean force balance.
D. Letessier, V. Roig and P. Ern

<figure><img src="image_13.png" /><figcaption>Figure 13. Mean balance in the transverse direction, averaged over a time period, for different elongation ratios and density ratios. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14. Mean contributions to the force balance in the transverse direction, measured forces and modelled ones, as a function of the elongation ratio of the cylinders, for two density ratios. Filled symbols indicate present study; open symbols indicate data from Toupoint (2018).</figcaption></figure>

### 3.3.2. *Instantaneous force balance experienced by the body along its path*

We now turn our attention to the instantaneous force balance experienced by the body along its path. We first compare the temporal evolution over a fluttering period of the vortical force $F_\omega$ obtained from the experimental measurements to that obtained from the model. Figure 15 presents this comparison for aluminium cylinders featuring five different elongation ratios. As a first point, we can see from both the experimental data and the model that the lowest magnitude of the vortical force is experienced by the body close to an extremum of its path, whatever $\xi$. While differences between the experiments and the model are visible (the largest difference being observed for the highest elongation ratio of the study, $\xi = 20$), we can also see that the model is able to stay close to the experimental determination all along a fluttering period for aluminium cylinders with different elongation ratios. For the shortest cylinder investigated, $\xi = 3$ (not
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_15.png" /><figcaption>Figure 15. Comparison of the temporal evolutions over a fluttering period of the measured vortical force to those obtained with the model for aluminium cylinders of various elongation ratios. Ten cylinders are drawn in black at regular time intervals. The trajectory of the centre of gravity of the cylinder is plotted in blue.</figcaption></figure>

presented here), the model does not provide satisfactory results. Also, the instantaneous force balance is not verified satisfactorily for the fluttering titanium cylinders, though the model provides, for all density ratios and elongation ratios investigated here, a good closure of the mean force balance (figure 13). To extend the range of validity of the model, more sophisticated dependences of the model parameters on the control parameters could be elaborated, for instance by improving the dependence of $C_d$, $C_T$ and $C_R$ on Re (the decrease as $Re^{-1}$ requiring a correction when Re is sufficiently large). Also, different formulations for the forces could be considered to improve the present modelling, especially for the history force, which plays a crucial role in the instantaneous force balance.

In spite of its limitations, the model provides a rather good prediction of the mean components of the hydrodynamical forces acting on cylinders featuring different materials and elongation ratios, falling in a fluttering or rectilinear regime. It also provides a vision of the different forces exerted instantaneously on fluttering aluminium cylinders along their paths. Figure 16 displays the instantaneous contributions of the three forces $F_D$, $F_Γ$ and $F_H$ involved in the model for the vortical force $F_ω$. In magnitude, the circulation-related forces are stronger than the quasi-steady drag term. For all the elongation ratios considered, the translation-induced lift is weaker than the rotation-induced lift. However, both contribute in the same direction, except when the cylinder inclination reaches a maximum. At this stage, the translation-induced lift changes sign, slightly before the rotation-induced
<figure><img src="image_16.png" /><figcaption>Figure 16. Illustration of the different contributions involved in the model for the vortical force for aluminium cylinders with different elongation ratios.</figcaption></figure>

lift. By definition, the relative weight between the drag and lift contributions is related to the body degrees of freedom and to the control parameters. As shown in figure 7(b), the amplitude of oscillation of the inclination decreases for aluminium cylinders when the elongation ratio $\xi$ increases from 5 to 20. The inclination angle of the velocity vector also decreases, reducing the lift contribution relative to that of the drag. The important contribution of the history force is also conspicuous.

In addition, figure 17 displays, over a fluttering period, the instantaneous contributions of the inertial forces in translation and rotation. As expected, these become weaker as the amplitudes of oscillation of the body degrees of freedom decrease ($\xi$ increases). We can also observe in this representation that the mean contribution of the nonlinear inertial force due to rotation $F_{IR}$ has opposite sign relative to the buoyancy force, as already seen in figure 10. The inertia term coupling the rotation of the body with its longitudinal velocity therefore behaves as an upward driving force on the body.

Finally, a more detailed illustration of the force balance experienced by the fluttering cylinders is presented in figure 18. The figure shows the temporal evolution of the different forces, measured and modelled, acting along the longitudinal and transverse directions for an aluminium cylinder with $\rho_c/\rho_f = 2.7$ and $\xi = 10$. Both longitudinal and transverse directions display complex fluctuating dynamics, and complex shapes for the different force components. For both directions, the sum of the different force
Fluttering motion of a confined freely falling cylinder

<figure><img src="image_17.png" /><figcaption>Figure 17. Illustration of the evolution over a fluttering period of the inertial force in translation and rotation, of the buoyancy force (represented only for the cylinder at the highest position) and of the vortical force for aluminium cylinders with different elongation ratios.</figcaption></figure>

components (expected to be zero) is plotted with a black dashed line. Measured forces are
those provided by the Kelvin–Kirchhoff equations, namely the buoyancy force F<sub>B</sub> and the
inertial forces in translation and rotation, respectively F<sub>IT</sub> and F<sub>IR</sub>. These exhibit important
fluctuations along the longitudinal direction, particularly the buoyancy force and the force
corresponding to translational inertia. We can also see that in this acceleration term, higher
frequencies emerge. Yet a rather satisfactory equilibrium of forces is reached with the
model for F<sub>ω</sub>. The fluttering dynamics is well captured by the model, since the remaining
fluctuations are essentially at a much larger frequency than for fluttering. The different
contributions introduced to model F<sub>ω</sub> have amplitudes comparable to one another. Good
agreement is in particular reached for the transverse force balance, for which both mean
and fluctuating components have to be compensated. We can observe that the introduction
of the history term is necessary to complement the fluctuations of the quasi-steady drag
term in this direction. Interestingly, the projection on both the transverse and longitudinal
directions of the circulation-induced lift, in particular the contribution associated with
translational effects, yields strong longitudinal fluctuations. At variance with the transverse
direction, a history force has to be introduced in the longitudinal direction to compensate
partly for these fluctuations.
<figure><img src="image_18.png" /><figcaption>Figure 18. Temporal evolutions of the (measured and modelled) forces experienced by an aluminium cylinder with $\xi = 10$ along its fluttering path. (a) Projection on the longitudinal axis of the cylinder. (b) Projection on the transverse axis of the cylinder.</figcaption></figure>

4. Conclusion

We first investigated experimentally the planar paths displayed by single cylinders falling freely in a thin-gap cell containing liquid at rest. For a significant range of control parameters (elongation ratio of the cylinder, solid-to-fluid density ratio, and Archimedes number), we determined the characteristics of the cylinder kinematics. Their analysis revealed that high-amplitude oscillations of the order of magnitude of the mean vertical velocity of the cylinder are observed for both translational and rotational velocities. The experimental results further showed that the mean fall velocity $\overline{u}_v$ does not scale with the gravitational velocity, which overestimates $\overline{u}_v$ and is, in particular, unable to capture the influence of the density ratio on it. We therefore developed a model able to account for the cylinder dynamics under these kinematical conditions, also encountered for distinct freely falling or rising bodies in confined or unconfined configurations. To model the body behaviour, we proposed a force balance that includes proper and added inertia terms, the buoyancy force and vortical contributions. Averaging the equations over a temporal period provides a mean force balance that governs the mean fall velocity of the cylinder, showing that the coupling between the translational and rotational velocity fluctuations induces a mean upward inertial force responsible for the decrease of $\overline{u}_v$. The consideration of this contribution enables us to take into account in the force balance the effect of the density ratio where it is expected physically, in an inertia term and not in a drag
Fluttering motion of a confined freely falling cylinder

term. This mean force balance also provides a scaling for the frequency of oscillation of heavier cylinders in agreement with the experimental measurements of the present study and with those of previous studies in confined and unconfined environments (Marchildon et al. 1964; Chow & Adams 2011; D'Angelo et al. 2017; Toupoint 2018). The mean force balance presented here may therefore be generally relevant for bodies displaying velocity fluctuations comparable to their mean fall velocity in confined and unconfined situations. We then considered the instantaneous force balance experienced by the body, and proposed three contributions for the vortical force modelling. These can be interpreted as drag, lift and history forces, and expressions for their dependence on the control parameters were adjusted on the basis of the experimental evolutions of the body’s degrees of freedom. The drag and lift formulations developed here extend ideas conveyed in particular by Andersen et al. (2005) (and references therein) for falling plates, enabling us to satisfactorily close the mean force balance acting on the cylinders. However, an additional zero-mean force, proposed under the form of a history force that breaks the instantaneous dependence of the vortical force, is necessary in our configuration to close the instantaneous force balance experienced by the cylinder.

**Acknowledgements.** We acknowledge financial support from ANR under the project ‘Muscats’ (ANR-19-CE05-0010). We also thank C. Toupoint for sharing his experimental data.

**Declaration of interests.** The authors report no conflict of interest.

Author ORCIDs.

Dylan Letessier https://orcid.org/0009-0006-2675-9631;

Véronique Roig https://orcid.org/0000-0002-9315-0467;

Patricia Ern https://orcid.org/0000-0002-6431-9083.

REFERENCES

ANDERSEN, A., PESAVENTO, U. & WANG, Z.J. 2005 Unsteady aerodynamics of fluttering and tumbling plates. *J. Fluid Mech.* **541** (1), 65–90.

ASSEN, M.P.A., WILL, J.B., NG, C.S., LOHSE, D., VERZICCO, R. & KRUG, D. 2024 Rising and settling 2-D cylinders with centre-of-mass offset. *J. Fluid Mech.* **981**, A7.

BASSET, A.B. 1888 *A Treatise on Hydrodynamics: with Numerous Examples*, vol. 2. Deighton, Bell and Co.

BATECHELOR, G.K. 1970 Slender-body theory for particles of arbitrary cross-section in Stokes flow. *J. Fluid Mech.* **44** (3), 419–440.

BELMONTE, A., EISENBERG, H. & MOSES, E. 1998 From flutter to tumble: inertial drag and Froude similarity in falling paper. *Phys. Rev. Lett.* **81** (2), 345–348.

BLEVINS, R.D. & PLUNKETT, R. 1980 Formulas for natural frequency and mode shape. *J. Appl. Mech.* **47** (2), 461–462.

CABRERA, F., SHEIKH, M.Z., MEHLIG, B., PLIHON, N., BOURGOIN, M., PUMIR, A. & NASO, A. 2022 Experimental validation of fluid inertia models for a cylinder settling in a quiescent flow. *Phys. Rev. Fluids* **7** (2), 024301.

CHOW, A.C. & ADAMS, E. 2011 Prediction of drag coefficient and secondary motion of free-falling rigid cylindrical particles with and without curvature at moderate Reynolds number. *J. Hydraul. Engng* **137**, 1406–1414.

CLIFT, R., GRACE, J.R. & WEBER, M.E. 1978 *Bubbles, Drops, and Particles*. Academic Press.

D'ANGELO, M.V., CACHILE, M., HULIN, J.-P. & AURADOU, H. 2017 Sedimentation and fluttering of a cylinder in a confined liquid. *Phys. Rev. Fluids* **2** (10), 104301.

ERN, P., RISSO, F., FABRE, D. & MAGNAUDET, J. 2012 Wake-induced oscillatory paths of bodies freely rising or falling in fluids. *Annu. Rev. Fluid Mech.* **44** (1), 97–121.

ERN, P., RISSO, F., FERNANDES, P.C. & MAGNAUDET, J. 2009 Dynamical model for the buoyancy-driven zigzag motion of oblate bodies. *Phys. Rev. Lett.* **102** (13), 134505.

GAI, G. & WACHS, A. 2024 On the dynamics and wakes of a freely settling Platonic polyhedron in a quiescent Newtonian fluid. *J. Fluid Mech.* **985**, A3.

GIANORIO, L., D'ANGELO, M.V., CACHILE, M., HULIN, J.-P. & AURADOU, H. 2014 Influence of confinement on the oscillations of a free cylinder in a viscous flow. *Phys. Fluids* **26** (8), 084106.
HOWE, M.S. 1995 On the force and moment on a body in an incompressible fluid, with application to rigid bodies and bubbles at low and high Reynolds numbers. *Q. J. Mech. Appl. Maths* **48**, 401–426.

INOUE, O. & SAKURAGI, A. 2008 Vortex shedding from a circular cylinder of finite length at low Reynolds numbers. *Phys. Fluids* **20** (3), 033601.

JAYAWEERA, K.O.L.F. & COTTIS, R.E. 1969 Fall velocities of plate-like and columnar ice crystals. *Q. J. R. Meteorol. Soc.* **95** (406), 703–709.

KHARROUBA, M., PIERSON, J.-L. & MAGNAUDET, J. 2021 Flow structure and loads over inclined cylindrical rodlike particles and fibers. *Phys. Rev. Fluids* **6** (4), 044308.

LAMB, H. 1932 *Hydrodynamics*. Cambridge University Press.

LAWRENCE, C.J. & WEINBAUM, S. 1986 The force on an axisymmetric body in linearized, time-dependent motion: a new memory term. *J. Fluid Mech.* **171** (1), 209.

LI, H., GOODWILL, T., WANG, J.Z. & RISTROPH, L. 2022 Centre of mass location, flight modes, stability and dynamic modelling of gliders. *J. Fluid Mech.* **937**, A6.

LOEWENBERG, M. 1993 Stokes resistance, added mass, and Basset force for arbitrarily oriented, finite-length cylinders. *Phys. Fluids* **5**, 765–767.

LORITE-DIEZ, M., ERN, P., CAZIN, S., MOUGEL, J. & BOURGUET, R. 2022 An experimental study of flow-structure interaction regimes of a freely falling flexible cylinder. *J. Fluid Mech.* **946**, A16.

MARCHILDON, E.K., CLAMEN, A. & GAUVIN, W.H. 1964 Drag and oscillatory motion of freely falling cylindrical particles. *Can. J. Chem. Engng* **42** (4), 178–182.

MORICHE, M., FLORES, O. & GARCÍA-VILLALBA, M. 2017 On the aerodynamic forces on heaving and pitching airfoils at low Reynolds number. *J. Fluid Mech.* **828**, 395–423.

MORICHE, M., UHLMANN, M. & DUŠEK, J. 2021 A single oblate spheroid settling in unbounded ambient fluid: a benchmark for simulations in steady and unsteady wake regimes. *Intl J. Multiphase Flow* **136**, 103519.

MOUGIN, G. & MAGNAUDET, J. 2002 The generalized Kirchhoff equations and their application to the interaction between a rigid body and an arbitrary time-dependent viscous flow. *Intl J. Multiphase Flow* **28** (11), 1837–1851.

PESAVENTO, U. & WANG, Z.J. 2004 Falling paper: Navier–Stokes solutions, model of fluid forces, and center of mass elevation. *Phys. Rev. Lett.* **93** (14), 144501.

PIERSON, J.-L. 2023 Inertial settling of an arbitrarily oriented cylinder in a quiescent flow: from short-time to quasi-steady motion. *arXiv:2303.15909*.

PIERSON, J.-L., AUGUSTE, F., HAMMOUTI, A. & WACHS, A. 2019 Inertial flow past a finite-length axisymmetric cylinder of aspect ratio 3: effect of the yaw angle. *Phys. Rev. Fluids* **4** (4), 044802.

PIERSON, J.-L., KHARROUBA, M. & MAGNAUDET, J. 2021 Hydrodynamic torque on a slender cylinder rotating perpendicularly to its symmetry axis. *Phys. Rev. Fluids* **6** (9), 094303.

PRUPPACHER, H.R., CLAIR, B.P. & HAMIELEC, A.E. 1970 Some relations between drag and flow pattern of viscous flow past a sphere and a cylinder at low and intermediate Reynolds numbers. *J. Fluid Mech.* **44** (4), 781–790.

SAXBY, J., BECKETT, F., CASHMAN, K., RUST, A. & TENNANT, E. 2018 The impact of particle shape on fall velocity: implications for volcanic ash dispersion modelling. *J. Volcanol. Geotherm. Res.* **362**, 32–48.

SEMIN, B., DECOENE, A., HULIN, J.-P., FRANÇOIS, M.L.M. & AURADOU, H. 2012 New oscillatory instability of a confined cylinder in a flow below the vortex shedding threshold. *J. Fluid Mech.* **690**, 345–365.

TAGLIAVINI, G., MCCORQUODALE, M., WESTBROOK, C., CORSO, P., KROL, Q. & HOLZNER, M. 2021 Drag coefficient prediction of complex-shaped snow particles falling in air beyond the Stokes regime. *Intl J. Multiphase Flow* **140**, 103652.

TCHOUFAG, J., FABRE, D. & MAGNAUDET, J. 2014 Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders. *J. Fluid Mech.* **740**, 278–311.

TOUPOINT, C. 2018 Path and wake of cylinders falling in a liquid at rest or in a bubble swarm towards the hydrodynamical modeling of ebullated bed reactors. PhD thesis, INPT, Toulouse, France. https://theses.fr/2018INPT0126.

TOUPOINT, C., ERN, P. & ROIG, V. 2019 Kinematics and wake of freely falling cylinders at moderate Reynolds numbers. *J. Fluid Mech.* **866**, 82–111.

WILL, J.B., MATHAI, V., HUISMAN, S.G., LOHSE, D., SUN, C. & KRUG, D. 2021 Kinematics and dynamics of freely rising spheroids at high Reynolds numbers. *J. Fluid Mech.* **912**, A16.

YANG, Y., FENG, Z. & ZHANG, M. 2022 Onset of vortex shedding around a short cylinder. *J. Fluid Mech.* **933**, A7.

ZHOU, W., CHRUST, M. & DUŠEK, J. 2017 Path instabilities of oblate spheroids. *J. Fluid Mech.* **833**, 445–468.