Numerical analyses of the flow past a short rotating cylinder

Abstract This work studies the three-dimensional flow dynamics around a rotating circular cylinder of finite length, whose axis is positioned perpendicular to the streamwise direction. Direct numerical simulations and global stability analyses are performed within a parameter range of Reynolds number $Re=DU_\infty /\nu <500$ (based on cylinder diameter $D$, uniform incoming flow velocity $U_\infty$), length-to-diameter ratio ${\small \text{AR}}=L/D\leq 2$ and dimensionless rotation rate $\alpha =D\varOmega /2U_\infty \leq 2$ (where $\varOmega$ is rotation rate). By solving Navier–Stokes equations, we investigated the wake patterns and explored the phase diagrams of the lift and drag coefficients. For a cylinder with ${\small \text{AR}}=1$, we found that when the rotation effect is weak ($0\leq \alpha \lesssim 0.3$), the wake pattern is similar to the unsteady wake past the non-rotating finite-length cylinder, but with a new linear unstable mode competing to dominate the saturation state of the wake. The flow becomes stable for $0.3\lesssim \alpha \lesssim 0.9$ when $Re<360$. When the rotation effect is strong ($\alpha \gtrsim 0.9$), new low-frequency wake patterns with stronger oscillations emerge. Generally, the rotation effect first slightly decreases and then sharply increases the $Re$ threshold of the flow instability when $\alpha$ is relatively small, but significantly decreases the threshold at high $\alpha$ ($0.9<\alpha \leq 2$). Furthermore, the stability analyses based on the time-averaged flows and on the steady solutions demonstrate the existence of multiple unstable modes undergoing Hopf bifurcation, greatly influenced by the rotation effect. The shapes of these global eigenmodes are presented and compared, as well as their structural sensitivity, visualising the flow region important for the disturbance development with rotation. This research contributes to our understanding of the complex bluff-body wake dynamics past this critical configuration.


Introduction
The study of flows around rotating bluff bodies, including cylinders and spheres, constitutes a fundamental problem in fluid dynamics.It provides critical insights into vortex formation and wake dynamics, which are relevant to various natural phenomena and engineering applications.For instance, these flows are critical in circulation control on an airfoil (Tennant et al. 1976), heat exchangers (Roslan et al. 2012), laminar/turbulence separation (Afroz et al. 2017) and design of guided rockets (de Celis et al. 2017).Besides, the use of rotating effects to control the wake flow past bluff bodies has attracted much attention (Gad-el Hak & Bushnell 1991;Modi 1997), including its applications such as Flettner rotors utilising the Magnus effect (Seifert 2012).However, there is a research gap in the study of flow past a rotating cylinder with two free ends.Experimental studies and real-world applications typically involve finite-length rotating cylinders.In contrast, theoretical and numerical studies traditionally consider infinitely long rotating cylinders.This discrepancy highlights the need for further investigation on the effect of aspect ratio and free ends.To address this research gap, our study aims to examine the threedimensional (3-D) flow past a short rotating cylinder, an area that previous researchers have not explored.By investigating this problem, we hope to contribute to the knowledge base on the flow around rotating bluff bodies and provide insights that are relevant to various engineering applications.
For a non-rotating infinitely long circular cylinder, it is well known that the wake first experiences a Hopf bifurcation at Re ≈ 46.7 and then a 3-D wake transition at Re ≈ 190 (Williamson 1996a,b).The non-dimensional number Re quantifies the ratio between inertia and viscosity.The rotation effect enriches the flow dynamics of the wake flow.Kang et al. (1999) conducted a two-dimensional (2-D) numerical study of the flow past a rotating circular cylinder and showed that rotation could effectively suppress the vortex shedding (mode I) found in the stationary cylinder at α > α c (where α c is the critical dimensionless rotation rate), and the relationship between the lift/drag and rotation rate in the range of 0 ⩽ α ⩽ 2.5 is significantly different from that predicted by potential flow theory.Subsequently, the numerical studies of 2-D rotating cylinders by Stojković et al. (2002Stojković et al. ( , 2003)), Mittal & Kumar (2003) and Mittal (2004) revealed that when the rotation rate α is relatively large, there is a secondary instability phenomenon (model II) characterized by low-frequency vortex shedding.Especially, Mittal & Kumar (2003) brought to light this instability mechanism of 2-D perturbations by the global stability analysis, which will be extended in the current study, focusing on a rotating finite cylinder, to account for 3-D perturbations.Built upon the previous works, El Akoury et al. (2008) extended the neutral stability curves for these wakes in the Re-α plane by direct numerical simulations (DNS) and Landau model.The experiments by Kumar et al. (2011) provided evidence of the existence of mode II at Re = 200, 300, 400 and 0 < α < 5.The experimental study of Linh (2011) also reported observations of the lowfrequency mode II vortex.Their experimental Strouhal number and wake patterns agree well with numerical data of Mittal & Kumar (2003).Later, Pralits et al. (2010Pralits et al. ( , 2013) ) conducted an extensive study of the linear global dynamics of the 2-D rotating cylindrical wake flow.The authors explored neutral stability curves on the (Re, α) plane, providing a comprehensive understanding of this phenomenon.They also observed multiple steady solutions at high α, explaining the decay of the secondary shedding wake.More recently, Sierra et al. (2020) fully described the bifurcation, neutral curves and global instability modes in the parameter space (Re, α) ⊂ [0, 200] × [0, 10], exploring the relations among Takens-Bogdanov bifurcations, cusps and generalized Hopf bifurcations when varying the parameters in the rotating cylinder wake flow.To sum up, for the infinitely-long rotating cylinders, mode I and mode II are fundamentally different flow phenomena.The mode I undergoes a supercritical Hopf bifurcation and becomes linearly unstable at 0 ⩽ α ⩽ 2, which is related to the classical Bénard-von-Kármán vortex street, characterized by alternating vortices with opposite signs of spanwise vorticity.On the other hand, the physical mechanism of the linear instability mode II (4.5 ⩽ α < 6) is featured by low-frequency vortex shedding with the same vorticity sign.
Further research has shown that rotation can lead to complex 3-D instabilities.Numerical investigations by Rao et al. (2013a,b) demonstrated several 3-D modes becoming unstable to spanwise perturbations in the steady and unsteady regimes of Re = 400 flows.Five 3-D modes were identified to be unstable in the mode I shedding regime, while four 3-D modes were observed in the steady flow regimes for α ⩾ 2. Radi et al. (2013) proved experimentally the existence of the above numerically predicted 3-D modes.They additionally showed a highly 3-D wake and the absence of 2-D periodic shedding at high α (i.e.Re = 200, α = 4.5) previously reported in Mittal & Kumar (2003).Navrose et al. (2015) conducted 3-D numerical studies and found that the span length of the rotating cylinder plays an important role in the evolution of the wake with Re ∈ [200,350], α ∈ [0, 5].Specifically, only linear global modes with wavelengths that are integer multiples of the cylinder span are selected for growth in nonlinear DNS.It is thus necessary to consider 3-D configurations that take into account the spanwise length in studies of rotating bluff-body flows.
Researchers have also studied the flow past a rotating bluff body of other forms.It is instructive to review relevant works on these flows as they will also be discussed in this study.The flow past a rotating sphere, along either the transverse axis or the streamwise axis, received considerable attention.Citro et al. (2016) applied the global linear stability analysis (LSA), adjoint-based structural sensitivity analysis and weakly nonlinear analysis (WNL) to reveal the mechanism of flow instability around a rotating sphere around the transverse axis.They characterised the evolution processes of the first (at low α) and second (at high α) instability modes.Fabre et al. (2017) further improved the WNL from than that in Fabre et al. (2012) to achieve a better comparison between the WNL expansion result and the DNS.Namely, the comparison demonstrates that Fabre et al. (2017)'s ϵ expansion (ϵ = √ Re − Re c /Re c , where Re c is the critical Reynolds number of a pitchfork bifurcation in their work) provides a better reproduction of the DNS results for both angular velocity and associated lift forces, compared to Fabre et al. (2012)'s ω expansion (where ω represents the dimensionless rotation rate normalizing the actual rotation with U 0 /D).The ω expansion failed to predict the DNS results for Re around and beyond Re = 212, whereas the ϵ-expansion accurately reproduces the DNS results up to Re ≈ 225.For the flow past a sphere rotating along the streamwise direction, Lorite-Díez & Jiménez-González (2020) conducted a DNS study on the wake evolution of a strongly rotating sphere and observed a sequence of continuous bifurcations from periodic, quasi-periodic, and irregular states to chaos, over the parameter range 0 < α < 3 and Re = 250, 500, 1000.Later, Sierra-Ausín et al. (2022) employed global LSA to determine the neutral curves of three non-zero frequency global modes on the Reα plane, and used the normal form expansion to reveal the nonlinear interactions among the global modes.Their predictions of the normal form analysis were satisfactory and close to the DNS results, which led to a more detailed phase diagram of the nonlinear patterns.Besides, Jiménez-González et al. (2014) carried out global LSA of the wake flow past a streamwise rotating bullet-shaped body and plotted the neutral curves on the Re-α plane.Their work indicated that the streamwise rotation can also delay the Hopf bifurcation of a bluff body with a large aspect ratio (ar = 2) when increasing Re, which is different from a sphere with an aspect ratio of 1.To sum up, in addition to the sphere rotating along the streamwise direction, the aforementioned rotating bluff bodies of various shapes and aspect ratios exhibit a moderate rotation regime in the Re − α plane where the neutral curve of the Hopf bifurcation is significantly shifted to higher Re values.Therefore, the aspect ratio plays a significant role in the wake transition of a rotating bluff body, which should be further researched.
Our literature review has identified a research gap in the understanding of the dynamics of uniform flow past a finite-length rotating cylinder along its cylinder axis.This flow configuration is common in nature and engineering applications, but its instability mechanism, bifurcation properties and transition path are still unclear.In a recent study, we have conducted a detailed investigation of the wake flow around a non-rotating finite-length cylinder (Yang et al. 2022).Building on this research, we aim to extend our investigation to the wake flow around a rotating finite-length cylinder to explore its 3-D effects, in line with other similar works focusing on the rotation effect in the wake flow such as Pralits et al. (2010); Citro et al. (2016); Sierra-Ausín et al. (2022); Zhao & Zhang (2023), reviewed in Rao et al. (2015).The primary objective of this study is to clarify the effect of rotation on the unsteadiness of finite-length cylinder wakes using the (nonlinear) DNS method.Additionally, we aim to determine the instability threshold at which the unsteadiness occurs using the global stability approach.To identify the instability region responsible for the unsteadiness, we will also probe the structural sensitivity of the flow.Our work will contribute to a deeper understanding of the wake flow around rotating finite-length cylinders and provide insights into the dynamics of this flow configuration, which has practical implications for various engineering applications.
The paper is organised as follows.Section § 2 introduces the configuration of a 3-D finite rotating cylinder flow, the boundary conditions, the governing equations (i.e.nonlinear Navier-Stokes equations and their corresponding linearised direct and adjoint equations) and the numerical methodology.In section § 3, we show the results and discuss the base states (time-averaged flow or steady flow), nonlinear wake patterns, global eigenmodes, neutral curves in parametric plane Re-α, bifurcations in this flow.Finally, the results are summarised in section § 4 and conclusions are provided.In the appendices, we provide additional results of the Stuart-Landau model, the global modes at different aspect ratios ar and a verification step of the numerical codes.

Flow configuration and governing equations
We study the 3-D stability of the flow around a finite-length rotating cylinder of length L, diameter D and aspect ratio ar = L/D, subjected to a uniform incoming flow in a Cartesian coordinate system.As shown in figure 1, the origin of the coordinate system is located at the center of the cylinder, the x axis points in the flow direction, the y axis represents the transverse direction and the z axis extends along the center line of the cylinder.The nondimensional Navier-Stokes (NS) equations for the unsteady Newtonian incompressible flow read where U = (U x , U y , U z ) is the velocity vector and P is the pressure.The Reynolds number Re = DU ∞ /ν is defined based on cylinder diameter D, the velocity of the uniform incoming flow U ∞ at infinity and the kinematic viscosity coefficient ν.Strouhal number St = f D/U ∞ is defined based on the frequency f of vortex shedding.The dimensionless rotation rate α = ΩD/2U ∞ , where Ω is the rotating angular speed of the cylinder along the z axis, as shown in figure 1. Setting ρU 2 ∞ as the reference dynamic pressure, the drag and lift coefficients are defined respectively as where ρ is the fluid density, F dp = Sc P x dS and F dv = Sc τ wx dS are the pressure drag and friction drag on the cylinder surface S c along the streamwise direction, F lp = Sc P y,z dS and F lv = Sc τ wy,wz dS are the pressure lift and wall shear stress lift acting in either y axis or n z (defined as C ly and C lz , respectively) and A is the reference area A = LD.Here (P x , P y , P z ) are the components of the pressure acting on the cylinder surface along the x, y and z axes, respectively.τ w is wall (surface) shear stresses.Furthermore,  x-axis, negative x-axis, negative x-axis and positive x-axis, respectively.The vector e b points in the direction that is perpendicular to both the normal vector and the tangent vector.Here I is the identity tensor.

Linearisation
The global linear stability/instability of the flows past the finite-length rotating cylinder will be studied.Reynolds decomposition U = U b + u, P = P b + p will be substituted into the nonlinear governing equations 2.1.The base-state terms (U b , P b ) satisfying the steady Navier-Stokes equations and the nonlinear terms are neglected, yielding the linearised equations for the infinitesimal perturbations (u, p) residing on these base states, i.e., where u is the 3-D perturbation velocity vector u = (u x , u y , u z ) and p is the perturbation pressure.Homogeneous boundary conditions are applied for the perturbed variables as follows, u = 0 on S c and S in , (2.5a) (2.5d) Linear equation 2.4 is rewritten in matrix form with q = (u, p) T as where L U b is the linearised Navier-Stokes operator depending on the base states U b .The elements of mass matrix M and the Jacobian matrix L U b are As the considered base flow states are steady, we seek the wavelike solution q(x, y, z, t) of the form q(x, y, z, t) = q(x, y, z)e λt , where λ = σ + i2πω.
(2.8) Substituting this form (Eq. 2.8) into equation 2.6, we can get the following eigenvalue problem where the stability of the base state U b is dictated by the eigenvalues λ in the linearised problem with σ being the temporal growth/decay rate of perturbations and ω the eigenfrequency.The flow is linearly unstable if σ > 0; stable otherwise.The eigenfrequency ω of the most unstable eigenvalue determines whether the base state U b experiences a regular bifurcation (ω = 0) or a Hopf bifurcation (ω > 0).Note that the flow problem considered in this work is not spatially periodic or homogeneous in either x, y, z directions, and q depends on all the three coordinates, leading to a global stability problem (Theofilis 2011).

Sensitivity analysis
Sensitivity analyses based on the adjoint approach (Luchini & Bottaro 2014) will be conducted to identify the instability mechanism responsible for the unsteadiness.Following Giannetti & Luchini (2007), the adjoint equations of the linearised Navier-Stokes equations read (2.10) where u + and p + are the adjoint vector of perturbation field u and p, respectively.Following Giannetti & Luchini (2007); Marquet et al. (2008); Citro et al. (2016), the boundary conditions of the adjoint equations are set as u + = 0 on S c , S in , S xz and S xy , (2.11a) The identification of the core region of the instability can help to understand the instability mechanism (Giannetti & Luchini 2007;Luchini & Bottaro 2014).According to Giannetti & Luchini (2007), the sensitivity wavemaker ζ can be identified by overlapping the direct eigenvector u and adjoint eigenvector u + , (2.12)

Numerical method
In order to obtain the accurate wake pattern and the base states of the fully 3-D flow past a short rotating cylinder at medium and low Reynolds numbers, we adopt the highorder parallelised open-source code Nek5000 (Fischer et al. 2020) (version 19.0), which is based on the nodal spectral element method (SEM) originally proposed by Patera (1984).
Hexahedral elements with a polynomial order N = 7 are used, which indicates that there are eight points in each spatial dimension of the element (Fischer et al. 2020).The time step ∆t is determined by the Courant-Friedrichs-Lewy (CFL) condition with the target Courant number ≲ 1.0.Following the practice in our previous work Yang et al. (2022), the boundary layer elements in the vicinity of the rotating cylinder have been refined by the O-type mesh (see figure 1b).The only difference lies in the enlarged computational domain and increased number of elements to accommodate the rotational effects.We focus on studying two types of base states (both of which are denoted as Q b = (U b , P b ) T in the text to follow as long as there is no confusion) to analyse their global dynamics.The flow will become unstable in a certain region of the parametric space.The unstable base flow Q b in this case cannot be obtained directly by time-evolving the NS equations.For these unstable base flows, the selective frequency damping (SFD) method proposed by Akervik et al. (2006) was used to obtain an equilibrium solution to the NS equations 2.1.This type of base state will be called (SFD) base flow.Another base state of interest is the mean flow, which is obtained by time averaging the periodic flow with vortex shedding.For the present global stability analysis, at least ten vortex shedding cycles will be used in the time-average procedure.
For a non-parallel 3-D flow past the finite rotating cylinder, the numerical discretisation of linearised NS equation 2.4 will result in a large-scale Jacobian matrix A in the generalised eigenvalue problem 2.9.It is impractical to solve a large-scale eigenvalue problem for its whole eigenspectrum in a 3-D flow.Based on Nek5000 solver and the ARPACK package (Lehoucq et al. 1998), the matrix-free time-stepper method (Theofilis 2011;Doedel & Tuckerman 2012) will be adopted in the present work, implementing the Implicitly Restarted Arnoldi Method (IRAM) (Radke 1996;Lehoucq et al. 1998).
The validation of the nonlinear DNS code and the linear stability code is provided in Appendix C together with a convergence study on the size of the computational domain.In this section, we present the base state of the flow past a short rotating cylinder.We focus on the case of ar = 1 and vary the values of Re and α.Both temporally mean flow and the steady base flow (solved using the SFD method if unstable) will be displayed and discussed.By extracting pressure field contours and streamlines, we show in figure 2 the typical unstable steady base flow obtained by the SFD method.First, we review the flow past a finite-length cylinder without rotation (α = 0) in panels (a, b).For Re < Re c = 172.2(figure 2a), the wake flow is steady, which has two mutually Besides, the positions of the two separation points slightly change, and the pressure contours indicate that the pressure near the upper part of the cylinder arc surface is larger than that of the lower part, which results in a negative lift coefficient C ly .The sign of the C ly in the case of no rotation is unimportant because of the symmetric setting.Figure 2(c, d) show respectively the base flow structures of cases Re = 160 (stable) and Re = 330 (unstable) at a low rotation rate α = 0.1.Even for a low Re = 160, it can be seen that the rotation breaks the symmetry of the wake, leaving the flow only symmetric with respect to the oxy plane.Compared with the non-rotation case (figure 2a), the positions of the upper and lower separation points shift significantly along the rotation direction.Overall, the asymmetric recirculation region generated by the weak rotation is similar to the asymmetric wake generated by regular bifurcation in the non-rotated case.The counterclockwise rotation makes the flow velocity near the cylinder arc surface on the lower side increase, and the velocity near the upper wall surface decrease.It can be deduced that the pressure near the lower side wall will decrease, and the pressure near the upper side wall will increase, so a negative lift C ly will be obtained, i.e. the Magnus effect.It can be seen from figures 2 (c − h) that the numerical results conform to the classical Magnus effect.In the parameter space studied in this paper (especially relatively low Re), there is no boundary layer transition from a laminar flow to turbulence, that is, there is no inverse Magnus effect (Kim et al. 2014).As mentioned above, the symmetry breaking effect caused by the regular bifurcation resembles that due to the rotation (see also the rotating sphere by Citro et al. 2016).Therefore, the rotation 'strengthens' the symmetry breaking caused by the inherent wake mechanism observed in the regular bifurcation without rotation, resulting in the lift force in the rotating cases in figures 2 (c − h) being larger, to be shown in figure 4(b).On the other hand, the drag coefficient C d is not very sensitive to the rotation when at low rotation rate α < 0.3, as shown in figure 4 (a).Continuing to increase the rotation rate, the figures 2(g) (Re = 160) and 2(h) (Re = 330) show the cases of a high rotation rate (α = 1.2), and their wake structure and pressure distribution are similar to those of figures 2(e) and 2(f ), but the wakes are swung further upwards.The hyperbolic stagnation point (see the white dot) on the upper side moves forward along the direction of tangential velocity, which is similar to the results of Citro et al. (2016) on a rotating sphere at a high rotation rate.

Separation bubble
Next, we discuss figure 3 where the streamlines are further analysed.Here, the separation point is defined as the intersection between the envelope of the separation bubble (indicated by the green curve with U x (x, y, z = 0) = 0 in figure 3a) and the surface of the cylinder, and thereby, the separation angle θ s is defined as the angle from the negative x-direction to the separation point.In the present rotating cylinder, the values of parameters Reynolds number, rotation rate, and aspect ratio have no effect on the separation angle.The separation bubble is always confined between angles θ s = 0 and θ s = π, as shown in figure 3(a).Note that in the case of non-rotation cylinder (Yang et al. 2022), the separation angles may be varying depending on Re and aspect ratio.We define the separation bubble length as . Panel 3(b) shows x s as a function of Re for the ar = 1 cylinder, compared with the separation length of a fixed sphere (Johnson & Patel 1999) and a non-rotating finite-length cylinder (Yang et al. 2022).The effect of Reynolds number on x s is different for the high-and low-rotation-rate cases.When the rotation rate is low, the increase of Reynolds number makes x s increase monotonically, which is similar to the (mean flow) cases of fixed finite-length cylinder (Yang et al. 2022) at 0.5 ⩽ ar ⩽ 2 prior to Hopf bifurcation, sphere (Johnson & Patel 1999), and rotating sphere (Sierra-Ausín et al. 2022).The increase of α decreases the value of x s at 0.025 ⩽ α ⩽ 1.2, which means that the rotation shortens the recirculation length.It has been proposed that shortening the recirculation length stabilises the flow; for example, Sierra-Ausín et al. (2022) discussed such a stabilising control strategy.
When the rate α is high, x s slightly decreases with the increase of Re, for example, see α = 0.6, 1.2 in figure 3(b).According to the above discussion on the relation between rotation, recirculation length and flow stability, we cannot make the present rotating cylinder flow more stable at high α ⩾ 0.6-1.2 by shortening the base flow recirculation region along the streamwise direction because the value of x s bounces back when α ⩾ 0.6-1.2.We provide a further mechanistic explanation of this result in the structural sensitivity analysis section 3.2.3.Besides, the higher the rotation rate, the less sensitive the value of x s is to the variation of the Reynolds number.
For fixed Re, the effect of rotation rate α on x s is not monotonous, unlike the sphere rotating along the streamwise axis (Sierra-Ausín et al. 2022).Sierra-Ausín et al. (2022) also suggested that x s does not have to increase monotonically as α increases after the bifurcation, which can also be seen in the time-averaged recirculation length of the sphere rotating along the streamwise (Kim & Choi 2002; Lorite-Díez & Jiménez-González 2020).

Drag and lift coefficients
Now, we turn to figure 4 to discuss the drag and lift coefficients in the short rotating cylinder flow.Starting from the fixed finite-length cylinder case (Yang et al. 2022), panel (a) shows that the drag coefficient increases as the rotation rate increases, which is caused by the fact that the rotation makes the negative pressure in the leeward area smaller, but has little effect on the maximum positive pressure at the stagnation point on the windward side.These can be seen from the pressure contour in panels a, b and c, d of figures 2, and more quantitative results on the lift and drag coefficients varying with α are shown in figures 28(a) and (b) in Appendix C. A comparative analysis reveals that the augmentation of α increases the lift coefficient of both the 2-D infinite cylinder and present 3-D finite-length cylinder.However, in the 3-D case, the lift coefficient ceases to increase further when α approaches 2, whereas that in the 2-D case increases monotonically with approaching 2.5 (Stojković et al. 2003;Mittal & Kumar 2003).The effect of α on drag shows contrasting behaviors between the 2-D infinite cylinder and present 3-D finite-length cylinder cases within α < 2.5.Note that our analysis is based on the drag and lift in the SFD base flow, whereas the drag and lift for the 2-D infinite cylinder are time-averaged (Stojković et al. 2003;Mittal & Kumar 2003).Despite the presence of rotation, it is evident from figure 28 (a) that the predominant contributor to the drag is the pressure differential.As the rotational speed increases, the proportion of pressure drag C dp to total drag C d becomes greater.Similarly, the amplitude of the time-averaged lift coefficient |C lp | generated by the pressure differential remains a predominant component of the total lift coefficient |C l |.Furthermore, emphasis is placed on the primary determinant of drag, namely, the pressure differential.As shown in figure 5, the distribution of spanwise-averaged pressure P s (panel a) acting on the upper surface, as well as its component P sx = P s cosθ (panel b) along the x-axis, is characterized with respect to the angle θ.As an example, the areas corresponding to the shaded regions in panel (b) are denoted as Σ 1 , Σ 2 and Σ 3 for the case α = 2. So, the drag resulting from the pressure differential can be regarded as the net area enclosed between the blue solid line and the P sx = 0 axis, specifically Σ 1 − Σ 2 + Σ 3 .Consequently, it can be observed that the increase in drag is primarily attributed to a significant reduction (see panel a) in pressure P s on the back surface of the cylinder due to rotation, resulting in a substantial enlargement (see panel b) of the area Σ 3 .
As with a fixed infinite cylinder (Schlichting & Gersten 2016), a fixed sphere (Johnson & Patel 1999) and a fixed finite-length cylinder (Yang et al. 2022), within the parameter range shown in figure 4

Global stability analysis
In this section, we will discuss the Hopf bifurcation diagram for the ar = 1 case; the effect of aspect ratio will be discussed in section 3.2.3.As the low-and high-rotation-rate flows present different behaviours, we will discuss them separately.
3.2.1.At low rotation speed 0 ⩽ α ⩽ 0.3 We discuss figure 6 and figure 7 collectively, showing the eigenspectra and the eigenfunctions respectively.Figure 6 illustrates the influence of parameters (Re, α) on the leading eigenmodes at low rotation rate α ⩽ 0.3.The global LSA of the flow past the short rotating cylinder shows that there exist two linear unstable modes (LA & LB) for Reynolds numbers Re < 360.The linear growth rates σ of modes LA and LB increase linearly with Re in the vicinity of the instability, as shown in Figure 6(a).When increasing Re for the cases (α = 0, 0.025 and 0.1), both mode LA and mode LB become unstable through a Hopf bifurcation because the frequency at the neutrally stable condition is non-zero, as shown in figure 6(b).Mode LA undergoes the Hopf bifurcation prior to mode LB at a lower Re.At α = 0.15, the critical Re for mode LB is smaller.
From panel b, one can also see that the frequency of global mode LA is not sensitive to α or Re with an approximate value around 0.14, while the frequency of global mode LB increases rapidly with the increase of rotation ratio α.Combining the growth rate σ and the frequency St, we plot in panel 6(c) the modes LA and LB for the representative points LP1-LP5 and the LT point (which can be found in figure 8 for their Re, α values).The implication of panel c is similar and, thus, we will not discuss it further.
Figure 7 plots four eigenmodes LA, LB, LC, LD at α = 0.1 and Re = 290.Only mode LA is unstable; see also the solid green lines in figure 6(a) and one of the unnamed hollow diamonds on the Re-α plane in figure 8(a).The structure of the low-frequency global mode LA is shown in figure 7(a), corresponding to the periodic wake LA (to be discussed in figure 13a obtained by DNS) and is mainly caused by the vortices shedding    7 maintains the same symmetry as the base flow, namely, symmetry with respect to the xy plane, consistent with that of the non-rotating short cylinder (Yang et al. 2022).
Examining closely the solid and dashed lines in figure 6(a), when α is relatively small, mode LA becomes unstable before mode LB when increasing Re.In the range of larger α, mode LB starts to become unstable before mode LA.The intersection point of the σ-Re lines of mode LA and mode LB denotes the condition where σ LA = σ LB , as represented by the black stars in figure 6(a).This competitive phenomenon also exists in spheres rotating along the streamwise direction (Sierra-Ausín et al. 2022), but differs from the single-mode instability observed in spheres rotating along the transverse direction (Citro et al. 2016).Note that our cylinder is rotating along its axis which is in the transverse direction.
The neutral curves in the Re − α plane associated to the two most unstable modes LA and LB are displayed in figure 8.The present asymmetric steady state is linearly stable in the blank region and linearly unstable in the shaded region, as shown in panel (a).With the rotation speed increasing, the threshold Reynolds number for the instability first decreases slightly and then increases.The critical Reynolds number approximately reads Re c = 279 at α = 0.0375.As the rotational speed continues to increase, the threshold Reynolds number becomes less sensitive to the rotation speed.Thus, the result indicates that the rotation of a short cylinder can influence and control the stability properties of the flow.
The present work designates the overlapping area between the unstable regions of modes LA and LB in figure 8(a) as the bi-unstable region, where points LP2, LP3 and LP4 are located.In this study, the intersection of the neutral curves of mode LA and mode LB is denoted as the codimension-two point LT, where two different bifurcations (due to mode LA and mode LB) occur simultaneously.Similar to the rotating sphere (Sierra-Ausín et al. 2022), the present point LT, as the organizing center of the linear system, represents a turning point of competition between different modes and results in the generation of three distinct wake patterns around it in the nonlinear system (to be discussed in figure 13).The dotted line in figure 8(a) passes the codimension-two point LT in the bi-unstable region, on which the growth rate σ LA = σ LB .Above this dotted line, the growth rate of mode LB is greater than that of mode LA; vice versa.
The frequencies St 0 corresponding to the neutral conditions in panel (a) are reported in figure 8(b) as a function of α.As the frequencies are non-zero, the unstable flows will undergo Hopf bifurcation, that is, the unstable mode will oscillate at a certain frequency.Similar to the flow past a sphere rotating along the transverse (Citro et al. 2016) and streamwise directions (Sierra-Ausín et al. 2022), the St 0 of present rotating cylinders increase rapidly as α increases in the regime of low rotation rates α < 0.3.

3.2.2.
At high rotation speed 0.9 < α < 2 We consider relatively low Re in our work.Consequently, in the medium range of the rotation rate α from 0.3 to 0.9 and Re < 360, the wake is steady without vortex shedding.Thus, we will not discuss this range of parameters.
We continue to study high rotation speeds in the range of 0.9 < α < 2 and Re < 360, where the flow may become linearly unstable as shown in figure 9  one can observe that the global LSA based on the asymmetric SFD steady state (see figure 2) indicates three unstable modes in this space of parameters, namely HA, HB and HC, whose critical conditions are depicted as three lines in the figure.Four typical pairs of Re and α in this figure, denoted as HP1 to HP4, are highlighted and their HA, HB and HC eigenvalues are quantified in table 2. The four probed points will be further analyzed in figure 13.From panel a, one can see that higher rotation speeds in the range of 0.9 < α < 2 can significantly lower the critical Reynolds number for instability, which is different from the low-rotation-rate cases 0 ⩽ α ⩽ 0.3.
As shown in panel 9b, the three eigenmodes HA, HB, HC are characterised by different frequencies St 0 , namely, low frequency for HA, immediate frequency for HB and high frequency for HC.The non-zero frequencies of these modes again indicate that the unstable wake flow experiences Hopf bifurcation.By comparison to figure 8(b), the influence of the rotation on the eigenfrequencies of HA, HB and HC is less significant compared with the mode LB at low rotation rate.This indicates that in the high-rotation-speed regime, changing the value of α affects less the frequency in the flow.For the following discussions, when a nonlinear wake flow possesses multiple characteristic frequencies of these global modes at the same time, we will name the wake flow by combining the modes; for example, if both the eigenfrequencies of HA and HB modes are observed in a nonlinear wake, we will call it HAB (see figure 13 to be discussed).
Besides, the three unstable modes can also interact with each other, resulting in the three turning points HT1, HT2, HT3 as shown in figure 9(a).In the parameter ranges of 0.9 < α < 1.64 and 203 < Re < 335, the neutral curves of global modes HA and HB almost coincide, that is, the low-frequency mode HA and mid-frequency mode HB simultaneously become linear unstable.For α > 1.64 (see e.g. the point HT1 in figure 9a), the Hopf bifurcation in the flow begins to be dominated by the high-frequency mode HC.In general, the differences of the Hopf bifurcation in our case and that in the rotating sphere are summarised as follows.
(1) When α > 0.7, Citro et al. (2016) reports only one unstable mode in the flow past a rotating sphere.However, there are three unstable modes in our case, indicating that the Hopf bifurcation process in the present finite cylinder may be more complex.
(2) The eigenfrequency of a sphere is in general higher than that of the rotating cylinder (see the black line with circles in panel 9b).
In figure 10, we show the flow structure of the global eigenfunctions for the three eigenmodes HA, HB, HC at the point HP3 in figure 9.It can be seen that global mode HA (figure 10a) and mode LA (figure 7a) have a similar wake structure, and the difference is that the greater rotation rate makes the transverse offset of mode HA larger.The mediatefrequency mode HB and the high-frequency mode HC have smaller flow structures than those in mode HA.

Structure sensitivity
The global modes detailed above characterise the perturbation growth.The sensitivity of the flow to the perturbation can be studied via the structural sensitivity analysis.To further study the destabilization mechanism, such an analysis for the selected lowrotation and high-rotation flows is conducted.Figure 11 depicts the case of Re = 290, α = 0.1, ar = 1 for a low-rotation case.The transparent red and opaque blue isosurfaces, computed according to ζ in equation 2.12, delimit the wavemaker region, which is the superposition of the leading global mode and its adjoint mode.For both the unstable LA and LB modes, we compute their adjoint modes in the global LSA based on the SFD base flow and time-mean base flow, respectively.In general, the wavemaker region is located in the near wake region of the cylinder (in the top-right region of the xy plane).Its structure remains the same symmetry as the base state, being symmetric with respect to the Oxy plane.The spatial distributions of wakemakers based on mean flow and base flow are similar.Since the wavemaker region indicates the most sensitive region in the flow, one can infer from these observations that (i) the region responsible for the instability is located in the recirculation region behind the cylinder, and (ii) the instability mainly amplifies the perturbations near the cylinder surface (Citro et al. 2016).
Figure 12 presents the flow sensitivity for the high-rotation case of HP3, showing the superposition of the global modes HA, HB, HC based on the SFD base flow with their respective adjoint modes.One can see that the three wavemaker regions are differ sigificantly, in contrast to the low-rotation case.The wavemaker region for the HB mode is located further downstream of the rotating cylinder compared to the other two modes, whose sensitivity regions are close to the cylinder.This implies that the control of the unstable modes in the high-rotation case can be treated separately; especially, control of the HB mode may be achieved more easily as it is not mingled with other modes in the spatial distribution.

Comparison with nonlinear results
In the previous sections, we have identified the global modes in the linearised wake flow past a short rotating cylinder.The relevance of these global modes in the nonlinear simulations of the flow should be established and confirmed.Thus, in this section, we will analyze the results of 3-D nonlinear simulations of the flow past a short rotating cylinder and find the trace of the identified global modes therein.Both high and low rotating speeds are considered.

Wakes behind the short rotating cylinder
Some representative spatial structures and the phase diagrams of the lift-drag coefficients are depicted in figure 13 for the short rotating cylindrical wake flow, obtained by the nonlinear DNS.To analyze the wake structure and compare with the results of the global stability analysis, the frequency in the nonlinear saturated system is also computed, i.e., the power spectral density (PSD) in figure 14.The spectra are obtained by calculating the oscillatory part in the time series of the drag coefficient, e.g.
We calculate the PSD from the drag coefficient, instead of the lift coefficient, because the vortex shedding can take place from the end plates and also the curved surface of the short cylinder, as depicted in figure 1.When the vortices shedding from the end plates (causing spanwise oscillations) are much weaker than those shedded from the curved surface, we noticed that the corresponding oscillation frequencies cannot be observed clearly in the fast Fourier transform (FFT) spectra of the lift coefficient C ly .Compared to the lift coefficient, the drag coefficient seems to be a more robust option for analysing the time-history data in our case.When the vortices shed alternately, the frequency of the drag coefficient is twice that of the lift coefficient.
At low rotation, panels 13(a, b) show that the saturated states at points LP2 and LP4 are a limit cycle.The wakes LA and LB represent the wake dominated by vortices shedding from the cylinder' flat ends and the circular arc surface, respectively.Both types of wake structures LA and LB are also observed in the fixed cylinder flow (see wake patterns P3-2 and P3-1 in figure 10 of Yang et al. 2022, respectively).The difference is that, due to the rotation effect, wake LB undergoes a Hopf bifurcation and becomes a saturated wake state.In the non-rotating cylinder flow, wake LB is only an intermediate transitional state.The results of high rotation speeds are shown in panels (c-f ).Three frequencies of oscillations, including low (HA), medium (HB), and high (HC), are identified in the present nonlinear wakes, see previous discussions on figure 9.The wake HA pertaining to the point HP1 (figure 13c) is characterised by a low-frequency oscillation in spanwise direction and its higher-order harmonics (see figure 14c).Its phase diagram of the lift- drag coefficient is also a limit cycle.The wake HAB (figure 13d) is identified with both a low and a medium frequency oscillation at point HP2.The phase diagram indicates a limit torus with two incommensurate frequencies, see also the PSD result in figure 14(d).
From the perspective of the flow structure, the medium frequency oscillation (HB) is caused by the vortices shedding from the cylindrical arc surface, while the low frequency (HA) is associated to the oscillation of the vortices in the z-direction.Currently, we are not able to identify a monochromatic wake HB with only a medium frequency.From the vortex street structure of the wake HAB at point HP2, it can be seen that the medium-frequency content exists (referring to the gray-coloured structure in the upper part).In the next subsection 3.3.2,we will further use the DMD method to decompose the main components in HAB wake to understand the wake HB.The wake HAC (figure 13e) is identified with a high frequency and a low frequency at point HP3, and its phase diagram of the lift-drag coefficient is also a limit torus (see figure 14e for the PSD result).Finally, in figure 13(f ), a main low-frequency oscillation and its higher harmonics are identified at point HP4, along with a broadband of high-frequency oscillations, leading to a very chaotic signal.HP4 is located at the region where the modes HA and HB are both unstable from the linear analysis.The corresponding phase diagram is also more chaotic compared to the previous cases.

Comparison with DMD modes
According to the global LSA results of a non-rotating 3-D finite-length cylinder flow (Yang et al. 2022), at the Hopf bifurcation point, we can accurately predict the vortex shedding frequency through the linear stability analysis of the time-mean flow.Moreover, the eigenfrequency of the steady base flow (solved using the SFD method if unstable) does not differ too much from the nonlinear vortex shedding frequency.In the present work, we want to establish a qualitative/quantitative relationship between nonlinear and linear systems by comparing the frequency and shape of the dynamic mode decomposition (DMD) modes with the linear global modes.
The DMD (Rowley et al. 2009;Schmid 2010Schmid , 2022) and its extensions have been applied extensively in flow analyses, and have been tested and proven useful to identify the See figure 8(a) and table 1 for the definitions of points LP2 and LP4.Panels (a, c) are on the unit circle and panels (b, d) are on the growth-rate-St plane, where the DMD modes 0, 1, 2 and 3 are marked in red, to be discussed in figure 16.
spatiotemporal patterns of nonlinear flow associated with periodic (Bagheri 2013(Bagheri , 2014) ) and quasiperiodic oscillations (Sierra-Ausín et al. 2022), transitional regimes (Le Clainche & Vega 2017) and turbulent channel flows (Le Clainche et al. 2020).The DMD analysis enables a better understanding of the influence of the linear instability on the onset of vortex shedding of the short rotating cylinder flow.For the data to be processed in the DMD analyses, we typically utilised 150 snapshots over five periods of vortex shedding.We have made sure that the simulations reached a steady-state vortex shedding state before capturing the data.In cases where the wake exhibited multiple cycles, we considered the longest cycle when determining the period of interest.
We first recall the results in figure 8(a) that below the codimension-two point LT for α < 0.1294, the steady-state flow transitions supercritically to a wave LA as Re exceeds Re c (represented by the LA curve in this range of α).To demonstrate the supercriticality, we have calculated the Landau coefficient c 1 and c 3 in Appendix A, see figure 24.Above the codimension-two point, i.e. α > 0.1294, the steady-state flow transitions to a wave LB with high frequencies.That is, the codimension-two point corresponds to a double Hopf bifurcation, which characterizes the interaction between mode LA and mode LB.In the overlap shaded area shown on the Re − α plane in figure 8(a), where the modes LA and LB both are linearly unstable in the rotating cylinder flow, we cannot obtain single-periodic state corresponding to LA or LB separately for low α < 0.3 by using different initial conditions.That is, regardless of the initial conditions, our nonlinear flow in the low-rotation-rate regime always converges to the most unstable mode.This is in contrast to the results of Sierra-Ausín et al. (2022) on a rotating sphere, where they have identified a bi-stable region where a single-mode state can be obtained separately and different initial conditions may lead to different flow modes.In figure 15, the DMD spectra of the flows with low rotation rates corresponding to the points LP2 and LP4 in figure 8 are displayed.The results feature fully saturated modes located on the unit circle.
Figures 16 and 17 show the comparison of the DMD modes with the global modes based on the SFD base flow and time-mean base flow at selected points LP2 and LP4, respectively.Figure 16 shows that DMD mode 0 is the time-averaged flow filed.The DMD mode 1 presents alternated flow structures downstream.According to the imaginary part of the DMD eigenvalues, DMD modes 2, 3 are the second and third harmonics of mode 1, respectively.At point LP2, the frequencies of the leading DMD mode (panel b) and the leading global eigenmode based on the SFD base flow (panel e) are 0.1436 and 0.1411, respectively, with a difference of 1.7%.But there is no DMD mode with a frequency close to ω = 0.1782, which is also a linearly unstable mode (figure 16f ).In the global stability analysis, the (unstable) base flow is solved by the SFD method, whereas our DMD method analyses the saturated flow regime, which hints that the difference in the frequencies in the two methods may be reduced if we use the time-mean flow of the saturated regime in the global stability analysis.This has been carried out and the results are shown in the last row of figure 16 with the subscript MF .Now we can see that the frequency of the The above results pertain to the low-rotation-rate cases.Next, we compare the linear and nonlinear results for the high-rotation-rate flows.We plot the DMD eigenspectra in figure 18 for the selected HP2 and HP3 cases.We discuss the DMD modes labelled in red in the figures.For HP2, the DMD modes are plotted in figure 19  the discussions are similar, we will not go into detail about them.To sum up, through the comparison of the frequencies and the shapes of the linear global modes and DMD modes, we can establish a connection between the linear and nonlinear systems.

Effect of ar
In this last section, we will discuss the effect of ar on the global modes in a short rotating cylinder wake flow, as this information seems to be scarce in the literature on the short (rotating) cylinder flows.
Figure 21 shows the neutral stability curves in panel (a) and the shedding frequency in panel (b) for the case of ar = 0.75.The results of the flow past a sphere (Citro et al. 2016) are also shown for a comparison.Similar to the ar = 1 results in figures 8,9, the low and high rotation flows at ar = 0.75 present dissimilar behaviours.In the case of low rotation speeds, two modes LA, LB undergo Hopf bifurcation due to linear instability Figure 22.Neutral stability curves and the corresponding frequencies for the flow past a short rotating cylinder at ar = 2, compared to a vertically rotating sphere (Citro et al. 2016).
around Re = 330, 340, respectively.The linear unstable region of the mode LB is very small as shown in panel 21(a), almost completely inside the unstable region of mode LA.In the case of high rotation speeds, four global unstable modes are identified up to α = 2, successively experiencing Hopf bifurcation, resulting in three turning points (TP2 to TP4).All these modes are characterised by different frequencies as shown in panel (b).The frequencies in the low-rotation-rate flows are close to the frequencies of the slowly-rotating sphere, whereas the frequencies of the high-rotation-rate modes are smaller than those of the corresponding sphere.The structures of the global modes at low and high rotation rates are shown in figures 25,26 and discussed in Appendix B.
The results of ar = 2 are shown in figure 22.In this case, the calculation seems to be more difficult to converge.We will interpret the with caution.The shapes of the corresponding global modes are shown in figure 27 in Appendix B. When the rotation rate is small in this case, we can identify an unstable mode appearing similarly to the LD mode in ar = 1; this mode will be similarly called LD.Further increasing α, the of the most unstable global mode changes, see the colour transition from orange to purple in panels 22(a, b) and the comparison between figure 27(a) and figure 27(b).In the higher α regime, the most unstable global mode changes abruptly, denoted by green and blue lines in figure 22.We tried to converge as many unstable modes as we could, but some calculations were not converged.Thus, we will not go into details for the high-rotation-rate cases in the ar = 2.Such difficulty in converging the 3-D wake flow is not uncommon, reflecting the complex nature of these flows and highlighting more research efforts to decipher their dynamics.
In the end, we consolidate and compare all the significant results in this work spanning a large parameter space with Re ∈ [100, 500], α ∈ [0, 2] for ar = 0.75, 1, 2 in figure 23.To place our results in a more general context, we also compare our results with other rotating bluff-body flows such as 2-D rotating cylinders (Pralits et al. 2010;Rao et al. 2015), spheres (Citro et al. 2016;Sierra-Ausín et al. 2022) and bullet-like body (Jiménez-González et al. 2014).It can be seen from figure 23 that in the low-rotation-rate regime, larger ar renders the flow more unstable as the critical Re decreases from Re c = 340 for ar = 0.75 to Re c = 50 for ar = ∞ (2-D case, Pralits et al. 2010).In the high-rotationrate end (with the maximum rotation rate being α = 2 in our work), it is difficult to summarise a trend of increasing ar from our finite-length cylinders to the infinitely long cylinder.As discussed above, the computations in this regime are more difficult, calling for more research efforts to elucidate the difference.The 3-D instability in the infinitelylong cylinder (Rao et al. 2015) presents a smoother transition in the range of α ∈ [0, 2], different from our 3-D results where distinct behaviours in the low-(α ≲ 0.6) and highrotation rate (0.6 ≲ α ≲ 2) cases can be identified.This is because α = 2 is not a high rotation rate for infinitely-long cylinder flows; in Fig. 1 of Rao et al. (2015), dissimilar wake behaviours in this flow are separated by α ≈ 2.5.The dynamics of the rotating sphere wake flow along the transverse direction (Citro et al. 2016) is similar to that of the short rotating cylinder with ar slightly larger than 1, whereas the rotating sphere wake along the streamwise direction (Sierra-Ausín et al. 2022) looks more dissimilar than ours.This is likely because our cylinder is also rotating along a transverse axis.In the end, the flow past a spinning bullet-shaped bluff body (Jiménez-González et al. 2014) is also shown for a comparison and its low-rotation-rate behaviour appears similarly to our flow.By studying the effect of ar, we can qualitatively connect our results with those for the sphere, the cylinders of infinite length, etc. and explain the difference between these benchmark flows in a large parameter space.

Conclusions
In this work, a 3-D flow stability problem past a short rotating cylinder has been studied.The motivation for considering this flow configuration is due to its applications in various engineering settings and its relevance to flow control strategy by rotation.New flow modes have been identified in our direct numerical simulations and global stability analyses of this flow.The linear results have also been compared with the nonlinear results to find their traces in real flows.The wavemaker region responsible for the instability generation has been delimited.We have also studied the effect of aspect ratio ar to understand how the 3-D flows change with its geometry parameters.
Firstly, for a cylinder with ar = 1, when the rotation rate α is slower than 0.3, the rotation effect only trivially affects the flow past a short cylinder.The two unstable global modes (see figure 8a) in the low-rotation-rate cases resemble those in the nonrotating flows, corresponding to the vortices shedding from the flat ends of cylinder and the vortices shedding from the curved surface, respectively.The rotating effect swings the recirculation region towards the rotating direction and in general decreases the separation bubble length, defined in our work.Besides, the rotation also casts its effect on the flow instability and bifurcation.For example, the rotation strengthens the symmetry breaking caused by the inherent wake mechanism observed in the regular bifurcation without rotation, leading to a similar asymmetric recirculation region in the non-rotating flows generated by the regular bifurcation.We have also investigated the lift and drag coefficients in the short rotating cylinder flow.Larger rotation rates both increase the absolute value of the drag coefficient and lift coefficient in the transverse direction.An interesting correspondence of the lift coefficients between the short rotating cylinder and the rotating sphere with a doubled rotation rate is observed when Re is relatively large.The global stability analyses reveal that the parameter space α can be divided into low and high regimes when Re is relatively low < 500.When the rotation rate is smaller than approximately 0.3, two unstable global modes exist with non-zero frequencies undergoing Hopf bifurcation, whose interaction and competition giving rise to a codimension-two transition state where the two Hopf bifurcations can occur simultaneously.The critical Re at a certain α slightly decreases and then obviously increase with increasing α.When the rotate rate is large, more unstable modes are observed experiencing Hopf bifurcations.The critical Re in this case decreases with increasing α, highlighting the different effects of Re in the two rotation regimes.The eigenvectors as well as their superposition with the corresponding adjoint modes have also been probed.Especially, we observed that the sensitivity region of the mode HB (with a high rotation rate) is distinguished from other modes close to the cylinder, indicating that its control may be achieved separately.
The comparison of the linear and nonlinear results aims to attach more physical significance to the linear analyses.The traces of the global modes are identified in the nonlinear simulations by comparing their frequencies (i.e., eigenfrequency in the linear analysis and the shedding frequency in the nonlinear DNS).The DMD method is employed to conduct the comparison with the global stability analysis based on the steady flow and the time-mean flow (averaged over several oscillating periods).In general, we can find better correspondence of the time-mean flow results with the DNS results, simply because both analyses were applied to the nonlinear saturated oscillation.In term of the phase diagram, the low-rotation-rate cases characterise limit cycles whereas high-rotation-rate cases present an increasing degree of complexity, encompassing limit cycles, limit torus and chaos, reflecting more complex flow structures and dynamics when the cylinder rotates faster.
Then, the effect of the aspect ratio ar has also been investigated.The aim of this investigation is to compare our flow configuration with other bluff body wake dynamics in a large parameter space.Even though more data points are needed, we can find the trend of how the increasing ar renders the short rotating flow more unstable in the low-rotation-rate cases α < 0.6.When the rotation rate is large but less than α = 2, the critical Re decreases with increasing α.Compared with other bluff body dynamics, we found that the dynamics of the rotating sphere wake flow along the transverse direction is similar to that of the short rotating cylinder with ar slightly larger than 1, and the rotating sphere wake along the streamwise direction differs more significantly than our results due to the different flow configuration.
The current work focuses on the first instability in the short rotating cylinder wake flow.The flow dynamics already presents a high degree of complexity.In order to further understand the underlying mechanism, as a future direction, the subsequent flow bifurcations in this flow can be studied in detail by employing the global linear stability and weakly nonlinear stability analyses.

Appendix C. Validation of numerical codes
This appendix demonstrates the validation of the numerical codes used in the present work.The results of α = 0 connects smoothly with the non-zero-α results in the main text as qualitative proof of the accuracy of our codes.To further test the accuracy in a quantitative manner, we compare our simulated results of the lift and drag coefficients in the 2-D rotating cylinder flow with those in Kang et al. (1999); Stojković et al. (2002).
Figure 28 displays a good comparison of our results with theirs at Re = 60 and 100, indicating the accuracy of the used nonlinear numerical code in the current work.For the verification of the linear code, figure 29 presents the wavemaker region in the classical 2-D cylindrical wake flow in Giannetti & Luchini (2007); Marquet et al. (2008).We can see a very good comparison is achieved between our results (top) and theirs (bottom).This good comparison entails the linear code solving correctly both the global modes and the adjoint modes.
In the end, we furnish a test study on the size of the computational domain.In our previous work on the 3-D non-rotating short cylinder (Yang et al. 2021), we converged a suitable computational domain balancing the computational efficiency and accuracy.With the rotation effect, we found that the size of the computational domain should increase to accommodate the swinging effect brought by the rotation.The result is shown in table 3 for different values of L a and L o , see their definitions in figure 1.The symmetry boundary condition is imposed on the surfaces S xy,t , S xy,b , S xz,f and S xz,b , where the position is set to a range of ±10D to ±15D, that is L a = 10 − 15.
As shown in table 3, unlike the non-rotating cylinder (Cadieux et al. 2017), the drag coefficient is more sensitive to the domain size than St number.The relative errors of the drag, lift coefficients and the frequencies from mesh M2 to M4 are less than 1% compared with the reference results of mesh M5.Thus, in order to balance the efficiency and accuracy, the mesh M2 and the order N ord.= 7 are adopted in the present work to compute all computational instances.

Figure 1 .
Figure 1.The computational domain and boundary conditions (not to scale) (panel a) and mesh design (panel b).The red unit vectors (en, eτ , e b ) in panel (a) represent the directional vectors of the surface Sxy,t.The finite-length cylinder is rotating around its axis that is perpendicular to the incoming flow.
.3c) U • e n = 0, (∇U • e τ ) • e n = 0, (∇U • e b ) • e n = 0 on S xy , S xz , (2.3d) where êz = (0, 0, 1) is a unit vector aligned with the positive z-axis.The vector r = (x c , y c , z c ) is the position vector of a point located on the cylinder surface, i.e. the vector from the origin of the coordinate system to the point.Here e n , e τ , e b are the unit normal, unit tangent and unit bitangent vectors, respectively.As shown in figure 1(a), the vector e n of surfaces S xy,t , S xy,b , S xz,f , S xz,b and S out points out the computational domain.The directions of the vector e τ of surfaces S xy,t , S xy,b , S xz,f , S xz,b point along positive Base flows and pressure by SFD method

Figure 3 .
Figure 3. (a) A separation bubble under the effect of rotation in plane oxy, whose area is surrounded by the green curve (Ux(x, y, z = 0) = 0) and the wall of the rotating cylinder.The white dot shows the position of the hyperbolic stagnation point.The colour represents the pressure field.(b) Separation bubble length xs of the SFD base flow as a function of Re, compared with a fixed sphere (Johnson & Patel 1999) and a sphere (Sierra-Ausín et al. 2022) rotating along the streamwise (blue dashed lines).The blue and black dash-dotted lines represent the xs obtained from the neutrally stable SFD base flow of rotating sphere (Sierra-Ausín et al. 2022) and the present cylinder, respectively.The black and red thick lines represent the xs obtained from the SFD base flow and the mean flow of a fixed finite cylinder (α = 0) (Yang et al. 2022), respectively.
Figures 2(e) (Re = 160) and 2(f ) (Re = 330) show the base states at moderate rotation α = 0.6, and both their nonlinear saturation states are steady.So the base flow and mean flow are the same for these cases.Compared with the previous low rotation case, the flow topology further changes, i.e. the recirculation zone in the wake almost disappears, in either the xy or xz plane.But a stagnation point (white point in the panel) similar to that observed in the 2-D rotating cylinder (see figure 3b in Sierra et al. 2020) appears, which is located in the second quadrant of the cylinder.

Figure 4 .Figure 5 .
Figure 4. Drag (panel a) and lift (panel b) coefficients of the steady SFD base flow as a function of Re at ar = 1.Comparison with a fixed sphere (Johnson & Patel 1999), a fixed finite cylinder (Yang et al. 2022), and a sphere rotating about the transverse direction (Citro et al. 2016).(a) (b) (a), an increase in Re generally leads to a decrease in drag coefficient.On the other hand, unlike the drag coefficient C d , the values of the lift coefficient C ly in figure 4(b) present more variation when Re changes.At low rate α ⩽ 0.3, C ly first decreases with the increase of Re and then increases; at high rate α ⩾ 0.6, C ly increases monotonically as Re increases from 50 to 330.Besides, in the interval 0.2 < α < 0.6 of figure 4(b), C ly is relatively insensitive to Re, and the wake is steady (see the neutral stability curves in the next section).It is worth noting that the lift coefficient C ly | α=c of the rotating cylinder and the C ly | α=2c of the rotating sphere gradually coincide for Re ≳ 270 (with c ⩽ 0.1 in present works), see the top-right corner in panel (b).For example, the values of C ly | α=0.025 of a rotating cylinder and C ly | α=0.05 of a rotating sphere, or C ly | α=0.05 of a rotating cylinder and C ly | α=0.1 of a rotating sphere, are approximately the same for Re ≳ 270.

Figure 6 .
Figure 6.The growth rate σ (panel a) and frequency (panel b) of the leading global modes (of the SFD base flow) as a function of Re at ar = 1.The position of the black star symbol '✩' in panel (a) indicates that σLA = σLB.(c) The eigenspectra of global modes LA (gray shaded area) and LB at points LP1, LP2, LT, LP3, LP4 and LP5 in figure 8 (a).

Figure 7 .
Figure 7. Four representative eigenmodes with weak asymmetry in the flow past a short rotating cylinder at α = 0.1 and Re = 290.Mode LA (in panel a) and mode LD (in panel d) are similar to the mode A and mode B reported by Yang et al. (2022), respectively.The new mode LB (in panel b) with high frequency is more asymmetric caused by rotation.The mode LC (in panel c) has zero frequency.The Q-criterion isosurfaces Q = 0 are colored by the x-component of the vorticity ranging from −2 × 10 −3 to 2 × 10 −3 .

Figure 8 .Figure 9 .
Figure 8.(a) Neutral stability curves undergoing the Hopf bifurcation for the flow past a finite rotating cylinder at ar = 1 and the rotating sphere (the mode I by Citro et al. 2016).The specifications of the points LP1-5 and the codimension-two point LT are shown in table 1.The hollow symbols represent the corresponding nonlinearly saturated wake looks similarly to the global mode LA, whereas the solid symbols denote the nonlinearly saturated wake resembles more the global mode LB.(b) Plot of critical Strouhal numbers St0 on the neutral stability curves against rotating ratio α.(a) (a).From this figure,

Figure 10 .
Figure 10.The global modes with strong asymmetry for the flow past a rotating finite cylinder at high rotation rate at Re = 170, α = 1.8 (point HP3).Where panels (a), (b), and (c) represent the modes HA, HB, and HC at point HP3.The Q-criterion isosurfaces Q = 0 are colored by the x-component of the vorticity ranging from −2 × 10 −3 to 2 × 10 −3 .The cylinder centroid is located at (0, 0, 0).

Figure 11 .Figure 12 .
Figure 11.Flow sensitivity in a low-rotation-rate case.The wavemaker isosurfaces are plotted for the first two unstable modes LA (a, c) and LB (b, d) at Re = 290, α = 0.1, ar = 1.Transparent red is for ζ = 0.2 and opaque blue is for ζ = 0.4.Panels (a, b) show the results of the global LSA based on the SFD base flow and panels (c, d) on the time-mean flow.(a) (b) (c)

Figure 13 .
Figure 13.Left panel: the nonlinear wakes spatial structure obtained by DNS.The Q = 0 isosurfaces are colored by the streamwise vorticity ranging from -0.1 to 0.1.Right panels: the corresponding phase diagrams of d − C l for rotating cylinder.(a) wake LA at point LP2, (b) wake LB at point LP4, (c) wake HA at point HP1, (d) wake HAB at point HP2, (e) wake HC at point HP3 and (f ) wake HAC at point HP4.The corresponding Re, α and eigenvalues for each case are shown in tables 1 and 2.

Figure 14 .
Figure 14.The power spectral density (PSD) for the cases (a) LP2, (b) LP4, (c) HP1, (d) HP2, (e) HP3 and (f ) HP4, The PSD is calculated based on the oscillatory part of the time series of the drag coefficient.

Figure 15 .
Figure 15.The DMD spectra of point LP2 in panels (a, c) and of point LP4 in panels (c, d).See figure8(a) and table 1 for the definitions of points LP2 and LP4.Panels (a, c) are on the unit circle and panels (b, d) are on the growth-rate-St plane, where the DMD modes 0, 1, 2 and 3 are marked in red, to be discussed in figure16.

Figure 17 .
Figure 16.Comparison of the DMD modes with the global modes at the point LP2 (α = 0.05, Re = 305).Four DMD modes and the corresponding DMD eigenvalues are shown in the panels (a − d).The global modes LA (e, g) and LB (f, h) in the global stability analysis based on the SFD base flow (e, f ) and mean flow(g, h).All figures colored by streamwise vorticity.The eigenfrequencies of the DMD modes are in good agreement with those obtained from FFT method (figure14 a), with a relative error of 0.48%.Moreover, the eigenfrequencies of the DMD modes match well with the eigenfrequencies of the mean flow global modes (panel g and f ), with a relative error of 0.83%.However, the difference between the characteristic frequencies of the DMD modes and the eigenfrequencies of the SFD base flow global mode is slightly larger, with a relative error of 2.22%.Additionally, the topological structure of the DMD modes is almost identical to that of the mean flow global modes.

Figure 18 .Figure 19 .
Figure 18.The DMD eigenvalues spectrum of points HP2 (a, b) and HP3 (c, d) on the unit circle and on the growth-rate-St plane.
(a, b) along with the global modes based on the SFD base flow (c, d) and the time-mean flow (e, f ).Compared with the case of a low rotation rate in figures 16,17, the similarity of the DMD modes with the global modes in the high-rotation-rate case is greater; for example, the first column in figure 19 shows that the three modes look similar, and the global mode based on the mean flow is again slightly better compared to the DMD mode.In the second column, we can also find DMD mode 2 that resembles the global modes in the global stability analyses, where the global mode based on the time-mean base flow looks closer to the DMD mode.The same conclusion can be drawn for the HP3 point in figure 20.As

Figure 20 .Figure 21 .
Figure 20.Comparison of the DMD modes (a, b) with the global SFD base flow modes (c, d) and mean flow modes (e, f ) at point HP3.The corresponding eigenvalues are shown in the panels (a − f ).All figures colored by streamwise vorticity.The eigenfrequencies of the DMD modes are in good agreement with those obtained from FFT analysis (figure14 e), with a relative error of 0.7%.However, the frequencies of DMD modes are not in complete agreement with those of mean flow and SFD base flow, with errors of 5.8% and 6.6%, respectively.Additionally, the topological structure of DMD modes is nearly identical to that of mean flow global modes, as compared to the SFD base flow.

Figures 25 ,
Figures 25,26 display the structures of global linear modes for ar = 0.75 at low and high rotation rates, respectively.In figure25, the modes LA and LD are presented.They are named as such because they resemble the LA and LD modes in the ar = 1 flow as shown in figure7(a, d).Figure26features the global modes HA1, HA2, HB and HC.Again, comparison can be made to the global modes for the ar = 1 flow in figure10at a similar rotation rate.The structures of some selected global linear modes for ar = 2 are displayed in figure27.Panel a shows an unstable mode at a low rotation rate α = 0.1.It is called LD Figures 25,26 display the structures of global linear modes for ar = 0.75 at low and high rotation rates, respectively.In figure25, the modes LA and LD are presented.They are named as such because they resemble the LA and LD modes in the ar = 1 flow as shown in figure7(a, d).Figure26features the global modes HA1, HA2, HB and HC.Again, comparison can be made to the global modes for the ar = 1 flow in figure10at a similar rotation rate.The structures of some selected global linear modes for ar = 2 are displayed in figure27.Panel a shows an unstable mode at a low rotation rate α = 0.1.It is called LD

Figure 28 .Figure 29 .
Figure 28.Validations by comparing the lift and drag coefficients of 2-D rotating cylinder flow between the present DNS code results with those in Kang et al. (1999); Stojković et al. (2002).(a) time-averaged drag coefficient; (b) time-averaged lift coefficient; (c) amplitude of C l ; (d) amplitude of C d ; (e) Strouhal number calculated using the C l signal.Besides, the drag and lift coefficients of steady SFD base flow are added in panels (a) and (b) for comparative analysis, and decompose them into Cl = Clp + Clv , Cd = Cdp + Cdv for the case (ar = 1, Re = 160).(a)

Table 1 .
Yang et al. (2022) of the five typical points in the parameter space (Re, α) (see also figure8a) for the cylinder ar = 1.Subscripts LA and LB represent global modes LA and LB (figure 7 a and b), respectively.The eigenfrequencies marked in bold are close to the frequency StDNS in DNS (e.g., see also figure14) of the saturated nonlinear wake.someothercases,such as LT, may take this form.LT, as shown in figure8(a), represents a low-rotation-rate transition state where the LA and LB modes are both neutral, to be discussed further below.Mode LD in panel d appears similar to the mode B identified inYang et al. (2022)for the non-rotating short cylinder.Each global mode in figure

Table 3 .
Mesh La × Lo Ntot.(N ord. ) C d A grid sensitivity test for the nonlinear DNS case Re = 290, α = 1.2, ar = 1, ∆t = 10 −3 .As shown in figure1, La is length from surfaces Sin, Sxz and Sxy to cylinder centre; Lo is length from surface Sout to cylinder centre.Ntot. is the total number of hexahedral elements inside the computational domain.N ord. is the polynomial order of each hexahedral element.