Higher-Order Spectral Analysis to Identify Quadratic Nonlinearities in Fluid-Structure Interaction

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.


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][2][3][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.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 highlift 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  13 between the lift components at   and 3  is around 90 ∘ .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   = 2  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 that Equation ( 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 twoand three-dimensional flows past stationary and oscillating cylinders at Re = 525 and 1000.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 crossbicoherence 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.

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.

Continuity Equation
Momentum Equation where ,  = 1, 2, 3;   represents the Cartesian velocity components (, V, );  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 Re  =  ∞ /], where ] is the kinematic viscosity.The nondimensional governing equations in the generalized coordinates are written in a strong-conservative form as where the flux (  ) is defined as In (6),  −1 = det(  /  ) is the inverse of the Jacobian or the volume of the cell;   =  −1 (  /  )  is the volume flux normal to the surface of constant   ; and   =  −1 (  /   )(  /  ) 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 twodimensional 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.
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.

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.
Re D = 525.We perform two-and three-dimensional direct simulations of the flow over a stationary circular cylinder at Re  = 525 with 192×256 and 128×192×32 grid points, respectively.In the simulation, the outer domain is 30 while nondimensional spanwise length (  ) is  for the threedimensional simulation.We use 8 × 1 and 8 × 2 processors for two-and three-dimensional simulations, respectively.A nondimensional time step of Δ = 2 × 10 −3 is selected to keep the CFL (Courant-Friedrich-Levy) number ≈ 0.2 for stable temporal advancement.The Strouhal number and the meandrag 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 threedimensionality in the flow structures.At the wake of the cylinder, we observe clockwise-rotating (CR) and counterclockwise-rotating (CCR) 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.Re D = 1,000.Similar to the Re  = 525, we simulate this Reynolds number past a stationary cylinder using two-and three-dimensional grid structure.For the twodimensional grid, we employ a 192 × 256 mesh with a domain size of 30D.The grid is distributed over 8 processors such that each processor has a work load of 192 × 32 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 192×256×192 distributed over 64 (8×8) processors with a computational domain of 30 and spanwise length of 4 that matches that of the numerical study of Evangelinos and Karniadakis [12].The load per processor was 192 × 32 × 24 grid points.A spanwise grid spacing of Δ = 0.06545, inspired by the LES study of Kravchenko and Moin [26] at Re  = 3,900, 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.
Table 2 tabulates the mean-drag coefficient and the Strouhal number for the two-and three-dimensional flows at Re  = 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 (/, /) = (1.0,0.5) in the flow.In addition to the peak at the shedding frequency 0.20, the power spectrum exhibits a −5/3 law in the inertial range, which extends about half a decade in wave number.This proves the capability of the code in capturing turbulence.

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 threedimensional grids for flow past a stationary cylinder at Re  = 525.

Grid Independence Study.
For the two-dimensional geometry arrangement, we considered two grid resolutions at domain size of 30, where the first grid resolution consisted 192×256 grid points and the second grid was 288×384, 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, 128×192× 32 and 192×256×48 grid points in the (, , ) directions.The computational domain size was fixed to 30 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 meandrag coefficient and St is less than 5% for both two-and threedimensional 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.

Domain Independence Study.
In order to establish domain independence, we increased the domain size to 50, that is, 60% larger domain, with 216 × 256 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 50, 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 216 × 256 × 32 grid points.Table 4 shows that a domain size of 30 is sufficient for the current study with the maximum percentage difference less than 5% for the meandrag coefficient and Strouhal number.Thus, for the current study a domain size of 30 meshed using 192 × 256 for 2D simulations and 128 × 192 × 32 for 3D simulations is deemed sufficient.

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.

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 Equation where ,  = 1, 2, 3.In an inertial reference frame, the structure moves in response to the force per unit length  exerted on it by the fluid according to where  is the mass per unit length,   is the natural frequency, and  is dimensionless structural damping and ẋ = k is the structure velocity.At the domain boundary, the velocity boundary condition is modified to include the effect of moving body, such that where u  is the velocity in the inertial reference frame.On the structure surface, the velocity boundary condition is typically   = 0. Pressure boundary conditions are obtained by dotting the domain unit outward normal  into (9), and using the vector identity ∇ 2  = ∇(∇ ⋅  − ∇ × ∇ × ).The result is In the far field boundary,   / = − ẍ  , while on the surface   / = 0 and   = 0. Similar methodology has also been used for free vibrations on the cylinder [29,30].

Forced Oscillations.
In the current study, the equation governing the nondimensional displacement (forcedoscillation) 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 twodegree-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 Re = 200, 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  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].

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 wavewave interaction exists.For any wave-wave interaction to be considered nonlinear, it must satisfy the following selection rules or resonance conditions:  1 +  2 =  3 and  1 +  2 =  3 [35], where   is the wave frequency and   is the wave number.However, there is always the possibility that the peak at  3 is not a result of the nonlinear interaction waves, if  3 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 ( 1 ,  2 ).
The cross-bispectrum can be expressed as where () = fft(()), () = fft(()), 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 crossbispectrum, we see that it measures the statistical dependence of three waves.That is, if  1 ,  2 , and  3 =  1 +  2 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 crossbispectrum 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 crossbicoherence is bounded within [0, 1].A value of crossbicoherence 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 0 ≤  2 ≤ 0.5  and  2 ≤  1 ≤   − 2 , and Region D, which defines the difference interaction region and is calculated on −  ≤  2 ≤ 0 and  2 ≤  1 ≤   , as shown in Figure 9.

Results and Discussion
Simulations are performed for six different cases of stationary and crossflow oscillating cylinders at Re = 525 and 1,000 oscillating at   = 0.16 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 Re  = 525), 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, 3  , 5  ; 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 (Re  = 525 and   = 0.07), 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 2  −   , 2  +   , 4  −   , 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).magnitudes of the peaks appearing in the cross-bicoherence at (  ,   ) and (3  , −  ) are summarized below, along with the coupling relations between peaks in the spectra.
(i)   +   → 2  ,  2  = 0.99, (ii) 3  −   → 2  ,  2  = 0.97, 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 crossbicoherence to identify the strength of quadratic coupling between lift and drag.Some of the levels of dominant interaction are as follows:  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 thirddimensionality 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.

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.

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.

Figure 1 :
Figure 1: A 2D layout of an "O"-type grid in the (, )-plane showing the inflow and outflow directions.

Figure 2 :
Figure 2: (a) A 2D layout of "O" grid distributed among 8 processors in the -direction.The grid is plotted only for the region of processor "1".(b) A 3D layout of the complete domain and the grid is plotted for only one processor, indicating the load per processor in a 24 (8 × 3) processor platform.The grid is plotted only for the region of processor "1".

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Spanwise vorticity contours of one vortex shedding cycle over a stationary cylinder at Re  = 1,000.

Figure 7 :Figure 8 :
Figure 7: A schematic of the possible stationary and moving boundary cases.

Figure 9 :
Figure 9: Symmetry properties of bispectrum reduce the region to sum (S) and difference (D).

Figure 10 :
Figure 10: Power spectra of lift coefficients   for the flow past a cylinder.