2.1. Mathematical formulation

The mathematical formulation and the numerical method are basically the same as that of Chrust et al. (Reference Chrust, Bouchet and Dušek2013) or Zhou et al. (Reference Zhou, Chrust and Dušek2017). The details can be found in the PhD thesis by Chrust (Reference Chrust2012). Here we sum up the main features and present the minor improvements implemented to tackle the configuration of falling (rising) elongated cylinders.

The fluid equations are solved in a vertical cylindrical domain $\varOmega$ accompanying the translation motion of the cylinder with respect to a fixed frame $(0_{fix},x_{fix},y_{fix},z_{fix})$ (see figure 1a). The centre $O$ of the moving frame $(O,x_c,y_c,z_c)$ with axes parallel to the fixed frame moves with velocity $\boldsymbol {u}$ with respect to the fixed frame. The figure represents both the situation of falling and rising cylinders. In the latter case, the effective gravity is still represented as downwards pointing and the bottom cylinder face is always the inflow. Testings in previous work showed that the following dimensions make the results insensitive to the position of outer boundaries: $L_u=12 d$, $L_d=25 d$ and $R_c=8 d$. The flow velocity field $\boldsymbol {v}$ is measured with respect to the fixed frame. As the result, the inflow boundary condition imposed at $z=-L_u$ and at the outer cylinder boundary simulate a quiescent fluid of velocity $\boldsymbol {v}=0$. At the outflow face $z=L_d$ a Neumann no-stress condition is used. A spherical subdomain $\varOmega _s$ of radius $R_s$ is centred, again, at $O$ but the axis $O z_s$ is defined by the cylinder axis so that the whole subdomain rotates with the cylinder (see figure 1b). The whole domain $\varOmega$ is decomposed into $\varOmega _s$ and $\varOmega _{out}$ where $\varOmega _{out}$ is the outer complement of $\varOmega _s$ in $\varOmega$.

Figure 1. (a) Outer cylindrical domain. The circle of radius $R_s$ represents the boundary of the inner spherical domain of (b). Here $g_{eff}$ represents the direction and orientation of the vector of effective gravity, $(0_{fix},x_{fix},y_{fix},z_{fix})$ is a fixed frame with respect to which all velocities are expressed and $(O,x_c,y_c,z_c)$ is the local frame in which the cylindrical domain is discretised. (b) Inner spherical domain of radius $R_s$. Here $(O,x_s,y_s,z_s)$ is a rotating frame with $O z_s$ axis along the cylinder symmetry axis.

The flow equations are non-dimensionalised using the cylinder diameter $d$ as length scale and the velocity scale

(2.1)\begin{equation} U_{ref}=\sqrt{V^*g_{eff} d}, \end{equation}

where $V^*$ is the non-dimensionalised volume:

(2.2)\begin{equation} V^*= \frac{V}{d^3}= \frac{\rm \pi}{4} \frac{L}{d}. \end{equation}

This results in

(2.3)\begin{equation} \frac{\partial {\boldsymbol{v}}}{\partial t}+ \left[(\boldsymbol{v} - \boldsymbol{u}- \boldsymbol{\omega} \times \boldsymbol{r}) \boldsymbol{\cdot} \boldsymbol{\nabla} \right] {\boldsymbol{v}} + \boldsymbol{\omega} \times \boldsymbol{v} ={-}\boldsymbol{\nabla} p + \frac{1}{Ga} {{\nabla}}^2 {\boldsymbol{v}} \end{equation}

and

(2.4)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\boldsymbol{v}} = 0, \end{equation}

where $\boldsymbol {u}$ is the translation velocity of the centre of the cylinder and $\boldsymbol {\omega }$ is the angular velocity of the rotating frame in $\varOmega _s$ and zero in $\varOmega _{out}$. The velocity scale (2.1) yields the Galileo number defined, in the case of cylinders, as

(2.5)\begin{equation} Ga= \frac{\sqrt{\dfrac{\rm \pi}{4} L d^2 g_{eff}}}{\nu}. \end{equation}

It should be noted that the so-defined Galileo number differs only by the factor of $\sqrt {{\rm \pi} /6}$ from the Archimedes number of Toupoint et al. (Reference Toupoint, Ern and Roig2019) ($Ga=\sqrt {{\rm \pi} /6} Ar$). The no-slip boundary condition on the cylinder surface amounts to setting

(2.6)\begin{equation} \boldsymbol{v}|_S = \boldsymbol{u} + \boldsymbol{\varOmega} \times \boldsymbol{r}|_{S}, \end{equation}

where $\boldsymbol {r}|_{S}$ is the position vector of a point of the cylinder surface in the frame of the spherical domain and $\boldsymbol {\varOmega }$ the angular velocity of the cylinder. The cylinder is let to rotate with respect to the frame $(O x_s y_s z_s)$ so that $\boldsymbol {\varOmega }$ is only parallel to $\boldsymbol {\omega }$ but not equal.

The flow equations are coupled with the solid body motion equations solved in the spherical subdomain in the frame defined by the cylinder axis:

(2.7)$$\begin{gather} m^* \left(\frac{{\rm d} \boldsymbol{u} }{{\rm d}t}+ \boldsymbol{\omega} \times \boldsymbol{u} \right) = \boldsymbol{F}_{fl} (\boldsymbol{v},p) - \boldsymbol{k}_{fix}, \end{gather}$$

(2.8)$$\begin{gather}\boldsymbol{I}^* \frac{{\rm d} {\boldsymbol{\varOmega}} }{{\rm d}t} + \boldsymbol{\omega} \times (\boldsymbol{I}^* \boldsymbol{\varOmega}) = \boldsymbol{M}_{fl}, \end{gather}$$

where $m^*$ stands for the non-dimensionalised mass

(2.9)\begin{equation} m^*=\frac{\rho_s}{\rho} V^*, \end{equation}

$\boldsymbol {I}^*$ for the diagonal non-dimensionalised moment of inertia

(2.10)\begin{equation} \boldsymbol{I}^*= \text{diag} (I^*_\perp, I^*_\perp, I^*_\parallel); \quad I^*_\perp = \frac{1}{16} m^*\left(1+\frac{4}{3} \left(\frac{L}{d}\right)^2\right),\quad I^*_\parallel = \frac{1}{8} m^* \end{equation}

and $\boldsymbol {k}_{fix}$ is the unit vector pointing opposite to the effective gravity. $\boldsymbol {F}_{fl}$ and $\boldsymbol {M}_{fl}$ are the hydrodynamic force and torque obtained by integration of the pressure and shear stress on the body surface.

In what follows, the Galileo number (2.5), the density ratio $\rho _s/\rho$ and the aspect ratio $L/d$ are used as parameters of the investigation.

2.2. Numerical method

The spectral–spectral-element discretisation consisting of the spectral element discretisation of the axial–radial plane combined with Fourier azimuthal expansion (Ghidersa & Dušek Reference Ghidersa and Dušek2000) is applied separately in both subdomains reconnected at the spherical surface of radius $R_s$ by using the spherical function expansion and the Wigner matrices of the rotation group representation. The truncation of the spherical function and azimuthal expansions is defined by the highest degree $\ell _{max}$ of spherical functions $Y_{\ell,m}$ (Hecht Reference Hecht2000). If no information is to be lost, the azimuthal expansion is to be truncated also at the wavenumber $m_{max}=\ell _{max}$. An example of the axial–radial spectral-element discretisation is presented in figure 2 for the cylinder of aspect ratio $L/d=3$.

Figure 2. Enlarged view of the spectral element mesh of the spherical subdomain delimited by the interface represented by the thick half-circle of radius $R_s$ (here $R_s=2.5$). The cylinder (its half-section $-1.5 \le z \le 1.5, r \le 0.5$) is represented by the filled rectangle. The immediate neighbourhood outside of the gliding interface is also visible. Individual spectral elements are discretised by $6 \times 6$ collocation points.

Among the numerical parameters needed to be optimised, there appeared to be the radius $R_s$ of the spherical subdomain enclosing the cylinder, the number of elements along the interface $N_i$ (uniform distribution is adopted), the total number of elements $N$ and their distribution and the truncation of the spectral expansions $\ell _{max}$. Since each shape requires a specific discretisation a variety of meshes were implemented (see table 1).

Table 1. Tested and mostly used discretisation parameters in the simulations: $R_s$, radius of the gliding interface; $N_i$, number of elements at the interface; $N$, total number of spectral elements; $\ell _{max}$, maximum degree of spherical function expansion (equal to the largest wavenumber of azimuthal expansion).

The major issue of the present numerical study was the treatment of the largest aspect ratios of 3 and 5. The method was developed and used, so far, for configurations involving axisymmetry breaking. As a result, cases when the axisymmetry axes ($O z_c$ and $O z_s$) of the subdomains were mutually perpendicular rarely occurred and, when they occurred, then only for short times. In the present case, the axes are constantly perpendicular in the equilibrium position below the threshold of the onset path instabilities. In addition, the rotating spherical subdomain must be large enough to contain the cylinder. This situation was already identified as unfavourable by Chrust (Reference Chrust2012) who showed, however, that increasing the spectral accuracy (enlarging $\ell _{max}$) solves the problem. The necessity to account for high degrees $\ell$ of spherical functions made us check two identities: the orthogonality of Wigner matrices $d^\ell _{m m'}({\rm \pi} /2)$ involved in the formula for rotation matrices

(2.11)\begin{equation} \sum_{m'={-}\ell}^\ell d^\ell_{m m'} d^\ell_{m'' m'} =\delta_{m,m''} \end{equation}

and the numerical accuracy of the orthogonality of Legendre polynomials

(2.12)\begin{equation} \frac{2 \ell+1}{2} \frac{(\ell-m)!}{(\ell+m)!} \int_0^{\rm \pi} P_\ell^m(\cos \varphi) P_{\ell'}^m(\cos \varphi) \sin \varphi \,{\rm d}\varphi = \delta_{\ell,\ell'}. \end{equation}

Having replaced the basic formula of Hecht (Reference Hecht2000) involving factorials by the recursion of Blanco, Flórez & Bermejo (Reference Blanco, Flórez and Bermejo1997), we satisfy (2.11) up to $\ell =63$ practically with machine precision. As for (2.12), the tests yielded the results summed up in table 2.

Table 2. Results of testing of the numerical accuracy of discretisation of the orthogonality relation (2.12). ‘Mesh no.’ refer the meshes described in table 1.

The results on lines 2 and 4 of the table might cast doubt if the accuracy is sufficient. For this reason, several tests have been executed to compare results of simulation of the same oscillating regimes for meshes $2, 3$ and $5,6$. The largest variation of vertical velocity was 0.16 % and 0.3 % of the amplitude of oscillations when mesh 2 was compared with mesh 3 (much less for meshes $5,6$). Other tests concerned the spectral-element resolution. A reliable test (testing also the spectral element distribution) consists of increasing the number of collocation points. An $8\times 8$ collocation point refinement was tested with a less than 0.2 % effect on the vertical velocity.

The parameter to which the results appeared to be the most sensitive was the truncation $\ell,m \le \ell _{max}$. To avoid too tedious testing and to benefit from the used fast-Fourier-transform (FFT) algorithm optimised for powers of two, the truncations were varied by a factor of two, i.e. with taking account of the numbering from 0, $\ell _{max}=15$, $31$ and 63. Many results have been checked with two different truncations. In particular, most of the presented results concerning cylinders with aspect ratio $L/d=2$ and 3 have been calculated both with truncation $\ell _{max}=15$ and 31 (the more accurate being selected). The average variation of the vertical velocity was of the order of tenths of a per cent. For $L/d=5$, the (expensive) truncation at $\ell _{max}=63$ has been compared with $\ell _{max}=31$ (meshes 4 and 5). The difference being of less than 0.1 %, the truncation $\ell _{max}=31$, was mostly kept as sufficient. In spite of these quantitatively satisfactory results any unexpected behaviour, especially bi-stabilities, was checked for numerical robustness by varying the spectral truncation.

3.1. Flutter

Flutter, in the sense we take up from Toupoint et al. (Reference Toupoint, Ern and Roig2019), is the most conspicuous and most frequently observed behaviour in literature. It involves periodic oscillations of the cylinder axis in a fixed vertical plane, selected by initial conditions. The whole motion is thus planar, with the cylinder centre describing a zig-zagging (almost vertical in the mean) planar trajectory. In our simulations, the flutter sets in mostly via a supercritical bifurcation from the steady vertical state but, for aspect ratios $L/d=3$ and 5 and large density ratios ($\rho _s/\rho =5,10$), subcritical bifurcations with coexistence of large flutter with non-oscillating or weakly oscillating states have been observed.

In the simulations, we observe the trajectory to settle to a fixed plane even if the initial condition does not correspond to a symmetric flow. Here, we define the horizontal $x$-axis in the direction of the plane and the $y$-axis as normal to the trajectory plane. In this reference frame, $u_x$ and $n_x$ denote the horizontal velocity and horizontal projection of the unit vector along the symmetry axis of the cylinder in the trajectory plane, $u_y, n_y$ are the projections normal to the plane and $u_z,n_z$ are the corresponding vertical components. Here $n_z$ represents the sine of the inclination angle of the cylinder axis with respect to the horizontal orientation. With this definition of the reference frame, the main features of this regime are:

1. (i) planar symmetry with respect to the trajectory plane $xOz$;

2. (ii) zero $u_y$ velocity;

3. (iii) significant periodic oscillations of the cylinder (more exactly of its symmetry axis) in the trajectory plane characterised by oscillations of $u_x$ and $n_z$ about (almost) zero;

4. (iv) oscillations of the vertical velocity $u_z$ about a mean value (which will be taken to evaluate the mean Reynolds number and the mean drag).

These characteristics are in perfect agreement with those originating in experimental observation of Toupoint et al. (Reference Toupoint, Ern and Roig2019).

A complementary remark concerning the choice of the numerical resolution is to be made at this place. For this purpose, we take the case of the cylinder of aspect ratio $L/d=3$ and density ratio 1.16 (chosen to agree with the experiment of Toupoint et al. (Reference Toupoint, Ern and Roig2019). Table 3 presents results of the simulation of the flutter at $Ga=200$ with the truncation $\ell _{max}=15$ and 31.

Table 3. Flutter of a cylinder $L/d=3$, $\rho _s/\rho =1.16$ simulated at Galileo number 200 for two numerical resolutions given by truncation $\ell _{max}$. $A(n_z)$, amplitude of the vertical projection of the cylinder axis; mean($n_z$), its mean value; mean($u_z$), mean value of the vertical velocity; $A(u_x)$, amplitude of the horizontal velocity; mean($u_x$), its mean value and period of oscillations (in units $d/U_{ref}$).

It can be seen that the reported values are quite acceptable even for the lower resolution. The comparison shows that the small mean values are not to be assigned to an insufficient convergence of the spectral method. The axial–radial spectral element mesh of the spherical subdomain is probably not symmetric within machine precision The non-zero mean value of the horizontal velocity results in a slightly oblique zigzagging trajectory, albeit with a very small inclination mean($u_x$)/mean($u_z$), typically of 10$^{-4}$ and less (see figure 3a). This obliqueness is systematically neglected. The better resolution is required mainly by the necessity to avoid the instability breaking the planar symmetry as seen in figure 3(b,c). Whereas the lower resolution yields a growth rate of oscillation amplitude of the velocity normal to the trajectory plane of 0.016 leading to a spurious quasi-periodic three-dimensional trajectory, figure 3(c) shows that, with the better resolution, the symmetry is well preserved with numerical oscillations of the normal velocity of constant amplitude of order of $10^{-7}$.

Figure 3. (a) View of the trajectory of the cylinder centre with extremely enhanced horizontal scale, (b) amplitude of oscillations of the velocity $u_y$ normal to the trajectory plane at truncation $\ell _{max}=15$ and (c) velocity $u_y$ with truncation $\ell _{max}=31$. Here $L/d=3$, $\rho _s/\rho =1.16$ and $Ga=200$.

The main properties of the flutter, enumerated previously, are illustrated in figures 4–6 referring to the motion of a cylinder of aspect ratio $L/d=5$ and density ratio $\rho _s/\rho =10$ at Galileo number $Ga=150$. Figure 4(a) shows a kinogram of the cylinder motion over one period of oscillations. In this case of high density ratio, the maximum inclination is reached at the maximum deviation from the mean vertical line, i.e. when the horizontal velocity becomes zero. Figure 4(b) represents the vertical (streamwise) component of vorticity $\omega _z$ in the wake at two opposite levels (level 0.28 has been chosen). The wake structure is characterised by two exactly anti-symmetric vortex sheets shed from the cylinder surface. Note the secondary vortices of opposite sign shed from the cylinder tip. Figure 5(a) shows a common graph of the cylinder inclination $n_z$ and of the horizontal velocity (the latter has been rescaled to obtain comparable amplitudes). The phase shift is close to ${\rm \pi} /2$ meaning that the maximum velocity is reached when the cylinder inclination is zero (and conversely), in agreement with the statement concerning the kinogram. This phase shift was discussed by Toupoint et al. (Reference Toupoint, Ern and Roig2019) at a small density ratio (1.16) and was found close to $5{\rm \pi} /4$ (${\sim }3.93$). In the sense of the definition of Toupoint et al. (Reference Toupoint, Ern and Roig2019), our found phase shift is close to $3 {\rm \pi}/2$ (4.71), more accurately 4.80. However, the phase shift varies with the density ratio. For the case of figure 3 (the same density ratio as that of Toupoint et al. Reference Toupoint, Ern and Roig2019), we find the phase shift 4.09, whereas for the same aspect ratio $L/d=3$ and Galileo number (200) but density ratio of 10, the phase shift was found to be 4.53. The situation of the example in figure 5(a) can thus be considered as characteristic for dense and long cylinders. The stability of the symmetry is documented in figure 5(c). The simulation referred to by the graphs was run from an initial condition in which the cylinder was in an almost horizontal position. This allowed us to follow both the instability amplification and its saturation (see figure 5b). At $Ga=150$, the amplification presents a slightly subcritical growth (increasing growth rate; see figure 6). This is due to the proximity of a bi-stability of the flutter with a regime of small oscillations setting in close to $Ga=200$. The primary instability, yielding the flutter, sets in at $Ga \sim 100$ in this case. Closer to this threshold, the subcritical effect disappears (see the case of $Ga=120$ also represented in figure 6).

Figure 4. (a) Cylinder motion over a period of oscillation. (b) Vertical vorticity component $\omega _z$ in the wake represented by two iso-surfaces of opposite sign at $\omega _z=\pm 0.28$. Here $L/d=5$, $\rho _s/\rho =10$ and $Ga=150$.

Figure 5. (a) Comparison of the phase shift between the oscillations of the cylinder axis (represented in terms of the vertical component of its unit vector $n_z$) and those of horizontal velocity $u_x$. For better reading the velocity $u_x$ is rescaled by a factor of 5. (b) Growth and saturation of instability from an initial horizontal position. (c) Stability of the symmetry documented by the velocity $u_y$ normal to the trajectory plane. Here $L/d=5$, $\rho _s/\rho =10$ and $Ga=150$.

Figure 6. (a) Growth and saturation of the amplitude of oscillations $A(n_z)$ of the cylinder axis at $Ga=150$ and 120. (b) The growth rate $\gamma =\dot {A}(n_z)/A(n_z)$ vs time. Here $L/d=5$ and $\rho _s/\rho =10$.

In some cases of aspect ratio $L/d=3$ ($\rho _s/\rho =1.16, Ga=300$ and $\rho _s/\rho =2, Ga=190,200,250, 275$) we evidenced a quasi-periodic (maybe subharmonic) modulation. It was found to disappear with increasing Galileo number and to preserve the flow symmetry. For this reason we do not distinguish this exceptional behaviour from the strictly periodic flutter.

3.2. Fluid modes

The flutter described in the previous section is characterised by a strong solid–fluid interaction. If the motion of the cylinder is inhibited, this instability disappears. Tchoufag et al. (Reference Tchoufag, Fabre and Magnaudet2014) introduced the distinction between such ‘solid’ modes and ‘fluid’ modes arising due to instabilities of the fluid flow (wake). Coexisting solid and fluid modes have been evidenced for flat falling objects (e.g. Dušek, Zhou & Chrust Reference Dušek, Zhou and Chrust2021). Fluid modes, i.e. regimes where the cylinder weakly responses to the wake oscillations were evidenced by Toupoint et al. (Reference Toupoint, Ern and Roig2019) and all the three transitional regimes investigated by Lorite-Diez et al. (Reference Lorite-Diez, Ern, Cazin, Mougel and Bourguet2022) can be qualified as such. Oblong cylinders we focus on in this paper also present oscillations, of much smaller amplitude than the flutter, and are fully due to instabilities of the wake.

3.2.1. Wake of a fixed cylinder

To substantiate this statement, let us start with presenting the transition scenario in the wake of a fixed cylinder placed perpendicularly to a uniform flow. This configuration was already investigated by Inoue & Sakuragi (Reference Inoue and Sakuragi2008). Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) presented additional details on the wake of a cylinder of aspect ratio $L/d=3$ for the yaw angle of 90$^\circ$ in the appendix of their paper. They evidence four regimes that can be roughly identified with the following description summarising our results obtained for a cylinder of aspect ratio $L/d=5$. In agreement with the choice of frame of § 3.1, the $z$-direction is oriented streamwise, the $x$-axis points spanwise in the direction of the cylinder axis and the $y$-axis is oriented crosswise perpendicularly both to the flow and the cylinder axis. We evidenced:

1. (a) a steady doubly symmetric (with respect to $xOz$ and $yOz$ planes) regime with no lift at $Re \lesssim 120$;

2. (b) a periodic doubly symmetric regime with no lift and periodic drag for $120 \lesssim Re \lesssim 140$;

3. (c) a periodic regime with only spanwise symmetry (with symmetry with respect to $yOz$ plane), i.e. with a periodically oscillating crosswise lift $L_y$, for $140 \lesssim Re \lesssim 175$;

4. (d) a chaotic regime with both symmetries broken at $Re \gtrsim 170$.

The quantitative characteristics are presented in table 4.

Table 4. Results of simulation of hydrodynamic forces acting on a fixed cylinder of aspect ratio $L/d=5$ placed perpendicularly to an incoming uniform flow. The column ‘regime’ identifies one of the enumerated regimes. $Re$, Reynolds number; $C_D$, drag coefficient of the vertical force (mean value in unsteady regimes); $A(C_D)$, drag coefficient amplitude ($^{a}$ standard deviation); $A(C_{L_y})$, amplitude ($^{a}$ standard deviation) of the crosswise lift coefficient; $St_z, St_y$ Strouhal number based on the cylinder diameter and on the frequency of the drag and lift oscillations, respectively ($^{b}$ main peak of FFT).

3.2.2. Fluid modes of a freely falling cylinder

3.2.2.1 Symmetric periodic state

The symmetric periodic fluid mode can easily be identified as the (b)-regime of the wake of a fixed cylinder. Since the oscillations of the vertical velocity are extremely small (of the order of $10^{-4}$ of the mean vertical velocity), the threshold is practically independent of the density ratio. Table 5 lists some properties of the symmetric periodic states close to their threshold for cylinders of aspect ratio $2$, $3$ and 5. As can be seen, the case $L/d=5$, $Re=130$ of table 4 is in agreement with the case $L/d=5$, $Ga=200$ of table 5. In figure 7, we present the vortex structures in the wake of cylinder of aspect ratio $L/d=3$ falling at $Re=165.8$. Exactly spanwise and crosswise views are selected to demonstrate both symmetries. It is interesting to note that, in spite of the very small amplitude of the velocity oscillation, the vortex shedding is far from weak. Two vortex pairs are shed from the ends of the cylinder. Each pair is perfectly anti-symmetric crosswise and both pairs are perfectly synchronised satisfying accurately the spanwise anti-symmetry. This makes both components of the horizontal velocity (i.e. of the lift if the cylinder is fixed) vanish. The unsteady vortex shedding develops fully in the far wake and, thus, only weakly influences the hydrodynamic force acting on the cylinder. The symmetry implies also a zero torque so that the cylinder axis remains exactly horizontal.

Table 5. Properties of the symmetric periodic state slightly above its threshold (i.e. no oscillations exist at Galileo numbers $200$, $170$ and 190). $u_z$, mean vertical velocity; $A(u_z)$, amplitude of its oscillation (non-dimensionalised by $U_{ref}$ defined by (2.1)); $Re$ and $St$, Reynolds and Strouhal numbers based on the cylinder diameter and on the mean value of the vertical velocity ($Re=Ga \, u_z$, $St=1/(T \, u_z)$); $C_D$, drag coefficient based on the cross section $Ld$ and the mean vertical velocity ($C_D= 1/(0.5 ({L}/{d}) \,u_z^2$)).

Figure 7. Example of a symmetric periodic state: $L/d=3$, $\rho _s/\rho =5$, $Ga=190$. Spanwise (a) and crosswise (b) view of iso-surfaces of streamwise vorticity at levels $\pm 0.4$.

3.2.2.2 Crosswise-oscillating, spanwise ($yOz$) symmetric regime

The mode labelled (c) in § 3.2.1 generates a regime characterised by a crosswise periodic oscillation when the cylinder is free. As already noted by Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019), the reason is the loss of symmetry with respect to the plane defined by vertical direction and the cylinder axis. The loss of the crosswise symmetry is visible in figure 8(a). The characteristics of the cylinder motion are reported in the first line of table 6. The amplitude of the crosswise velocity is about a factor of 10 larger than that of the vertical oscillations but remains 0.36 % of mean vertical velocity, which is considered very small. Again, this is due to the fact that the symmetry breaking occurs rather far in the wake. Figure 8(a) shows that, up to about 6$d$ downstream of the cylinder, the vortex structure does not present visible asymmetry and alternate vortex shedding clearly appears only about 10$d$ downstream of the cylinder.

Figure 8. Iso-surfaces of vorticity $\omega _z$ for (a) $L/d=5$, $\rho _s/\rho =1$, $Ga=250$, $\omega _z=\pm 0.45$ (crosswise-oscillating, spanwise ($yOz$) symmetric regime), (b) $L/d=3$, $\rho _s/\rho =5$, $Ga=195$, $\omega _z=\pm 0.45$ (weak quasi-periodic flutter), (c) $L/d=5$, $\rho _s/\rho =0$, $Ga=300$, $\omega _z=\pm 0.5$ (asymmetric chaotic regime).

Table 6. Properties of secondary fluid modes. The additional columns as compared with table 5 show the amplitudes of $u_y$, $u_x$ and $n_x$.

$^{a}$Standard deviation.

$^{b}$Highest and second highest peak of frequency.

$^{c}$Based on the period of oscillations of the vertical projection of the cylinder axis $n_z$.

3.2.2.3 Weak quasi-periodic flutter

At aspect ratio $L/d=3$ and density ratios 5 and 10, a small periodic ‘flutter’, i.e. oscillations of the cylinder axis in a fixed vertical plane, sets in (with critical Galileo number between 190 and 195 at $\rho _s/\rho =5$ and between 220 and 230 at $\rho _s/\rho =10$) witnessing of a preserved crosswise ($xOz$) symmetry. The feature distinguishing this state from a genuine flutter is the coexistence of these oscillations with those of the vertical velocity characteristic for the symmetric periodic state. Since the frequency of the flutter is not synchronised with that of the vertical velocity, the mode is, strictly speaking, quasi-periodic.

At a density ratio of 5, three features distinguish this regime from a genuine flutter: (i) the vortex shedding characteristic for the symmetric fluid mode subsists; (ii) the vertical velocity oscillations and those of the cylinder axis are not synchronised, making the regime quasi-periodic even though the vertical oscillations of the cylinder axis are periodic; (iii) the weak quasi-periodic mode coexists (starting from $Ga=198.5$) with a genuine strong flutter. Figure 8(b) corresponds to the same aspect and density ratio as figure 7 showing the effect of the instability, setting in between $Ga=190$ and 195, on the wake structure. Quantitative characteristics can be found on the second line of table 6.

At the density ratio 10, the vertical oscillations of the cylinder axis are quasi-periodic themselves and the $xOz$ symmetry is not quite exact. The amplitude of oscillations is about three times smaller than at $\rho _s/\rho =5$. Since the branches do not coexist in the same interval, the comparison applies to $Ga=220$ for $\rho _s/\rho =5$ ($A(n_z)=0.17$) and $Ga=250$, $\rho _s/\rho =10$, for which the amplitude is indicated in table 6 (line 3). Unlike all other regimes, the dominant frequency does not correspond to that of the flutter (for the same density ratio) but is 30 % to 40 % higher. The frequency of the flutter is still present causing the quasi-periodicity (see figure 9a). (The coexisting flutter has an amplitude of more than a factor of 10 larger and a frequency of 0.02.) The wake structure is more complicated than in figure 8(b) but the $xOz$ symmetry is still well preserved (see figure 9b).

Figure 9. (a) The FFT of the vertical projection of the cylinder axis. (b) Streamwise vorticity field $\omega _z=\pm 0.5$. Here $L/d=3$, $\rho _s/\rho =10$ and $Ga=250$.

A similar state was also observed for $L/d=2$ and density ratio 5 at $Ga=230$ and 233. It seems rather to be an exceptional case of more ordered states between asymmetric quasi-periodic states described in the following.

3.2.2.4 Asymmetric quasi-periodic regime

No symmetry subsists in this regime (see figure 10a). As a result, the horizontal velocity and the cylinder axis oscillate both crosswise and spanwise, moreover with different frequency (figure 10b,c). However, as shown in figure 10(d), the dynamics is still quasi-periodic and not yet chaotic. The spanwise velocity $u_x$ (as well as the not represented vertical projection of the cylinder axis $n_z$) is practically periodic with period $T=18.84$ yielding a Strouhal number $St=0.0466$ (reported as the characteristic Strouhal number of this regime in table 6). The figure represents the frequencies as multiples of $f_1=1/T$. This shows that the third peak of the vertical velocity $u_z$ is situated at an incommensurate frequency $f_2=3.6/T$. Similarly, $u_y$ and $n_y$ present a subharmonic frequency of about $f_3=0.21/T$. The second highest peak of $u_y$ corresponds to $2 f_1 - f_3$.

Figure 10. Example of an asymmetric quasi-periodic regime, where $L/d=2$, $\rho _s/\rho =5$ and $Ga=220$. (a) Spanwise and crosswise profile of vertical flow velocity induced by the cylinder motion 2$d$ downstream of the cylinder centre. (b) Horizontal projection of the velocity of the cylinder centre on a time interval of $12 T$, $T$ being the period of oscillations of $u_x$ velocity and of $n_z$, the vertical projection of cylinder axis. (c) Projection of the unit vector of the cylinder axis onto the vertical and crosswise ($yOz$) plane (note that horizontal scale is a factor of 10 smaller than the vertical scale). (d) FFT of quantities specified in the legend. The horizontal axis is scaled in multiples of the frequency $1/T$.

3.2.2.5 Asymmetric chaotic regime

The fluid modes in all cases eventually become chaotic at Galileo numbers not much exceeding the threshold of the primary symmetric periodic state. As an example, figure 8(c) presents the vorticity in the wake of a cylinder of aspect ratio 5 and density ratio zero at $Ga=300$. Though the $n_z$-oscillations are still fairly periodic, the horizontal oscillations of the cylinder axis and those of the vertical and transverse cylinder velocity ($u_z, u_y$) are clearly chaotic (see figure 11).

Figure 11. Example of an asymmetric chaotic regime, where $L/d=5$, $\rho _s/\rho =0$ and $Ga=300$. (a) Horizontal projection of the velocity of the cylinder centre on a time interval of $10 T$, $T$ being the mean period of oscillations of the vertical projection of the cylinder axis $n_z$. (b) Projection of the unit vector of the cylinder axis onto the vertical and crosswise ($yOz$) plane. (d) The FFT of quantities specified in the legend. The horizontal axis is scaled in multiples of the frequency $1/T$.

4.1. $L/d=2$

Figure 12 represents the obtained results for this aspect ratio in the form of graphs of the square of the amplitude of oscillations of the vertical projection of cylinder axis vs the Galileo number for five density ratios.

Figure 12. Square of amplitude $A(n_z)$ of oscillation of the vertical projection of the cylinder axis as a function of the Galileo number. Here $L/d=2$. The symbols are listed in table 7. The meaning of colours is given in the legend.

For moderate density ratios (up to 2), the typical weakly nonlinear behaviour typical for a supercritical bifurcation arising from a steady state was evidenced. In this case, the square of the instability amplitude is expected to grow proportionally to ${Ga-Ga_{crit}}$. This behaviour allowed us to obtain the critical Galileo numbers both by interpolation of growth rates and by extrapolation of the square of the amplitude with a very good agreement (see table 8). As can be seen in figure 12, the linear growth of the square of the amplitude of oscillations holds on a very large interval of Galileo numbers. The last column of table 8 contains the corresponding slopes. They decrease from $\rho _s/\rho =0$ to 2.

Table 8. Critical Galileo numbers of onset of vertical oscillations of the cylinder axis. Here $L/d=2$. Base and new states: state from which the oscillations develop and the arising state. $Ga_{crit,1}$, value of the critical Galileo number obtained by interpolation of growth rates at Galileo numbers close to the threshold; $Ga_{crit,2}$, value obtained by extrapolation to zero of the trend of the square of amplitude; Slope $A^2$, coefficient $K$ of the linear fit $A(n_z)^2=K\,(Ga-Ga_{crit})$ in the interval of linear growth.

At higher density ratios (5 and 10), the inertia delays the onset of the flutter which opens the way to fluid modes. For the length to diameter ratio $L/d=2$, the threshold of the primary bifurcation leading to the symmetric periodic state (in which the direction of the cylinder axis does not move and only the vertical velocity oscillates) lies at about $Ga=205$. For the density ratio of 5, the oscillations of the cylinder axis were detected at $Ga=215$ leading to the asymmetric quasi-periodic regime, in which the cylinder centre oscillates both parallel and perpendicularly to its axis. The asymptotic state is reached progressively, evolving over a symmetric periodic and crosswise oscillating state. This is illustrated in figure 13 for a simulation starting with a cylinder released at rest in an almost vertical position. The axis is deviated by $10^{\circ }$ with respect to the vertical direction in the $xOz$ plane. Figure 13(a) shows that the initial flutter of the cylinder practically decays at $t \sim 600$, where, nonetheless, oscillations of the vertical component of the velocity ($u_z$) corresponding to the symmetric periodic state subsist (not represented). Later on, a new instability breaking the initial planar symmetry and yielding a crosswise oscillating state sets in. The cylinder centre starts to oscillate perpendicularly to the initial trajectory plane as seen in figure 13(b). This crosswise oscillating state is, however, itself unstable and excites vertical oscillations. In the process, a subharmonic modulation of the crosswise velocity develops. In the asymptotic regime, the cylinder axis oscillates both in the vertical and horizontal direction (figure 13c) and the trajectory is fully three-dimensional. At the density ratio 5, similar asymptotic states were evidenced up to $Ga=250$. In two cases, the asymmetry is non-existent ($Ga=230$ and 233), so that we classified these regimes as weak quasi-periodic flutter. This particularity had, however, practically no influence on the trend of growth of the oscillation amplitude. At the largest considered Galileo number (300) the behaviour was found chaotic.

Figure 13. (a) Vertical oscillation of the cylinder axis. (b) Horizontal velocity of the cylinder centre perpendicular to the initial trajectory plane. (c) Horizontal (blue line) and vertical (red line) oscillation of the cylinder axis in the asymptotic regime. Here $L/d=2$, $\rho _s/\rho =5$ and $Ga=220$.

The density ratio 10 was investigated in less detail concerning the onset of fluid modes because the latter are virtually independent of the density ratio. At $Ga \ge 250$ two new features were found. First, the quasi-periodic regime gives way to a chaotic state with a very small amplitude of oscillations. The latter, measured as standard deviation of $n_z$ from zero, is typically about 0.05 (less than $3^{\circ }$). The trajectory of the cylinder centre is fully chaotic and the orientation of the cylinder axis strongly drifts horizontally. (At $Ga=280$, we evidenced a horizontal rotation of the cylinder axis by 23$^{\circ }$ within 50 time units, so that it can be said that the horizontal drift largely dominates the vertical oscillations of the cylinder axis.) Second, two simulations (at $Ga=280$ and 290) attest a coexistence of the chaotic state with a strong flutter. The stability interval of this branch seems to be short, since simulations restarted at $Ga=270$ and 300 decayed to the chaotic regime with small amplitude. However, the existence of a stable strong flutter agrees with the trend found for aspect ratio 3, where subcritical bifurcations with branches of strong flutter were clearly evidenced.

4.2. $L/d=3$

In the investigation of cylinders of aspect ratio $L/d=3$, the density ratio 1 has been replaced by 0.5 and 1.16, the latter being that of cylinders used in experiments by Toupoint et al. (Reference Toupoint, Ern and Roig2019). The difference between $\rho _s/\rho =1$ and $\rho _s/\rho =1.16$ is, however, qualitatively inessential. Moreover, the experimental data (figure 4 of Toupoint et al. Reference Toupoint, Ern and Roig2019) is very coarse around the threshold of the loss of stability of the steady vertical fall located between $Ar=200$ and 400, i.e. $Ga= 145$ and 289. The value in table 9 falls well within this interval but the asymptotic state at $Ar=400$ was found ‘irregular’ and only at $Ar=500$ ($Ga=362$) a ‘fluttering’ state was reported.

Table 9. Critical Galileo numbers of the onset of vertical oscillations of the cylinder axis. The meaning of the columns is the same as in table 8 except for the subcritical branches where the lower limit of stability is indicated in the column labelled $Ga_{crit,1}$. Here $L/d=3$.

Figure 14 shows many common features with the cylinders of aspect ratio 2. For density ratios $\rho _s/\rho \le 2$ the primary bifurcation is, again, supercritical. The steady vertical states loses its stability in favour of the flutter. The square of the amplitude of oscillations grows linearly on a large interval of Galileo numbers. Unlike for $L/d=2$, the slopes decrease only very little.

Figure 14. (a) Amplitude $A(n_z)$ of oscillation of the vertical projection of the cylinder axis. The filled circles denote flutter starting to lose its perfect periodicity and symmetry in an initial stage of transition to chaos. (b) Enlarged view of the instability onset represented by the square of amplitude. Here $L/d=3$.

The most significant difference as compared with $L/d=2$ consists of the clear coexistence of two branches of asymptotic solutions evidenced at density ratios 5 and 10. The lower stability limits of branches of strong flutter could be determined quite reliably from the transients of the convergence to the attractor. It could be traced until the upper limit of investigation at $Ga=300$ and they can be expected to continue further. The branches of small oscillations arising from fluid modes can be expected to have upper limits of stability due to subcritical bifurcations. The solution then settles to the branch of strong flutter. For the density ratio 10, the threshold of this bifurcation has not been found within the limit of investigation, at $\rho _s/\rho =5$ it is situated between $Ga=250$ and 300. Figure 15(a) illustrates the loss of stability of the branch of strong flutter with decreasing $Ga$ and figure 15(b) illustrates the loss of stability of the asymmetric quasi-periodic branch with increasing $Ga$, both for density ratio 5.

Figure 15. (a) $Ga=195$: simulation restarted from the solution on the branch of strong flutter obtained at $Ga=200$. (b) $Ga=300$: simulation restarted from the solution on the quasi-periodic asymmetric branch at $Ga=250$. Here $L/d=3$ and $\rho _s/\rho =5$.

4.3. $L/d=5$

For cylinders of length-to-diameter ratio of 5, the scenario appears to be practically the opposite of that found for $L/d=2$ and 3. Small chaotic oscillations dominate for small instead of high density ratios and the standard supercritical bifurcation leading from the steady vertical fall directly to the flutter is only found at the highest considered density ratio of 10.

A closer look at figure 16 shows, however, that the fluid modes have been evidenced for all density ratios. For the density ratio 0, no flutter has been found in the investigated interval of Galileo numbers. For the density ratio 1, our simulations yielded a single case of flutter at the largest considered Galileo number of 450. This is in agreement with figure 4 of Toupoint et al. (Reference Toupoint, Ern and Roig2019) where ‘irregular’ behaviour was reported up to $Ar=600$ ($Ga=434$). In what follows, we first describe the density-independent scenario of fluid modes before focusing on states with vertical oscillations of the cylinder axis where the density dependence requires a more detailed descriptions as a function of the density ratio. Table 10 sums up the critical Galileo numbers delimiting the stability intervals discussed in the text.

Figure 16. Square of amplitude $A(n_z)^2$ of oscillation of the vertical projection of the cylinder axis. Here $L/d=5$. (a) Two-dimensional superimposition for all density ratios. (b) Expanded view of the amplitude $A(n_z)$ showing the individual states along the branches of small oscillations.

Table 10. Critical Galileo numbers of regimes of the cylinder of aspect ratio $L/d=5$: $Ga_1$, lower limit of stability; $Ga_2$, upper limit of stability.

4.3.5. $\rho _s/\rho =5$

A similar subcritical bifurcation sets in, but already in the steady vertical state at $Ga=172.5$. The subcritical branch of flutter remains stable down to $Ga=146.7$ and very likely far beyond the upper limit of investigation at $Ga=300$. This yields a bi-stability interval $146.7 \le Ga \le 172.5$. Another bi-stability has been evidenced for $Ga > 220$, i.e. at the onset of the crosswise oscillating state, which appears to be stable at this Galileo number. Starting from this threshold the branch of small oscillations follows the scenario of transition to chaos characteristic of fluid modes. Let us remark that the two last chaotic states at $Ga=300$ and 350 present random fluctuations that, in the long term, may happen to be strong enough to let them settle on the flutter branch. This can occur after a very long simulation time as can be seen in figure 17 where the cylinder rotates horizontally and weakly oscillates vertically before toggling to the flutter.

Figure 17. Azimuth $\theta$ and elevation $\psi$ in radians of the cylinder axis of a cylinder released in an initial position $\theta =0$, $\psi =0.087$ ($5^{\circ }$). Here $L/d=5$, $\rho _s/\rho =5$ and $Ga=350$.

5.1. Drag coefficient

The commonly used definition of the drag coefficient for oblong cylinders uses the constant projection $Ld$ normal to the axis (Clift et al. Reference Clift, Grace and Weber1978) as the section perpendicular to the flow even if the cylinders fall freely and oscillate. For freely falling cylinders, the balance between the hydrodynamic force, the gravity and buoyancy relates the drag coefficient $C_d$ to the vertical velocity. In dimensional units (Toupoint et al. Reference Toupoint, Ern and Roig2019),

(5.2)\begin{equation} C_d= \frac{\rm \pi}{2} \frac{g_{eff} d}{\overline{u_z}^2},\end{equation}

where $\overline {u_z}$ is the average (dimensional) vertical velocity. Given our non-dimensionalisation of the velocity by the velocity scale (2.1), the formula (5.2) becomes

(5.3)\begin{equation} C_d= \frac{2 d}{L} \frac{1}{\langle u_z\rangle^2}, \end{equation}

where $\langle u_z\rangle$ is the average non-dimensionalised vertical velocity. In spite of the definition accounting for the cylinder length, it cannot be expected that the drag coefficient is independent of the aspect ratio, especially for moderate aspect ratios. This effect is already accounted for in table 6.1 of Clift et al. (Reference Clift, Grace and Weber1978). For this reason, we investigate each of the three aspect ratios $L/d=2,3,5$, separately.

Let us present the methodology for the length to diameter ratio of 5.

Results of 81 simulations are available in this case. Figure 18(a) shows that values of the drag coefficient in regimes with small or non-existent vertical oscillations are distributed along a single line, which represents a master curve for non-oscillating regimes. The larger the amplitude (represented by the colour bar), the larger the deviation from this curve. We shall set a lower limit for the amplitude of flutter at $A(n_z)= A_{min}=0.15$. Out of the original 81 values, 54 lie below this limit. Lowering it from 0.15 to 0.05 reduces the number of ‘non-fluttering’ or ‘weakly fluttering’ regimes to 50, which does not significantly change the following analysis.

Figure 18. (a) Scatter plot of drag coefficient vs Reynolds number coloured as a function of the vertical oscillation amplitude $A(n_z)$. The marker types indicate the density ratio: $\rho _s/\rho =0$, triangle; $\rho _s/\rho =1$, square; $\rho _s/\rho =2$, diamond; $\rho _s/\rho =5$, cross; $\rho _s/\rho =10$, circle. The colours of the colour bar quantify the oscillation amplitude. Full line: fit (5.4) of non oscillating or weakly oscillating data. (b) Drag correction with respect to the fit in (a) as a function of the square of oscillation amplitude (for colours, see the legend). Dotted lines: fit to the function 5.6. Here $L/d=5$.

The regimes with amplitudes less than $A_{min}$ will be processed as ‘practically’ steady $A(n_z) \approx 0$. The aspect ratio having been fixed and steady states being independent of the density ratio, there remains a single variable function to be fitted to $Ga$ or, almost equivalently, $Re$. (Since $Re$ is related to $Ga$ by the vertical velocity, itself a dependent variable subject to variations, the fits are not exactly the same.)

Clift et al. (Reference Clift, Grace and Weber1978) provided correlations for cylinders of aspect ratios 1 through $\infty$ (in table 6.1 of Clift et al. Reference Clift, Grace and Weber1978) as well as curves of drag coefficient vs Reynolds number in the steady regime at Reynolds numbers up to about 100. Their correlation results from a third degree polynomial fit of the logarithm of the Reynolds number vs the variable $w=\tfrac {1}{3}\log _{10}(C_d \,Re^2)$. The coefficients are given in table 6.1 of Clift et al. (Reference Clift, Grace and Weber1978) as functions of the aspect ratio $E=L/d$. In what follows, we use a simple second degree fit to the ensemble of drag coefficient values of weakly oscillating regimes:

(5.4)\begin{equation} \log(C_{d,0})= a_2 \log^2(Re) + a_1 \log(Re) + a_0. \end{equation}

Figure 19 presents the comparison of three fits: (i) $C_{d,Clift}$ of Clift et al. (Reference Clift, Grace and Weber1978, table 6.1); (ii) $C_{d,fC}$ obtained by the method of Clift et al. (Reference Clift, Grace and Weber1978) applied to our data; (iii) $C_{d,0}$ given by (5.4). As can be seen, the difference between $C_{d,fC}$ and $C_{d,0}$ is of the order of $10^{-3}$ in the interval $50 \le Re \le 300$ and can be neglected. The original fit $C_{d,Clift}$ by Clift et al. (Reference Clift, Grace and Weber1978) agrees for $Re<50$ and overestimates the drag coefficient for higher Reynolds number values. The coefficients of the fit (5.4) are provided in table 11. The standard deviation from the fitted master curve is also presented and is less than 0.005.

Figure 19. Comparison of the correlation from table 6.1 of Clift et al. (Reference Clift, Grace and Weber1978) ($C_{d,Clift}$) and of the fits of the computed drag coefficients for regimes where $A(n_z)<0.15$ using the method of Clift et al. (Reference Clift, Grace and Weber1978) ($C_{d,fC}$) and the second-degree formula (5.4) ($C_{d,0}$). The graphs of $C_{d,0}$ and $C_{d,fC}$ are practically superimposed. The enlarged difference $C_{d,0}-C_{d,fC}$ is seen in the inset.

Table 11. Coefficients of the fits (5.4) and (5.7).

The effect of the density ratio can be assessed from the marker types in figure 18(a). The largest upward deviation from $C_{d,0}$ is observed for the density ratio 2 (diamonds), whereas the largest downward shift is obtained for the density ratio 10 (‘o’). At the density ratio 5, the drag coefficient remains close to the fitted curve. Qualitatively, it is clear that the deviation from the master curve depends on the oscillation amplitude and the density ratio.

It may, strictly speaking, vary also as a function of the Reynolds (Galileo) number. However, if the difference of the drag coefficient from the master curve $dC_{d}= C_{d}-C_{d,0}$ is plotted vs the square of the amplitude (figure 18b), it is seen that $dC_{d}$ is quite well proportional to the square of amplitude of oscillations for a given density ratio independently of the Reynolds (Galileo) number. A simple polynomial fit as a function of the density ratio on the interval $[0,10]$ yields a satisfactory result but, since the drag coefficient remains finite at both ends of the interval $[0,\infty ]$ of density ratios, a transformation to a finite interval of the independent variable is more inadequate and provides the possibility of extrapolation to infinite density ratios. We shall use the transformation $[0,\infty ] \rightarrow [0,1]$:

(5.5)\begin{equation} s = \frac{r}{r_0 +r}, \end{equation}

where $r$ denotes the density ratio $r \equiv \rho _s/\rho \in [0,\infty [$ and $s \in [0,1]$ is the new variable used for the regression. The parameter $r_0$ makes the distinction between ‘light’ and ‘dense’. We consider a global two-dimensional regression, with respect to both variables $s$ and $A^2(n_z)$ and linear in each of them:

(5.6)\begin{equation} C_d - C_{d,0} \approx dC_d = \left(c_1 s + c_2 \right) A^2(n_z)+ c_3. \end{equation}

The regression uses the whole sample (of 81 sets of simulation results for $L/d=5$). The value of $r_0$ can be used for minimising the error and constitutes a fourth parameter of the fit. Increasing the degree of the regression brings no significant improvement and tends to yield spurious oscillations. Another method to improve the fit is to take account of the Reynolds number variation. This brings only a 10 % reduction of the (already small) error. The weak effect of the Reynolds number is in agreement with the general observation that the drag is almost independent of the Reynolds number starting from $Re=150$ (see Toupoint et al. Reference Toupoint, Ern and Roig2019, figure 7). As a result, the simple bi-linear fit (5.6) is satisfactory. Let us remark, that the amplitude of oscillations of the cylinder axis can be replaced by the amplitude of oscillations of the horizontal velocity because both are tightly correlated.

Figure 20 presents the slopes of the $A^2(n_z)$ scalings as a function of the density ratio for all the three length-to-diameter ratios of cylinders. Note that, for $L/d=5$ the slopes become negative. For a comparison, a linear (one-dimensional) fit has been performed individually for each density ratio. The so-obtained slopes represent the exact trends uncorrelated with other density ratios. They are also represented in the figure. The coefficients $c_i$, $i=1,2,3$, of the least-squares fit as well as the optimal value of $r_0$ are provided in table 12. The small value of the coefficient $c_3$ correctly accounts of the fact that the lines in figure 18(b) should start at zero. Table 12 also contains the total error cumulating the errors of both fits (5.4) and (5.6). As can be seen, the difference between the computed and the fitted values is so small that it makes no sense to present a figure of the scatter of the fitted values to compare them with figure 18(a).

Figure 20. Fitted slopes $(c_1 s + c_2 )$ of the $A^2(n_z)$ scalings (5.6). Dashed lines, fit vs $Re$; dash-dotted lines, fit vs $Ga$. The colours correspond to the aspect ratio (see legend). Extrapolations to infinite density ratio (values in brackets: fit vs $Ga$): $L/d=2: -0.09 (-0.16)$, $L/d=3: -0.18 (-0.20)$, $L/d=5: -0.44 (-0.53)$. The markers indicate the slopes obtained separately for each density ratio by a one-dimensional linear fit. Circles (triangles) refer to regression with the Reynolds (Galileo) number as the independent variable.

Table 12. Coefficients of the fits (5.6) and (5.8). $^a$The minimum error is reached for a value larger than the maximal density ratio 10 (18.8 and 14.4, respectively). Here 10 is used with a negligible loss of accuracy.

The correlations of the drag coefficient as a function of the Reynolds number are most common in the available literature. However, unlike the Galileo number, the Reynolds number is not an external parameter of the problem. The scatter of the vertical velocity in (5.1) influences the correlation (5.4) and thus also the drag correction (5.6). In tables 11 and 12 we also present the result of corresponding fits using the Galileo number as the independent variable:

(5.7)\begin{gather} \log(C_{d,0,Ga})= a_2 \log^2(Ga) + a_1 \log(Ga) + a_0, \end{gather}

(5.8)\begin{gather} C_d - C_{d,0,Ga} \approx dC_{d,Ga} = \left(c_1 s+ c_2\right) A^2(n_z)+ c_3. \end{gather}

Though the difference of the coefficients $c_i$ in (5.6) and (5.8) seems to be non-negligible, the difference of resulting fits $C_{d,fit,Re}=C_{d,0}+dC_{d}$ and $C_{d,fit,Ga}=C_{d,0,Ga}+dC_{d,Ga}$ presents a smaller root mean square than the total error of the drag coefficient.

For $L/d=3$ and 2, the processing follows exactly the same lines. A set of 114 simulation results has been processed for $L/d=3$. Figure 21(a) shows that the oscillations systematically increase the drag. As shown in figure 21(b) the drag correction still scales as the square of the amplitude. This is not necessarily expected for density ratios 5 and 10. For these density ratios the branches of large flutter are subcritical and disconnected from branches of small oscillations (see figure 14a) and the scaling cannot be substantiated by weakly nonlinear arguments. In spite of that, the $A^2(n_z)$ scaling still approximately holds. For example, at $\rho _s/\rho =5$, the five points corresponding to large values of amplitude and those of the supercritical branch of small oscillations are well aligned along the same straight line. Similarly as for $L/d=5$ the slopes decrease with increasing density ratio without, however, becoming negative. The overall fit (including that of the drag coefficient without oscillations (5.4) and that of the drag correction due to oscillation (5.6) yields an error as small as 0.0046 (see tables 11 and 12 for the values of coefficients and the errors).

Figure 21. (a) Scatter plot of drag coefficient vs Reynolds number coloured as a function of the vertical oscillation amplitude $A(n_z)$. Full line: fit (5.4) of non-oscillating or weakly oscillating data. The marker types indicate the density ratio (see the caption of figure 18, except for density ratios 0.5 (squares) and 1.16 (asterisks)). (b) Drag correction with respect to the fit in figure (a) as a function of the square of oscillation amplitude. Dotted lines: fit to the function (5.6). Here $L/d=3$.

The case of $L/d=2$ has been investigated in less detail. Although 52 simulations are available for the processing only 20 satisfy the criterion $A(n_z)>0.15$ (see figure 22). Qualitatively, there is much less difference between $L/d=2$ and 3 than between 3 and 5. The fit results are, again, presented in tables 11 and 12. For $\rho _s/\rho =0$, the drag correction values at large amplitude are not quite aligned. Obviously, the Galileo number dependence starts to play a role. This is due to the long array of simulations at $\rho _s/\rho =0$ in figure 12 spanning a long interval $Ga \in [115,200]$. In spite of that, the final error of the fit is smaller than 1 %.

Figure 22. (a) Scatter plot of drag coefficient vs Reynolds number coloured as a function of the vertical oscillation amplitude $A(n_z)$. Full line: fit (5.4) of non-oscillating or weakly oscillating data. For the meaning of the marker types refer to caption of figure 18. (b) Drag correction with respect to the fit in (a) as function of the square of oscillation amplitude. Dotted lines: fit to the function 5.6. Here $L/d=2$.

All other quantities vary similarly with the growth of the square of amplitude. This is the case of, for example, the vertical velocity, the Reynolds number and the oscillation period. The latter is analysed using a similar method in the next subsection.

5.2. Strouhal number

In this section, we focus on the frequency of oscillations of the cylinder axis in the vertical direction. As a consequence, we discard the frequencies of the fluid modes where the direction of the axis remains horizontal (symmetric periodic and crosswise oscillating states). Though not all regimes, where vertical oscillations of the cylinder axis are observed, are pure flutter, the frequencies do not differ significantly except for the case of the weak quasi-periodic flutter for the aspect ratio 3 and density ratio 10 described in § 3.2.2. These frequencies are not considered either. This still leaves 231 out of the total of 247 sets of simulation results for the analysis.

In the available literature, there is no general consensus as for the non-dimensionalisation of the frequency of falling and fluttering oblong cylinders. Marchildon et al. (Reference Marchildon, Clamen and Gauvin1964) remarked that, unlike for the (fixed) infinite cylinder and its von Kármán vortex street, the cylinder diameter, as length scale in the definition of the Strouhal number, is to be replaced by $\sqrt {Ld}$. For the velocity scale the (mean) vertical velocity of the fall is used. This yields

(5.9)\begin{equation} St = \frac{f \sqrt{Ld}}{\overline{u_z}} = \sqrt{\frac{L}{d}} \frac{1}{T \langle u_z\rangle}, \end{equation}

where the first expression is written in dimensional units and the second in the non-dimensionalisation of the present paper. Alternatively, Chow & Adams (Reference Chow and Adams2011) and Toupoint et al. (Reference Toupoint, Ern and Roig2019) use the gravitational velocity (1.2) as velocity scale. In the paper by Toupoint et al. (Reference Toupoint, Ern and Roig2019), this results in the non-dimensionalisation

(5.10)\begin{equation} St^* = \frac{f \sqrt{Ld}}{u_{g,eff}} = \sqrt{\frac{\rm \pi}{4}} \frac{L}{d} \frac{1}{T}, \end{equation}

where the second expression relates this Strouhal number to the inverse of period expressed in the non-dimensionalisation of the present paper. Chow & Adams (Reference Chow and Adams2011) provided an account, moreover, of the effect of density ratio and suggest that the Strouhal number

(5.11)\begin{equation} St_{CA} = \sqrt{\frac{\rho_s}{\rho}} St^*\end{equation}

should assume a constant value of $1/7.91$ (0.126). Figure 23 compares the three different definitions of Strouhal numbers (5.9), (5.10) and (5.11) in light of our simulation data.

Figure 23. Scatter plots of the values of (a) $St$ (5.9), (b), $St^*$ (5.10) and (c) of $St_{CA}$ (5.11). The colour bar indicates the density ratio, the marker type the aspect ratio: o, $L/d=2$; x, $L/d=3$; $\cdot$, $L/d=5$.

The Strouhal numbers $St$ and $St^*$ yield similar plots because

(5.12)\begin{equation} St^*/St = \overline{u_z}/u_{g,eff},\end{equation}

i.e. the ratio $St^*/St$ is equal to the mean vertical velocity non-dimensionalised by the gravitational velocity scale. This ratio is situated between 1.2 and 1.5 at Galileo numbers $150 \le Ga \le 300$ for all aspect and density ratios. The value assumed by Chow & Adams (Reference Chow and Adams2011) is 1.33 and data of Toupoint et al. (Reference Toupoint, Ern and Roig2019) also agree with these values. The factor $\sqrt {\rho _s/\rho }$ suggested by Chow & Adams (Reference Chow and Adams2011) over-compensates the effect of the density ratio (and makes no sense for very light bodies). Here $St^*$ concentrates the best the values for all the three aspect ratios for density ratios $\rho _s/\rho \le 2$. The average of 0.1 is in agreement with that of Toupoint et al. (Reference Toupoint, Ern and Roig2019). For larger density ratios (5 and 10), the values are considerably scattered both due to different aspect and density ratios. For this reason, we shall proceed similarly as in § 5.1 and treat each of the aspect ratios separately. Similarly as for the regression of the drag coefficient vs the Galileo number and the Reynolds number, it is practically indifferent whether the analysis is performed for $St$ or $St^*$ and whether the Reynolds or the Galileo number is taken as the independent variable. We shall choose $Ga$ as an independent variable and $St^*$ as a dependent variable.

Figure 24 presents scatter plots of $St^*$ vs $Ga$ separately for each aspect ratio $L/d=2$, $3$ and 5. Markers identify the density ratios, the colouring indicates the amplitude of oscillations. The inertia plays a significant role for density ratios 5 (crosses) and 10 (circles) yielding significantly lower frequencies than for light bodies. The bi-stabilities emphasise the important role of the amplitude tending to lower the frequency (increase the period) in agreement with expectation namely if tumbling is approached for very large amplitude.

Figure 24. Scatter plots of the values of $St^*$ (5.10) vs $Ga$ for the three aspect ratios: (a) $L/d=2$; (b) $L/d=3$; (c) $L/d=5$. The density ratios are identified by symbols of the legends, the amplitudes of oscillations are indicated by the colour bars. (d) $L/d=5$, fit: circles, exact values of $L/d=5$; crosses, fit (5.13) with coefficients given in table 13.

To obtain an accurate fit, it was seen in the case of the drag coefficient, that higher degree regression was needed. Moreover, the regression procedure could be split in two stages: a single-parameter fit of the master curve for negligible oscillations (independent of the density ratio) and a subsequent two-parameter regression accounting for oscillations. In contrast, the frequency depends on the density ratio whatever the oscillation amplitude making such a two step approach of little interest. To avoid too complicated and statistically unreliable correlations, we limit the regression to a tri-linear (linear in all of the three variables) fit as a function of $G \equiv Ga$, $r \equiv \rho _s/\rho$ and $a^2 \equiv A^2(n_z)$. The main source of inaccuracy comes from the dependence on the density ratio which is not even monotonous, however we consider that more data are needed to account for this behaviour reliably. Testing of the error has shown that the transformation (5.5) yields monotonously decreasing error with increasing $r_0$, which means that a regression depending linearly directly on the density ratio $r$ performs the best. As a consequence, for each of the three aspect ratios, we determine the least-squares fit of

(5.13)\begin{equation} St^* \approx St^*_{fit} = c_1 r a^2 G +c_2 r G + c_3 r a^2 + c_4 r + c_5 a^2 G + c_6 G + c_7 a^2 + c_8 \end{equation}

and compute the relative root-mean-square error

(5.14)\begin{equation} \varepsilon = \frac{\sqrt{\sum_i (St^*_i - St^*_{fit,i})^2}}{\sqrt{\sum_i (St^*_i)^2}} \end{equation}

and maximum relative error

(5.15)\begin{equation} \varepsilon_{max} = \mbox{max} (|St^*_i - St^*_{fit,i}|/St^*_i). \end{equation}

The numerical values of the coefficients $c_i$ of the fit (5.13) are given in table 13. The somewhat over-simplified tri-linear fit yields considerably larger errors than those obtained for the drag coefficient in table 12. Nevertheless, the several percent inaccuracy may be acceptable as a comprehensive summary of the obtained data. The worst case of aspect ratio $L/d=5$ is, moreover, represented graphically in figure 24 (scatter plot $L/d=5$, fit). The extreme values (corresponding to $\rho _s/\rho =0$ and 10) are well reproduced, the inaccuracies concern the intermediate density ratios.

Table 13. Coefficients of the fits (5.13) for the Strouhal numbers (5.10): $\varepsilon$, relative root-mean-square error (5.14); $\varepsilon _{max}$, maximum relative error.