Abstract
Hydrodynamic forces on a structure are the manifestation of fluid-structure interaction. Since this interaction is nonlinear, these forces consist of various frequencies: fundamental, harmonics, excitation, sum, and difference of these frequencies. To analyze this phenomenon, we perform numerical simulations of the flow past stationary and oscillating cylinders at low Reynolds numbers. We compute the pressure, integrate it over the surface, and obtain the lift and drag coefficients for the two configurations: stationary and transversely oscillating cylinders. Higher-order spectral analysis is performed to investigate the nonlinear interaction between the forces. We confirmed and investigated the quadratic coupling between the lift and drag coefficients and their phase relationship. We identify additional frequencies and their corresponding energy present in the flow field that appear as the manifestation of quadratic nonlinear interaction.
1. Introduction
The fluid-structure interaction has its significance in flow physics and industrial applications. The flow behind a circular cylinder has become the canonical problem for studying such external separated flows [1–4]. When flow passes over a bluff body, an organized and periodic motion of a regular array of concentrated vorticity, known as the von Kármán vortex street, appears in the wake of the bluff body. This vortex shedding exerts oscillatory forces on the body, which are often decomposed into drag and lift components along the freestream and crossflow directions, respectively. If the body is capable of flexing or moving rigidly, these forces can cause it to oscillate and the classical vortex-induced vibration (VIV) problem takes place.
Many experimental and numerical studies have been performed to understand, model, and predict the phenomenon of VIV for fixed, excited, and elastic cylinders. The problem of a vibrating cylinder due to exerted forces goes back to the work of Strouhal [5] in the area of aeroacoustic and to the work of Rayleigh [6] on the oscillations of violin strings subject to the incoming wind. The first significant contribution to this problem is credited to Bishop and Hassan [1] who experimentally studied the flow past an externally oscillated cylinder over a Reynolds number range of 5,850 to 10,800 within the lock-in frequency range. The lock-in frequency range is the bandwidth where the lift frequency is entrained by the oscillating frequency of the cylinder. They reported that this bandwidth is a function of the Reynolds number and the amplitude of motion.
Williamson and Roshko [7] performed an experimental study of the flow over a cylinder undergoing simple harmonic motion in a low Reynolds number range, 300 to 1,000. They used aluminum particles on the fluid surface for visualization. They observed various vortex patterns and constructed a map of vortex synchronization regions using coordinates as a function of the nondimensional wavelength of cylinder motion (). They denoted certain vortex patterns by symbols, such as “P + S”, “S”, and “2P”. Here “S” denotes a single vortex, and “P” signifies a pair of vortices of opposite signs. Thus, the vortex pattern “P + S” is one in which a pair and a single vortex are shed during each oscillation cycle. Similarly, two vortex pairs are shed per cycle in a “2P” wake pattern, while two singlets are shed in a “2S” wake pattern.
Krishnamoorthy et al. [8] investigated the wake patterns of an oscillating cylinder at Reynolds number of 1,500. They varied the frequency of oscillation from one-third to three times the natural shedding frequency with a nondimensional amplitude of 0.22. They observed abrupt switching between some vortex patterns and a phase jump in the vortex shedding relative to the cylinder displacement.
Recent numerical studies on the flow around stationary and oscillating cylinders include the direct numerical simulations by Dong and Karniadakis [13], where they have simulated a cylinder excited at a nondimensional amplitude of 0.3 with at Re = 10,000. They used the spectral element method with polynomial orders ranging from 5 to 8 and 300 million degrees of freedom to resolve the flow at this high Reynolds number. The computational domain employed for the simulations extended 20 diameters upstream and 50 diameters downstream, with 40 diameters in the crossflow directions. The problem was solved in a coordinate system that is attached to the cylinder.
Al-Jamal and Dalton [14] performed a two-dimensional LES to compare the flow past a forced and an elastically vibrating cylinder within the lock-in region at Reynolds number of 8,000. They used a vorticity formulation instead of primitive-variable formulation with second-order central-difference scheme for all spatial derivatives. The problem was solved in an absolute coordinate system over an O-grid with a radius of five cylinder diameters. For the elastically mounted cylinder, spectral and complex demodulation analyses show that the phase angle between the lift and motion is strongly modulated, which is remarkably different from what is typically observed in an externally excited motion.
Kim and Williams [15] experimentally studied the nonlinear coupling of fluctuating drag and lift on a forced oscillating cylinder at = 15,200. They considered two cases of cylinder oscillations, crossflow and inline, at an oscillating frequency different from the von Kármán vortex shedding frequency. The nonlinear coupling between the von Kármán vortex shedding and the forcing frequency was investigated using higher-order spectral tools. Kim and Williams [15] reported quadratic nonlinear interaction between the vortex shedding modes and the oscillation frequency which produced sum and difference modes in the lift and drag coefficients spectra. However, their study is limited to small amplitudes of vibration.
Marzouk and Nayfeh [16] performed two-dimensional unsteady RANS simulations of the flow around stationary and crossflow oscillating cylinders with different frequencies and amplitudes in the lock-in regime at Re = 500. The problem was solved on a curvilinear coordinate system on an O-grid mesh with a radius of 25 cylinder diameters. They employed the artificial compressibility method to solve the governing equations with second-order accurate schemes in both time and space: central differences for the diffusion term and an upwind scheme for the convective terms. The one-equation Spalart-Almaras model is applied to represent the unresolved scales in the flow. They observed discontinuities between the frequency response and phase angle of the lift relative to the displacement and identified two distinctive modes, namely, a low-lift mode and a high-lift mode. In the low-lift mode, the lift and displacement are out-of-phase while, in the high-lift mode, they are in-phase. Moreover, they reported that the motion amplitude must exceed a certain threshold for these discontinuities to take place. For a comprehensive review of the VIV phenomenon, the reader is referred to Williamson and Govardhan [17].
Nayfeh et al. [18] investigated two wake-oscillator models to represent the lift, namely, the van der Pol and Rayleigh oscillators. Using a higher-order spectral moments analysis, they found that the phase angle between the lift components at and is around . Based on this finding, they concluded that the van der Pol oscillator is more suitable for modeling the steady-state lift coefficient; that is,The angular frequency in (1) is related (but not equal) to the angular shedding frequency and the parameters and represent the linear and nonlinear damping coefficients, respectively. The values of the parameters in (1) depend on the Reynolds number. Nayfeh et al. [19] modified the computational fluid dynamics (CFD) based analytical model [18] to include the transient flow over the cylinder. Marzouk et al. [20] added a Duffing-type cubic term to the van der Pol oscillator model of the lift to match the phase between different harmonics. Later, Akhtar et al. [21, 22] extended the model to a more general shape of elliptic cylinders with different eccentricities. Unlike the lift oscillator model, the drag coefficient is modeled as a function of the or depending upon their phases, such thatEquation (2) represents the coupling between the lift and drag coefficients.
Modeling of the lift and drag coefficients is complex and requires a thorough understanding of flow physics even in the stationary case. This fluid-structure interaction becomes more complex in the nonstationary case where the cylinder is either forced to oscillate or vibrates under the influence of hydrodynamic forces. Due to inherent nonlinearity in this phenomenon, one would expect interaction of multiple frequencies in the dynamical system. For three-dimensional flows, which is often the case in real applications, the hydrodynamic forces are affected by the turbulent structures formed in the wake and can not be modeled without including the effects of these structures.
The objective of the current study is to quantify the nonlinear coupling between the hydrodynamic forces on stationary and oscillating circulars [23]. In the current study, we perform numerical simulations of the flow past two- and three-dimensional flows past stationary and oscillating cylinders at and . Higher-order spectral techniques, such as auto-power spectrum, cross-bispectrum, and cross-bicoherence, are employed to identify the quadratic coupling between the lift and drag coefficients acting on the cylinder. In this study, modes refer to various frequencies that appear as the manifestation of quadratic coupling.
The manuscript is organized as follows. Section 2 presents the numerical methodology to solve the incompressible Navier-Stokes equations. The numerical results are then validated and verified. In Section 3, we discuss the Accelerating Reference Frame (ARF) methodology employed to simulate the moving structure, that is, oscillating cylinder. Again, the numerical results are compared with the experimental results for an oscillating cylinder with different frequencies and amplitudes. We define cross-bispectrum and cross-bicoherence in Section 4. In Section 5, we present the numerical results and analyze the nonlinear interaction for different flow configurations. We identify quadratic coupling between the lift and drag forces. This coupling provides an insight of the flow physics and is useful in developing reduced-order models for the hydrodynamic forces on the structure and computing the phase relationship between the forces.
2. Methodology and Numerical Scheme
The following section outlines the computational technique employed in the current study. Validation of the solver is performed by comparing the results with existing experimental and numerical data.
2.1. Numerical Methodology
The flow in the current problem is governed by the incompressible continuity and momentum (Navier-Stokes) equations, which can be represented as follows.
Continuity Equation
Momentum Equationwhere ; represents the Cartesian velocity components ; is the pressure; is the fluid density; and is the fluid kinematic viscosity. In order to broaden the spectrum of the application, we employ body fitted curvilinear coordinates for the governing equations. Thus, the irregular physical domain is transformed into a regular computational space. Therefore, the time-dependent incompressible continuity and Navier-Stokes equations are solved in a generalized coordinate system . Furthermore, we nondimensionalize the governing equations with the cylinder diameter (D) as the characteristic length and freestream velocity () as the characteristic velocity. The flow Reynolds number is defined as , where is the kinematic viscosity. The nondimensional governing equations in the generalized coordinates are written in a strong-conservative form aswhere the flux () is defined asIn (6), is the inverse of the Jacobian or the volume of the cell; is the volume flux normal to the surface of constant ; and is the “mesh skewness tensor.”
The governing equations are solved on a nonstaggered grid topology [24]. A second-order central-difference scheme is used for all spatial derivatives except for convective terms where a QUICK scheme is employed [25]. Temporal advancement is performed using semi-implicit predictor-corrector scheme, where an intermediate velocity is calculated at the predictor step and updated at the corrector step by satisfying the pressure-Poisson equation.
In this study, an “O”-type grid is employed to simulate the flow over a circular cylinder as shown in Figure 1. Inflow and outflow boundary conditions at the domain boundaries are simulated using Dirichlet and Neumann boundary conditions, and no-slip no-penetration boundary conditions are applied at the cylinder surface.
The computational domain is decomposed using two-dimensional decomposition, which allows us to take advantage of scalable parallel computers (e.g., distributed-memory clusters, message-passing programming model). In the domain decomposition technique, each processor gets a “slice” of the grid, as shown in Figure 2. The grid is distributed in the circumferential () and the longitudinal () directions and periodic boundary conditions are applied in both directions. This technique allows for a simple way to implement boundary conditions and keeps an equal load distribution for each processor.
(a)
(b)
The fluid force on the cylinder is the manifestation of the pressure and shear stresses acting on the surface of the cylinder, which can be decomposed into two components, namely, lift and drag forces. Lift and drag coefficients are expressed in terms of the pressure and shear stresses as follows:where is the spanwise vorticity component on the cylinder surface, is the axial length of cylinder, and is the pressure acting on the cylinder surface. It is important to note that all the flow quantities and governing equations are nondimensional.
2.2. Validation
To validate our numerical results, we perform two- and three-dimensional numerical simulations of the flow past a stationary circular cylinder at Reynolds numbers of 525 and 1,000. For the flow past a circular cylinder, the wake starts to get turbulent after Re = 500 and becomes fully turbulent at Re = 1,000. We compute the lift and drag coefficients using (7a) and (7b) and compare physical parameters, such as the Strouhal number () and the mean-drag coefficient , with other experimental and numerical results. We also indicate the number of processors () used in each simulation, where and are the number of processors in the - and -directions, respectively, in the domain decomposition topology.
2.2.1.
We perform two- and three-dimensional direct simulations of the flow over a stationary circular cylinder at = 525 with and grid points, respectively. In the simulation, the outer domain is while nondimensional spanwise length () is for the three-dimensional simulation. We use and processors for two- and three-dimensional simulations, respectively. A nondimensional time step of is selected to keep the CFL (Courant-Friedrich-Levy) number 0.2 for stable temporal advancement. The Strouhal number and the mean-drag coefficient are compared with the experimental results of Roshko [3], Weiselsberger [9], and numerical results of Mittal and Balachandar [10]. Our numerical results are tabulated in Table 1 and show a good comparison with the results of other studies. For the three-dimensional case, we compare various vortex structures present in the flow with the structures obtained in the numerical study of Mittal and Balachandar [10]. Figure 3 shows isosurfaces of the instantaneous absolute vorticity of the flow, which allows us to identify three-dimensionality in the flow structures. At the wake of the cylinder, we observe clockwise-rotating () and counter-clockwise-rotating () vortices in the von Kármán vortex street. In addition to the vortex street, we also observe hairpin vortices and the ribs (R) due to three-dimensional effects.
2.2.2.
Similar to the , we simulate this Reynolds number past a stationary cylinder using two- and three-dimensional grid structure. For the two-dimensional grid, we employ a mesh with a domain size of 30D. The grid is distributed over 8 processors such that each processor has a work load of grid points. Simulations are initiated with an initial guess for the velocity field and the flow is allowed to reach statistically stationary state under the prescribed boundary conditions. Figure 4 shows that the simulations were able to successfully capture the vortex shedding of the cylinder as we observe in the discrete vortices with opposite signs shed from the cylinder.
For the three-dimensional calculation, we performed a DNS on a grid size of distributed over 64 () processors with a computational domain of and spanwise length of that matches that of the numerical study of Evangelinos and Karniadakis [12]. The load per processor was grid points. A spanwise grid spacing of , inspired by the LES study of Kravchenko and Moin [26] at , is applied, which is fine enough to capture the fine length scales generated along the cylinder span. Simulations were successful in capturing the vortical structures such as ribs and hairpin vortices, which exist at this Reynolds number, Figure 5. Moreover, as observed in Figure 5, the wake is fully turbulent and the vortex shedding is no longer ordered.
(a)
(b)
(c)
Table 2 tabulates the mean-drag coefficient and the Strouhal number for the two- and three-dimensional flows at = 1,000 and those of Norberg [11] (experimental) and Evangelinos and Karniadakis [12] (numerical). The numerical variation is less than 10% for both the experimental and numerical results.
Figure 6 shows the 1D energy spectrum at in the flow. In addition to the peak at the shedding frequency 0.20, the power spectrum exhibits a law in the inertial range, which extends about half a decade in wave number. This proves the capability of the code in capturing turbulence.
2.3. Verification
In order to ensure the sufficiency of the grid resolution and domain size, we conduct grid and domain independence studies of our domain. Grid and domain independence studies are performed for two- and three-dimensional grids for flow past a stationary cylinder at .
2.3.1. Grid Independence Study
For the two-dimensional geometry arrangement, we considered two grid resolutions at domain size of , where the first grid resolution consisted grid points and the second grid was , which is refined by 50% both in the - and - directions. Similar to the two-dimensional arrangement, we considered two grid resolutions for the three-dimensional simulations, and grid points in the directions. The computational domain size was fixed to with spanwise length of . Table 3 summarizes mean-drag coefficient and Strouhal frequency (St) of the different grids considered. The results show that the percentage difference of the mean-drag coefficient and St is less than 5% for both two- and three-dimensional grids, which clearly indicates that the flow in the vicinity of the structure is virtually grid independent. Thus, we choose the former grid size for numerical simulations.
2.3.2. Domain Independence Study
In order to establish domain independence, we increased the domain size to , that is, 60% larger domain, with grid points. The slight increase in the grid points in the -direction on the larger domain is chosen to maintain comparable grid spacing in the two grids. Similarly, for three-dimensional simulations, the computational domain is increased to , keeping the spanwise length constant (). We have not changed the domain size in the -direction because increasing the domain size in the spanwise direction changes the geometry of the flow structures and thus can not be considered a similar flow. Again, the grid points in the larger domain are adjusted to keep comparable grid spacing in the -direction making it grid points. Table 4 shows that a domain size of is sufficient for the current study with the maximum percentage difference less than 5% for the mean-drag coefficient and Strouhal number. Thus, for the current study a domain size of meshed using for 2D simulations and for 3D simulations is deemed sufficient.
3. Moving Boundary Method
In a fluid-structure interaction, a structure may be elastic or excited externally. The behavior of a structure undergoing vibrations depends mainly on the natural frequencies of the structure and vortex shedding. If the frequency of vortex shedding is close to the natural frequency of the structure, it might oscillate with high amplitudes, eventually leading to its damage or complete structural failure.
In order to simulate the motion of the structure, three basic methods can be employed for the numerical discretization of the flow problem, namely, the Immersed Boundary (IB) method, the Arbitrary Lagrangian-Eulerian (ALE) method, and the Accelerating Reference Frame (ARF) method. In the Immersed Boundary (IB) method, the flow is simulated with immersed boundaries on a grid that does not conform to the shape of these boundaries. The advantage of IB method appears in the grid generation, which does not require coordinate transformation at every increment of body motion. However, applying appropriate boundary conditions is not straightforward and resolving the boundary layer at high Reynolds numbers increases the grid-size requirement faster than a corresponding body-conformal grid. For detailed description and comparison of this method, readers are referred to Mittal and Iaccarino [27].
In Arbitrary Lagrangian-Eulerian (ALE) method, the mesh local to the structure is distorted continuously in time as the structure moves, and the boundary conditions on the body and in the far field are usually fixed in time. In terms of computational cost, the ALE method is more computationally expensive than the IB method due to the continuous remeshing and the temporal changes of the mesh interpolation functions. Unlike ALE, the mesh in the Accelerating Reference Frame (ARF) method is fixed to the structure and the momentum equations and corresponding boundary conditions are modified to allow for the motion. In this way, the computational overhead associated with the coordinate transformation at every time step can be avoided; however, the ARF method is limited in its application.
3.1. Accelerating Reference Frame (ARF)
In our problem of the flow past an oscillating cylinder, the ARF method is a suitable candidate for simulating moving boundaries. The momentum equations can be directly coupled with the cylinder motion by adding a frame acceleration term and modifying the boundary conditions to accommodate the moving reference frame [28]. Thus, the governing equations describing the relative motion of incompressible fluid in the reference frame attached to the structure are given as follows.
Continuity Equation
Momentum Equationwhere . In an inertial reference frame, the structure moves in response to the force per unit length exerted on it by the fluid according towhere is the mass per unit length, is the natural frequency, and is dimensionless structural damping and is the structure velocity.
At the domain boundary, the velocity boundary condition is modified to include the effect of moving body, such that where is the velocity in the inertial reference frame. On the structure surface, the velocity boundary condition is typically . Pressure boundary conditions are obtained by dotting the domain unit outward normal into (9), and using the vector identity . The result isIn the far field boundary, , while on the surface and . Similar methodology has also been used for free vibrations on the cylinder [29, 30].
3.2. Forced Oscillations
In the current study, the equation governing the nondimensional displacement (forced-oscillation) of the cylinder is modeled by a simple harmonic motion in the inline and crossflow directions as follows:where and are the nondimensional amplitudes in the inline and crossflow directions, respectively, and is the nondimensional excitation frequency. Therefore, the CFD flow solver is capable of simulating single- and two-degree-of-freedom motions, as shown in Figure 7. Most of the experimental and numerical studies have ignored the streamwise component of the motion (inline motion) since the crossflow in VIV problems is much more influential, with the exception of cases where the density of the cylinder is small relative to the fluid density [31]. Moreover, this assumption reduces the parameters in the simulations and enables better understanding of the effect of this component alone.
To validate the moving boundary feature in the parallel CFD solver, we simulate different flow configurations for a stationary and crossflow oscillating cylinders as shown in Table 5. We simulate the flow past a stationary cylinder and compute the shedding frequency . The excitation frequencies for the crossflow oscillating cylinders (Cases B–H) are nondimensionalized with respect to . The results obtained are compared with the experimental results of Koopmann [32]. In his experiment at , Koopmann [32] identified the lock-in region bounded by a curve in frequency-amplitude plane. From our numerical simulations, we also identify the lock-in and lock-out regions for each case and plot them in Figure 8. For lower nondimensionalized amplitude of oscillation, we observe lock-in only when is close to , that is, Case B. On the other hand, Cases A and C correspond to the nonsynchronous regime. In Cases E–H, the amplitude of oscillation is increased to 0.5 and frequency ratio is varied from 0.6 to 1.2 with an increment of 0.2. Case E is nonsynchronous whereas, in all other cases, lock-in is observed. From Figure 8, it is evident that the flow regimes are correctly predicted for crossflow oscillations thus validating the solver. Further details on the subject can be found in Akhtar [33] and Akhtar et al. [34].
4. Higher-Order Spectral Analysis
The fluid-structure interaction is a highly nonlinear phenomenon. One of the important properties of the nonlinear systems is that different frequency components in the system can interact with each other and generate new frequencies. These new components typically appear as the sum and the difference of combining frequencies. The phase of the frequency provides the information if it is a result of a nonlinear interaction of other frequencies. Analysis of the time series of any fluctuating physical quantity, for example, lift and drag coefficients, starts by determining the frequency distribution of power of the fluctuation, which can be found by calculating the auto-power spectrum of the time series. The power spectrum helps in identifying the dominant frequencies of the domain and their harmonics, but it does not provide information about whether a nonlinear wave-wave interaction exists. For any wave-wave interaction to be considered nonlinear, it must satisfy the following selection rules or resonance conditions: and [35], where is the wave frequency and is the wave number. However, there is always the possibility that the peak at is not a result of the nonlinear interaction waves, if is a normal mode of the system. In order to distinguish both cases, we perform higher-order spectral analysis techniques to measure the degree of phase coherence between modes [36], namely, the cross-bispectrum . The cross-bispectrum can be expressed aswhere , , and denotes Ensemble averaging. Also, is the complex conjugate of . In the above formulation, and are the time signals under consideration. From the definition of the cross-bispectrum, we see that it measures the statistical dependence of three waves. That is, if , , and are independent frequencies, then the statistical averaging of the cross-bispectrum will return a value close to zero due to the random phase mixing effect [36]. Whereas if the three waves are nonlinearly (quadratically) coupled, the Ensemble average of the bispectrum will not return a zero value.
In the current study, we investigate the nonlinear (quadratic) coupling between the time histories of the lift, , and drag, , coefficients of stationary and oscillating cylinders. The value of bispectrum in the frequency plane determines the degree of quadratic coupling in the system, where near zero value indicates independent random phases between , , and , while higher bispectrum value identifies the quadratic coupling between the two signals. In order to remove the cross-bispectrum dependence on the frequency components amplitudes, we use the cross-bicoherence, which is a normalization of the cross-bispectrum by the amplitudes of the spectral components amplitudes [37]. The cross-bicoherence is defined as [36] Based on Schwarz’s inequality, the value of the cross-bicoherence is bounded within . A value of cross-bicoherence close to one indicates strong quadratic phase coupling, a value near zero indicates the absence of coupling, and a value between zero and one indicates partially coupled frequencies.
From the symmetry properties of the cross-bispectrum and cross-bicoherence, the calculation of these quantities is limited to two regions: Region S, which defines the sum interaction region and is calculated on and , and Region D, which defines the difference interaction region and is calculated on and , as shown in Figure 9.
5. Results and Discussion
Simulations are performed for six different cases of stationary and crossflow oscillating cylinders at and 1,000 oscillating at with varying amplitudes, as summarized in Table 6. In the current study, only the crossflow oscillation case is considered due to its pronounced physical effect on the dynamics of the system as discussed earlier. Numerical simulations are initiated with an initial guess for the velocity field and the simulations are run long enough for the flow to develop under the prescribed conditions. The simulations continued running after the flow reached statistically stationary behavior for enough time to establish a long enough time histories of the lift and drag coefficients of the flow. These histories are used to calculate the different statistical quantities. The data sampling frequency is 50 Hz for all cases under consideration.
Time histories of the lift and drag coefficients of the cylinder are divided into periodic segments, which are used to calculate the different statistical variables, that is, power spectrum, cross-bispectrum, and cross-bicoherence, and Ensemble averaging is used to reduce statistical variance of the results. The power spectra presented here are calculated using three time segments with 4096 samples each. For higher-order spectral moments, the number of realizations is increased 12 with 1024 data points for each realization.
First, we compute the auto-power spectrum of the lift coefficient for each case as shown in Figure 10. In Case J (stationary at ), we observe the most dominant peak occurs at the vortex shedding frequency . We also observe peaks at the odd harmonics of the shedding frequency, that is, ; however there is a logarithmic decrease in the relative amplitudes of these harmonics as compared to the shedding frequency. On the other hand, as the cylinder starts to oscillate, additional peaks appear in the power spectrum in addition to the shedding frequency and its harmonics observed in Case J. The power spectrum of Case K ( and ), Figure 10(b), shows a clear peak at the excitation frequency of the cylinder apart from the vortex shedding frequency. Moreover, we observe additional peaks resulting from the interaction between excitation and shedding frequencies; for example, we observe peaks at , and so on. As we increase the amplitude of oscillation in Case L, Figure 10(c), the power contained in the excitation frequency increases and more energy is contained within shedding and excitation bandwidth. Similar to Case J, additional peaks appear as a result of the interaction between the excitation frequency and the shedding frequency and its harmonics, Figure 10(c). As the Reynolds number increases to 1000 for Cases M, N, and P, the shedding frequency increases to 0.238 for the two-dimensional cases, M and N, and to 0.21 for the three-dimensional case P. Similar to the cases, odd harmonics are also present in the lift spectrum, and as the cylinder is excited, Case N, we observe similar trend as in Case K, but with more energy, whereas the power spectrum of Case P, which is a three-dimensional flow past a stationary cylinder, shows the effect of turbulent structures in the flow manifested in the noise in the lift spectrum.
(a) Case J
(b) Case K
(c) Case L
(d) Case M
(e) Case N
(f) Case P
Figure 11 shows the powers spectra of the drag coefficient for different cases. Even harmonics of the shedding frequency () appear in the drag coefficient.
(a) Case J
(b) Case K
(c) Case L
(d) Case M
(e) Case N
(f) Case P
In order to identify the quadratic coupling between lift and drag coefficients, we perform the cross-bispectrum/bicoherence analysis of the lift and drag signals. Figure 12 presents contours of the cross-bispectrum () of the different flow cases considered, where we observe high levels at the interaction of different frequencies as their sum and difference. For Case J (2D at ), the magnitudes of the peaks appearing in the cross-bicoherence at () and () are summarized below, along with the coupling relations between peaks in the spectra.(i), ,(ii), ,
(a) Case J
(b) Case K
(c) Case L
(d) Case M
(e) Case N
(f) Case P
which is consistent with the findings of Kim and Williams [15]. In Cases K and L, we observe additional contours corresponding to the quadratic coupling of the dominant frequencies in their lift and drag spectra. We compute cross-bicoherence to identify the strength of quadratic coupling between lift and drag. Some of the levels of dominant interaction are as follows:(i), ,(ii), ,(iii), ,(iv), ,(v), .
High values of indicate strong quadratic coupling. Similarly, as the Reynolds number increases in Case N, we also observe high contour levels. Due to low Reynolds number and two-dimensional geometry, we obtain high levels of cross-bicoherence for all the cases. It is however interesting to study the effect of third-dimensionality by comparing Case M with Case P. In Figure 12, we observe a larger coupling region for 3D because of sidebands in the dominant frequencies. Multiple frequencies present in the lift spectrum due to third-dimensionality interact to produce wider range of frequencies in the drag spectrum. However, it would be interesting to analyze the effect of third-dimensionality for the crossflow motions which would be the future direction for our analysis.
6. Conclusions
We numerically simulated the flow past stationary and oscillating cylinders and compute the hydrodynamic forces on it. We performed higher-order spectral analysis of the lift and drag coefficients and identified different frequencies generated due to nonlinear interaction. Bispectrum/bicoherence analyses show the energy content of the quadratic coupling between the modes present in the lift and drag. In the future, we would perform 3D simulations of the flow past oscillating cylinders and analyze the effect of third-dimensionality on the nonlinear interaction between different frequencies and how it differs from the two-dimensional flows.
Conflicts of Interest
All the authors declare that there are no conflicts of interest regarding the publication of this article.