Chaotic vortex induced rotation of an elliptical cylinder

Non-linear oscillations of an elliptical cylinder, that can rotate about an axis that passes through its symmetry axle due to a torsional spring and hydrodynamic torque produced by the flow of a Newtonian fluid, were analysed in terms of a single parameter that compares vortex shedding frequency with the torsional spring's natural frequency. The governing equations for the flow coupled with a rigid body with one degree of freedom, were solved numerically using the Lattice Boltzmann Method (LBM). The Reynolds number used was $Re=200$ which, in absence of torsional spring, produces chaotic oscillations of the elliptical cylinder. When the torsional spring is included, we identified three branches separated by transition regions when stiffness of the restorative torque change, as in the case of vortex induced vibrations (VIV's). However in this case, several regions presenting chaotic dynamics were identified. Two regions with stable limit cycles were found when both torques synchronized and when stiffness of the torsional spring is big enough so that the ellipse's oscillation is small.


I. INTRODUCTION
The study of fluid-solid interactions has been a problem of great interest for a long time from varying points of views.Vortex induced vibrations on solid structures as elastically mounted cylinders [1][2][3] , ellipses [4][5][6] , flexible cylinders and cables [7][8][9] are subjects of research in many fields of engineering due to their wide range of practical applications.In addition, from a mathematical point of view, the study of non-linear phenomena due to flow of a viscous fluid around rigid bodies when different numbers of degrees of freedom are considered [10][11][12] and chaotic dynamics of the rigid body due to the flow [13][14][15] , has attracted the attention of several researchers because of modern development in computers and numerical methods.
Movement of an elastically mounted cylinder is modeled using linear and non-linear forced oscillators, depending on the value of the so called reduced velocity 1,2 ; a parameter that compares vortex shedding frequency and the cylinder's oscillation frequency without flow.In such case when only translational degrees of freedom are considered, the cylinder's movement is modeled as a forced-damped oscillator in the lock-in region, where amplitudes of oscillation transversal to the flow direction are the largest 1,2 .
When rotational degrees of freedom are considered, movement of the rigid body is due to an auto-rotation effect that was originally studied by J. C. Maxwell in the 1890's 16 and was later characterized by H. J. Lugt 17,18 .Vortex shedding behind the body produces a torque that oscillates with time and, depending on the Reynolds number, behaves as a nonlinear function of angular position, angular velocity and time.Non-linear dynamics observed in 2D and 3D have been analyzed from different points of view; for instance in the design of power extraction systems 4 and the study of non-linear oscillators 13,14 .
In this work we present numerical solutions of the flow of a viscous fluid around an elliptical cylinder that can rotate due to external restorative torque, with stiffness k, and hydrodynamic torque due to the flow.Reynolds number was fixed at Re = 200 where oscillations of an elliptical cylinder without restorative torque are known to be non-linear 4,13 .
A detailed analysis of cylinder and flow dynamics was made as function of a single parameter, proportional to the natural frequency of restorative torque.The problem shows a similar behavior to that observed in cylinders with translational degrees of freedom: an excitation region where oscillation is periodic and a lock-in region where amplitude of oscillation reaches maximum values.A third region can be identified with a decoherence region, where amplitude of oscillation is diminished.These regions are separated by transition regions and both show different non-linear behavior and chaos.The transition between decoherence and lock-in regions is a chaotic region with an associated strange attractor.In the transition between excitation and lock-in regions, forced-damped periodic oscillations were observed when k → ∞.Then, as k is diminished, the cylinder's oscillations evolve into a limit cycle through a region where the ellipse's amplitude in each oscillation changes erratically.Vorticity generated in the flow depends on the elliptical cylinder's degree of confinement which is controlled with the spring's stiffness.This article is organized with the problem statement in section II including a brief description of the numerical method.In section III the most relevant results are presented and, in section IV, we present a discussion of our main findings.The problem consists on uniform flow of a Newtonian viscous fluid, with mass density ρ and kinematic viscosity ν, around an elliptical cylinder that can rotate around its symmetry axis due to hydrodynamic torque and is attached to a torsional spring with constant stiffness k.The elliptical cylinder has a mass density ρ s , a and b are the semi-minor and semimajor axes respectively and θ is angle of rotation, as shown in Figure 1.We did not include structural damping to enforce the ellipse's oscillations in this work.The torsional spring's equilibrium angle is fixed at θ e = 0.
In order to work with non dimensional variables, space is scaled with b, velocities with velocity far from the body U, and times with b/U.With this choice, the problem is characterized by four non dimensional parameters: Reynolds number Re = Ub ν , modified Strouhal number S t = k I b U that compares natural frequency of restorative torque with a scale of the vortex shedding frequency, ratio between semi-major and semi-minor axis α = b a , ratio between solid and fluid mass densities m * = ρ s ρ f .Considering u(r,t) and P(r,t) as non dimensional velocity and pressure fields (pressure is scaled with characteristic viscous pressure ρ f νU b ), the governing equations for this problem are given by which correspond to non dimensional Navier-Stokes equations for an incompresible flow.Boundary conditions in this case are where î is the unitary vector in the x direction, r s is a point in the ellipse's surface and v s (r s ,t) its velocity.The above conditions correspond to non-slip at the ellipse's surface and a uniform flow far from the body.Rotation of the ellipse is solved using the Newton equation for torques given by where I * = πm * α(1 + α 2 )/4 is a reduced moment of inertia and dimensionless hydrodynamic torque is where R and ω(t) are the position of the ellipse's center of mass and angular velocity of rotation, respectively; σ represents the stress tensor of a Newtonian fluid.

A. Numerical Method
To find numerical solutions for the above problem we used a two dimensional, nine neighbors (D2Q9) lattice-Boltzmann model 19 .The proposed algorithm, that includes moving immersed boundaries, has been validated in previous works 11,20 .
In this method space is discretized using a squared lattice.Lattice spacing as well as time steps can be conveniently set to unity.The fluid's state at the node with vector position r at time t, is described by the particle distribution function f k (r,t) that evolves in time and space according to where τ is a relaxation time related to the fluid's kinematic viscosity ν = (τ − 1/2)/3.Equilibrium distribution function f (eq) k is given by which corresponds to a discrete Maxwell distribution function for thermal equilibrium.In the above expressions, macroscopic density and velocity fields are computed using and microscopic velocities e k are given by where w 0 = 4/9, w k = 1/9 for k = 1, . . ., 4 and w k = 1/36 for k = 5, . . ., 8. Notice that, with this set of microscopic velocities, expression (7) is always evaluated at lattice points.It is well known that the above procedure approximates solutions to the Navier-Stokes equations in the limit of small Mach numbers 19 .Equation ( 7) provides an explicit algorithm for updating all distribution functions f k at a given node in the lattice, as long as its 8 nearest neighbouring nodes are inside the fluid's domain.For nodes adjacent to a solid wall, distribution functions coming from neighbouring nodes outside the fluid domain must be provided.We chose the set of boundary conditions proposed by Guo and Zheng in Ref. 21 for curved rigid walls, which approximately gives the no-slip boundary condition (3).The force and torque acting on the body were computed using the momentum-exchange method of Mei et.al. in Ref. 22 and then equation ( 5) is solved using a leap-frog method.
In order to simulate an infinite domain, a constant velocity was enforced at the inlet whilst in the rest of the boundaries a stress-free condition was implemented (see for details F. Mandujano and C. Málaga 23 ).The domain's size was chosen to be big enough such that effects due to lateral walls were minimized but the internal boundary could have enough number of points.The numerical scheme was implemented to run in parallel in Graphical Processor Units (GPU) due to the large number of nodes involved in these simulations.

III. RESULTS
Simulations started with a uniform flow and the ellipse fixed in space, until the fluid reached a time dependent flow with vortex shedding.Then the solid, initially at θ (0) = θ 0 , was allowed to rotate following equation (5).Aspect ratio was fixed at α = 2 with m * = 1 and, as mentioned before, the Reynolds number was fixed at Re = 200.Numerical experiments were performed varying the modified Strouhal number S t .In absence of restorative torque (S t = 0), the hydrodynamic torque oscillates with time in a non-periodic fashion with a non-zero mean.Hence, for any choice of initial condition θ 0 (initial angular velocity is set to zero in all cases), the ellipse's semi-mayor axis will oscillate around the vertical axis, as shown in Fig. 2-(a).When θ 0 is chosen in the neighborhood of θ 0 = 0, trajectories in phase space are expelled from the origin and, as time grows, are attracted making open orbits around either θ c = π/2 or θ c = −π/2 (see fig 2 -(a)).The elliptical cylinder is then in an unstable equilibrium position and starts to move due to auto-rotation 17,24 .We identified θ s = 0 as a fixed-point that behaves as a saddle point (see Fig. 2-(b)).
Similarly, when θ 0 = ±π/2 much like spiral points trajectories in phase space are expelled outwards and also make open orbits in phase space, as shown in In any case, observed orbits oscillate inside a region around the fixed-point.This region thins out as the Reynolds number decreases and converges into a stable limit cycle.Stability properties of these three fixed-points depend on the Reynolds number, θ s = 0 is a saddle point as long as there is vortex shedding behind the body.The other two fixed-points, θ c = ±π/2, have an associated stable limit cycle at lower Reynolds number values.
According to Williamson and Roshko's nomenclature 25 , the wake produced in this regime is a 2S pattern.The vortices shed in each cycle seem to organize in a von Kármán wake close to the body, downstream vortex pairs interact with each other changing their pattern, as shown in figure 2-(c) where the elliptical cylinder at two points in time separated by one cycle is shown.Even though the two points in time seem to be similar, in each cycle the vorticity produced on the cylinder's surface changes slightly modifying the pattern downstream.Similar results were found by Weymouth 13 , who studied a rotating ellipse towed with a constant velocity, when the axle of rotation passed through its geometrical center.In our case, the ellipse's response is a non periodic oscillation with an amplitude that varies in a complicated manner; showing that hydrodynamic torque is a non-linear function that plays the role of forcing and damping forces in the language of nonlinear oscillators 26 .In VIV's literature, hydrodynamic forces due to vortex shedding processes are modeled as oscillatory periodic functions of time 1,2 .
When torsional spring is included with its equilibrium position at θ e = 0, there is a competition between restorative and hydrodynamic torques.When 0 < S t < 0.8, we found a similar behavior to the one described above (see figure 3-(a)).All initial conditions drive the elliptical cylinder to oscillate around one of two fixed points ±θ c that have moved towards the point θ s which is still a saddle point, as shown in figure 3-(b).As in the previous region, orbits are confined in regions that thin out as the Reynolds number decreases.Mean amplitude remains approximately constant within this region, while the power spectral analysis of θ (t) shows several frequencies.The highest peak in power spectral density for each experiment (main frequency) is within the interval [U/b, 2U/b] and seems to change rather erratically due to the oscillation's non periodicity.In contrast with the case when S t = 0, symmetry properties of trajectories in phase space change.In the former case, phase space trajectories are symmetric with respect to the origin.When S t = 0, trajectories are skew-symmetric under a reflection in θ (t).The elliptical cylinder rotates faster when it moves away than when it moves towards the saddle point θ s .The vorticity field also shows different symmetrical properties, as shown in figure 3-(c), where the ellipse at approximately the same angular position after one cycle is shown.
The vorticity generated at the cylinder's surface is rather different from one cycle to the next.Form, behavior and the way in which vortices interact with one another also change.They are deflected in the cross-flow direction which results in a wider wake.When 0.5 ≤ S t < 1.1 the ellipse starts to oscillate around either θ c or −θ c depending on initial conditions as in the previous region.After a time interval, amplitude of oscillation increases resulting in the elliptical cylinder being attracted to the opposite fixed-point and starting to oscillate around it for another time interval, and so on (see figure 4-(a)-(c)).Residence times, where the ellipse oscillates around either of the fixed-points ±θ c , seem to be distributed randomly.The two fixed-points ±θ c approach to θ c = 0, as S t gets close to one, where they both merge with the saddle point to form what seems to be a strange attractor.Dynamics are very similar to the ones found with the Lorenz system and Duffing oscillator 26 .In figure 4-(c) the vorticity field is shown at two points in time after one oscillation.In this case, it is found that shed vorticity differs and vortex pairs are further deflected vertically downstream compared to the previous region.
For 1.1 < S t < 1.2 (shedding and spring frequencies are approximately the same) we found stable limit cycles, as shown in figure 5-(a) and (b).Observed θ (t) frequencies are approximately U/2b, half the shedding frequency, and amplitude of oscillation reaches its maximum values.Both lift and torque oscillate within the the same frequency as θ (t) and drag coefficient with shedding frequency U/b.Now, the only fixedpoint seems to be θ s = 0 with a skew-symmetric stable limit cycle.Vorticity field is a 2P+2S wake, as shown in figure 5-(c) where we show two points in time separated by half a cycle of the ellipse's rotation.
Notice that for each half a cycle, two structures are shed: The wake seems to be an internal inverted von Kármán street with an external one as shown in figure 5-(c).As S t is increased outside the previous region, orbits oscillate around a thicker region and the rotational angle loses its periodicity with amplitudes that vary from one cycle to the next by small amounts, as can be appreciated in figure 6-(a).In figure 6-(b) three types of trajectories are shown with different initial conditions.When the initial condition is ±θ 0 ( = 0) we found two skew-symmetric solutions and, when the initial condition is close to θ 0 = 0, solutions found seem to be a combination of the first two.Hence, even when it seems that fixed-points collapse at the saddle point position θ s , influence of the two fixed-points ±θ c can still be seen, which shows that basins of attraction are divided into at least three regions.
Vorticity fields are made by 2S patterns, as in von Kármán wake, however, when angular velocity has changes as θ (t) approaches the saddle point (see figure 6-(b)) vortices are shed with slightly different amplitudes that modify the pattern downstream, as shown in figure 6-(c).In contrast with the initial region, where the amplitudes of oscillation around either of the fixed-points ±θ c remain approximately constant, in this region amplitude of oscillation rapidly decreases and, at some value of S t , reaches a minimum value and a limit cycle is found.The cylinder's oscillations in this region are both very small and periodic with a frequency of two times the shedding frequency U/b.
Figure 7 shows the root mean square of the angle of rotation as a function of S t .As mentioned before, for small val-FIG.7. a) Root mean square of θ (t) as function of S t .b) Evolution of ±θ c as function of S t ues of S t , cylinder trajectories in phase space seem chaotic, variations in mean amplitudes of oscillation are small and remain approximately constant as S t grows.At some value of S t , there is a transition where bigger oscillations are found, which correspond to time intervals when the cylinder passes from oscillating around one fixed-point to the other (shaded region in figure 7-(a)).Liapunov coefficients in these two regions are positive and have approximately the same value of λ ∼ 0.1.The differences in initial conditions used to compute λ were 10 −7 .
Evolution of fixed-points is shown in figure 7-(b).The starting point is θ c = ±π/2 for S t = 0 and then θ c decreases monotonically as S t grows.In the region of transition where the cylinder changes from chaotic oscillations around either of the fixed-points ±θ c to chaotic oscillations changing form one to the other intermittently, the behavior seems to follow the same trend.When S t ∼ 1 fixed-points jump and merge with the saddle point at θ s = 0, marking the transition's boundary.
The cylinder's amplitude of oscillation reaches its maximum when 1 < S t < 1.2, a tiny region where we found a stable limit cycle.The only fixed-point is at θ s = 0, as shown in figure 7, the other two points collapse with θ s .However, their influence is still seen in the skew-symmetric phase space trajectories.As S t is increased, A θ starts to decrease monotonically and phase space trajectories change their geometry, showing the influence of two fixed-points associated with hydrodynamic torque.When S t ∼ 1.6, amplitude A θ decrease faster and becomes small but different from zero and remains constant.Trajectories in phase space are ellipses, oscillation is harmonic and, when S t → ∞, the only fixed-point becomes an attractor.

IV. DISCUSSION
Non-linear oscillations of an elliptical cylinder that rotates around its axle due to flow of a Newtonian viscous fluid and a restorative force are discussed.When the only torque on the elliptical cylinder is due to fluid flow, the ellipse's major-axis oscillates around the axis in the cross-flow direction.Without fluid flow, the ellipse rotates harmonically around θ = 0 with a frequency of k/I.Reynolds number used was Re = 200 where it is known that hydrodynamic torque produces nonlinear oscillations.We used S t as a central parameter, which results from a comparison between the characteristic vortex shedding U/b and spring k/I frequencies.Results found showed similar behavior to vortex induced vibrations (VIV) around a circular cylinder with two and three degrees of freedom 1,8,11,27 .Behavior of oscillation amplitudes in the cross-flow direction is classified into three regions as a function of reduced velocity, which is proportional to the inverse of S t .The excitation region is where the cross-flow amplitude of oscillation is small but starts to increase.The lock-in and desynchronization regions are regions where the cross-flow amplitude reaches its maximum and begins to decrease, respectively.
In our case, the excitation branch can be identified as the region when S t > 1.8, as shown in figure 8, where a bifurcation diagram for this problem is portrayed.The fixed-point θ s = 0 behaves as an attractor when S t → ∞ (not shown on bifurcation diagram).Hydrodynamic torque works as a damping term in equation ( 5) since every initial condition leads the ellipse to a rest, as in a damped pendulum.As S t decreases, the ellipse starts to oscillate periodically with a period of 2b/U and an amplitude that remains approximately constant, similar to a forced damped pendulum.
The excitation branch is followed by an abrupt growing of amplitudes of oscillation which lead to the lock-in branch, where the cross-flow amplitude reaches its maximum.In our case, as S t decreases, amplitude starts to grow abruptly and reaches a maximum where we found a stable limit cycle.Around the maxima, we found a region where oscillation is periodic, which we identified as the lock-in region.Transition between excitation and lock-in regions occurs in a region where the stable limit cycle evolves in a new structure influenced by two fixed-points at ±θ c .Amplitude reached at each oscillation varies, producing asymmetries in the branches of the parabola-like curve in figure 8.As S t decreases, observed variations in amplitude diminishes and reaches zero.As S t approaches one, amplitudes collapse to a single value with a more uniform growth than in the transition region.
After the lock-in branch, the stable limit cycle evolves in what seems to be a strange attractor (shaded region in figure 8).Even when it seems that amplitudes of θ (t) still explore a wide range of values, amplitude of oscillation decreased since the ellipse oscillates intermittently about either θ c or −θ c .After the transition region, amplitudes of oscillation remain approximately constant and fixed-points ±θ c increase towards ±π/2 as S t approaches zero.We identified this region as the desynchronization zone.However, in our case, the ellipse's oscillation is non-periodic and its amplitude is larger than in the excitation branch.
Within the transition region the ellipse oscillates intermittently around fixed-points ±θ c .In the language of VIV's on circular cylinders, this region of transition is identified with a hysteresis process, meaning that it is different to go from the excitation region to the lock-in region than the other way around 1,8,28 .There is also evidence in 2D numerical simulations that transition between these two regions can be made following different paths by changing the set of initial conditions, because there are several coexisting solutions between the excitation region and the lock-in branch 11 .
The system under consideration shows a close behavior with the problem of VIV's with 2 degrees of freedom in a sense that similar regions or branches can be identified.However, in this case, we only found periodic solutions in a small region (identified as the lock-in branch) and when S t → ∞ which corresponds to a fixed ellipse.When S t = 0 the ellipse's oscillation is chaotic due to the non-linear nature of hydrodynamic torque.As S t is increased, degree of confinement in cylinder movement is modified and shed vorticity changes, which is due to a competition between hydrodynamic and restorative torques.Therefore, only when both torques synchronize or the movement is confined to very small oscillations, obtained solutions are periodic.In the rest of analyzed regions, the hydrodynamic torque dominates the flow and the ellipse's oscillation shows chaotic dynamics.For a future works, we will be extending the study to ellipses with three degrees of freedom and will include other bodies, to study solid-fluid-solid hydrodynamic interactions.

FIG. 1 .
FIG. 1. Diagram of the problem.The flow around an elliptical cylinder that can rotate around its symmetry axis due to hydrodynamic and restorative torques.The flow is considered to be uniform, with velocity U far form the body, a and b are the minor and major axis of the ellipse respectively and θ = θ (t) is the angle of rotation.

FIG. 2
FIG. 2. a) Rotational angle as function of time for different initial conditions with S t = 0. b) Phase space trajectories correspond to those shown in (a), arrows indicate the trajectories' starting point.c) Vorticity field at two different times separated approximately by one cycle.
Fig 2-(a) and (b).Time signals for different initial conditions are qualitatively similar but are different showing sensitivity in initial conditions.

FIG. 3
FIG. 3. a) Rotational angle as function of time for different initial conditions with S t = 0.4.b) Phase space trajectories correspond to those shown in (a), arrows indicate the trajectories' starting point.c) Vorticity field at two different times separated approximately by one cycle.

FIG. 4
FIG. 4. a) Rotational angle as function of time for different initial conditions with S t = 0.84.b) Phase space trajectories correspond to those shown in (a), the arrows indicate the trajectories' starting point.c) Vorticity field at two different times separated approximately by an oscillation of the elliptical cylinder.

FIG. 5
FIG. 5. a) Rotational angle as function of time for different initial conditions for S t = 1.16.b) Phase space trajectories correspond to those shown in (a), arrows indicate the trajectories' starting point.c) Vorticity field at two different times separated approximately by half an oscillation of the elliptical cylinder.

FIG. 6
FIG. 6. a) Rotational angle as function of time for different initial conditions for S t = 1.68.b) Phase space trajectories correspond to those shown in (a), arrows indicate trajectories' starting point.c) Vorticity field at two different times separated approximately by an oscillation of the elliptical cylinder.