{'fraction_non_alphanumeric': 0.028469564990270722, 'fraction_numerical': 0.009703159076959242, 'mean_word_length': 4.670468009478673, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decomposition of wake dynamics in fluid-structure interaction via low-dimensional models

We present a dynamic decomposition analysis of the wake flow in fluid-structure interaction (FSI) systems under both laminar and turbulent flow conditions. Of particular interest is to provide the significance of low-dimensional wake flow features and their interaction dynamics to sustain the free vibration of a square cylinder at a relatively low mass ratio. To obtain the high-dimensional data, we employ a body-conforming variational fluid-structure interaction solver based on the recently developed partitioned iterative scheme and the dynamic subgrid-scale turbulence model for a moderate Reynolds number. The snapshot data from high-dimensional FSI simulations are projected to a low-dimensional subspace using the proper orthogonal decomposition (POD). We utilize each corresponding POD mode for detecting features: the vortex street, the shear layer and the near-wake bubble. We find that the vortex shedding modes contribute solely to the lift force, while the near-wake and shear layer modes play a dominating role to the drag force. We further examine the fundamental mechanism of this dynamical behavior and propose a force decomposition technique via low-dimensional approximation. We ascertain quantitatively that the shear layer feeds the vorticity flux to the wake vortices and the near-wake bubble during the wake-body synchronization. Based on the decomposition of wake dynamics, we suggest an interaction cycle for the frequency lock-in during the wake-body interaction, which provides the inter-relationship between the high amplitude motion and the dominating wake features. Through our investigation of wake-body synchronization below critical Re range, we discover that the bluff body can undergo a synchronized high-amplitude vibration due to flexibility-induced unsteadiness. The interaction cycle for the wake-synchronization is found to be valid for the turbulent wake flow.


Resonant dynamics of coupled fluid-structure system
Unsteady flows involving fluid-structure interactions are widespread in numerous engineering applications and their fundamental understanding poses serious challenges due to the richness and complexity of nonlinear coupled physics. Even a simple configuration of a coupled fluid-structure system can exhibit complex spatial-temporal dynamics and synchronization as functions of physical parameters and geometric variations. Synchronization is a general nonlinear physical phenomenon in fluid-structure systems whereby the coupled system has an intrinsic ability to lock to a preferred frequency and amplitude. For example, when a bluff body immersed in a cross-flow is flexible or mounted elastically, there exists a strong coupling between the bluff body and the vortices forming in its wake. In particular, as the natural frequency (f n ) of the bluff body approaches to the frequency of the wake system, typically the frequency of vortex shedding (f vs ), the wakebody frequency lock-in behavior is observed which plays a crucial role in establishing the synchronization. During this frequency lock-in, the bluff body experiences a large selflimiting vibrations (Khalak & Williamson 1999) and a dynamical equilibrium between the energy transfer and dissipation exists. This wake-body synchronization has been a major topic of research to understand the mechanism of this energy transfer and the sustenance of self-excited vibrations. In the present study, we consider a prismatic square geometry to understand the wake-body synchronization and to perform the decomposition of wake dynamics during the fluid-structure interaction.
The phenomenon of frequency lock-in is a major concern in offshore, marine and aeronautical engineering, whereby structures are designed to avoid the large-amplitude vibrations by selecting optimal system parameters (e.g., geometric dimensions, stiffness, damping) and/or installing active and passive devices to control the intensity of fluidstructure interaction. In particular, several studies have been conducted with the purpose of controlling the wake-body interaction via passive and active devices (Law & Jaiman 2017;Guan et al. 2017;Narendran et al. 2018) with the physical insight based on the reliance of frequency lock-in on the large-scale features of the wake. In fact, these studies were found to be remarkably successful in suppressing large-amplitude motion of the body by avoiding the interaction between the major organized features of the wake. However, the mechanism of the interactions among the wake features and their impact on the free motion of the bluff body is not properly explained. Moreover, the available experimental and numerical data can be used to provide a deeper understanding and a new insight into the kinematics and dynamics of synchronized wake-body interaction. This paper aims at explaining how different organized flow features (i.e., near-wake structures) amplify the bluff body motion and sustain the energy transfer from the fluid flow to the vibrating body. Specifically, we examine the formation of the dominant coherent structures and their nonlinear interactions during the wake-body synchronization.
The vortex shedding pattern is undoubtedly the most prominent wake feature behind a bluff body. It is omnipresent in almost all of the separated wake flows and has been studied extensively in the literature. This primary wake feature begins at a much lower Re: for example in a circular cylinder wake, at Re ≈ 49, exhibits a classical Kármán vortex street and develops the three-dimensional vorticity patterns when Re 190. In addition to the vortex street, a free shear layer (not attached to a solid surface) is an important dynamical feature that represents a separating high-gradient layer behind a bluff body and it arises between the higher free stream velocity and the smaller velocity occurring in the wake region. The shear layer behaves likes a perturbed vortex sheet and is highly sensitive and unstable to small disturbances, giving rise to alternating thickening and thinning of the vortex sheet. The characteristic vortex structures develop when the thickening of shear layer occurs. For the unsteady 2D regime (49 Re 190 for a circular cylinder) the roll-up of the shear layers with the formation of the vortex street can be observed (Williamson 1996). These shear layers are predominantly elongated in the streamwise direction and have a high-gradient in the cross-flow direction. Behind a moving or stationary bluff body, the region of a recirculating region with the rotational flow is present due to the fluid viscosity. Owing to nonlinear flow separation and turbulence, complex interactions occur in the mean recirculation region, which is also referred to as the near-wake bubble. Several previous studies have explored the dynamical features inside the wake region (with the vortex shedding and the shear layer) using experimental (Cantwell & Coles 1983), numerical (Braza et al. 1986) and both (Bearman 1997;Dong et al. 2006) techniques. In our present analysis, we consider the near-wake bubble as a distinct feature from the vortex street and the shear layer. The near-wake region accounts for the complex interactions of the mean circulation region, which can be considered as a general feature and can be identified separately from the other two features. Hence, we divide the wake into three dominant organized coherent structures: the vortex street, the shear layer, and the near-wake bubble. These organized features have an intrinsic dynamics of their own and influence each other in a nonlinear manner over a wide range of space and time scales. A primary goal of this paper is to employ low-dimensional models to extract the organized wake features and to examine their roles during the wake-body synchronization.

Low-dimensional models for wake features
To extract the large-scale organized/coherent wake features, it is required to decompose the dynamic flow fields by scales into different constituent kinematical regions. The concept of decomposition by scale has been prevalent in many fluid dynamics research ranging from a low-dimensional projection of flow field to the turbulence modeling by ensemble averaging, temporal or spatial averaging. A general decomposition technique can be considered to separate the space-time data for representing different characteristics of the field. For example, the proper orthogonal decomposition (POD) extracts the most energetic modes in an optimal way and provides structural information from the wake data. The POD is a popular method for constructing low-order modeling from the data Holmes (2012), and it is often referred to as the Karhunen-Loève expansion or the principal component analysis. The key idea behind the Karhunen-Loève expansion is to determine a low-dimensional affine subspace from the high-dimensional data while retaining the important dynamics of the full-order model. After the determination of the best approximating low-dimensional subspace, a Galerkin projection is employed to project the dynamics onto it. In this work, we will employ this low-dimensional subspace projection procedure for extracting the large-scale wake features from the highdimensional flow dynamics data.
In the context of the present study, the POD-Galerkin projection method is quite attractive to capture the synchronized dynamics such as the vortex shedding and the near-wake interactions (Rempfer 2003;Noack et al. 2003). In addition, it has been the prominent empirical model reduction technique incorporated for the standard flow around a stationary circular cylinder for the past few decades. For example, in the one of pioneering study, Deane et al. (1991) reproduced the flow dynamics of the laminar wake by employing merely an eight-dimensional model, which was further generalized to generate reduced spaces for 3D velocity field by Ma & Karniadakis (2002) using direct numerical simulation (DNS) data. In general, the empirical POD-Galerkin models are capable of reconstructing the reference dynamics with higher accuracy than the standalone mathematical or physical Galerkin methods, while capturing the physically most significant modes (Noack et al. 2003). With regard to the applications of the POD-Galerkin to bluff body wake flows, these modes correspond to the organized wake features such as the vortex street, the shear layer and the near-wake bubble. Although there exist a significant body of works on the wake modes for a stationary cylinder, they have not been examined in the context of wake-body interaction and the lock-in process. One of the contributions of the present study is to build some connections between the wake features and the lock-in process. During the lock-in/synchronization, the vibrating body undergoes a highly nonlinear-wake interaction with self-sustained oscillations.
In the early studies of POD application to fluid flows (Lumley 1967;Sirovich 1987), the dynamic flow field is reconstructed by a linear combination of the most significant modes. Hence, it has a considerable local error in the highly nonlinear regions of the organized wake motions and the evaluation of the projected nonlinear term has a direct dependence on the large dimension of the original system. This problem is mitigated to a certain extent by increasing the sampling frequency and/or refining the spatial discretization of the reference data. However, these temporal and spatial refinements increase the cost of model reduction without directly addressing the nonlinear nature of the flow. To introduce the nonlinearity, Petrov-Galerkin projections to the Navier-Stokes formulation or Koopman operators are incorporated in some studies . Instead of such explicit models, we employ the recently developed discrete empirical interpolation method (DEIM) (Chaturantabut & Sorensen 2009) for dynamical systems, which reconstructs the fields as a nonlinear combination of the POD modes. Apart from the POD basis subspace, the method relies on the additional POD basis to enrich the low-rank approximation of the nonlinear terms. In the POD-DEIM, a set of best points are selected using a greedy selection and the reconstruction is based on the time history of the field data of those points. This reduces the computational cost of the technique and further allows to capture nonlinearities during the reconstruction of highly nonlinear dynamic wake fields ).

Contributions and organization
For the past few decades, the studies on the low-dimensional decomposition of wake features have been primarily focused on flow past stationary bodies, particularly on a circular cylinder (Deane et al. 1991;Noack et al. 2003;Taira et al. 2017;. This may be due to the fact that the flow exhibits a diverse set of complex phenomena despite its simple geometry. However, very few studies (Liberge & Hamdouni 2010;Yao & Jaiman 2017) are found on unsteady fluid-structure interaction systems. Here, we provide a modal reduction study on the flow past a freely vibrating sharpcornered square cylinder with two-degree-of-freedom motions. We consider a configuration of a square cylinder for our numerical study of wake-body synchronization because: (i) this configuration has fixed and perfectly symmetric separation points at the leading sharp corners, (ii) entirely resonance-induced lock-in exists (Yao & Jaiman 2017). The physical investigation is general for any fluid-structure system involving the interaction dynamics of flexible structures with an unsteady wake-vortex system. We hypothesize that the solution space of wake-body interaction attracts a low-dimensional manifold, which allows building a set of basis vectors for a low-dimensional representation of the high-dimensional space. The low-dimensional subspace is constructed by means of the samples collected from the high-dimensional solutions via projection-based model order reduction. We utilize the linear and nonlinear POD-based reduced order reconstructions to understand the most significant features in the wake flow.
To extract the modes for the dynamics of wake-body synchronization, the POD in conjunction with the nonlinear POD-DEIM is applied to a set of samples collected from the full-order simulations. We exploit the obtained POD modes to answer the following interesting questions that are prevalent in the field of fluid mechanics: (i) How does the each of large-scale features contribute to the unsteady forces acting on the bluff body? (ii) How do the wake features interact when the structural frequency and vortex shedding frequency are locked-in, such that the vortices remain very energetic even the fluid has transferred energy to the structure? (iii) Will the wake and bluff body undergo synchronized motion below the critical Re due to the structural flexibility? (iv) What role does the wake turbulence play when we attempt to decompose the wake into its large-scale features? In relation to the first question, we quantify the force contribution from each wake feature mode to the streamwise (drag) and transverse (lift) forces and explain the observed variation. We further investigate the modal contribution of different wake features in the pre-lock-in, lock-in and post-lock-in regimes and propose a cycle explaining the sustenance of lock-in phenomena of the wake-body synchronization. We then explore the below critical Re flows to examine whether the bluff body and the wake can undergo synchronization via flexibility-induced unsteadiness. Finally, we apply POD decomposition to the three-dimensional flow at moderate Re = 22, 000, whereas the wake is fully turbulent after flow separation. A well-established dynamic large-eddy simulation (LES) is employed for generating full-order data for the turbulent wake. At this subcritical Reynolds number, we explore the role of turbulence during the reconstruction of flow-field data and extend the wake-body synchronization cycle to the turbulent flow.
The paper is structured as follows. In Section 2, we briefly review the full-order model (FOM) for the coupled fluid-structure system, which follows by the formulation of modal reduction via linear POD and nonlinear POD-DEIM. Section 3 discusses the problem setup and the mesh convergence study performed for the full-order analysis. In Section 4, the reduced order reconstruction of fluid fields using the linear and nonlinear POD methods are presented together with the analysis on the role of wake features in generating the forces. In Section 5, the mode energy contributions from different flow features under lock-in conditions are investigated and a self-sustaining cycle is proposed to explain the wake interaction with the bluff body. Section 6 investigates the wake-body synchronization phenomenon at below critical Re. Section 7 explores the application of modal decomposition for moderate Re flows and extends the proposed wake interaction cycle to the turbulent flow. Concluding remarks and the main results of the present study are provided in Section 8.

Numerical methodology
We first briefly summarize our high-dimensional full-order model to simulate the coupled fluid-body interaction using the incompressible Navier-Stokes equations and the rigid body dynamics.

Full-order model for fluid-body interaction
We employ a variational formulation based on the arbitrary Lagrangian-Eulerian (ALE) to solve the following coupled fluid-body system where subscripts f and s denotes the fluid and structural domains, and Ω f (t) and Ω s represent the fluid and solid domains, respectively. Here ρ f is the fluid density, u f and w are the fluid and mesh velocities at a spatial point x ∈ Ω f (t), and b f denotes the body force in the fluid domain. For the structural system, M, C and K are the mass, damping and stiffness matrices of the bluff body and F s is the external force acting on the body.
The function ϕ s (z 0 , t) maps the initial position vector of the center of mass (z 0 ) to its position at time t, and σ f is the Cauchy stress tensor for a Newtonian fluid given by: where p is the fluid pressure. In addition to the initial conditions and the standard Neumann/Dirichlet conditions, the coupled system incorporates the velocity and traction continuity conditions at the fluid-body interface Γ as follows: where n is the outer normal to the fluid-body interface. The above fluid-body interface conditions are satisfied by the body-conforming Eulerian-Lagrangian treatment, which provides accurate modeling of the boundary layer and the vorticity generation over a moving body. While Eqs. The weak variational form of Eq. (2.1) is discretized in space using equal-order isoparametric finite elements for the fluid velocity and pressure. In the present study, we utilize the nonlinear partitioned staggered procedure for the full-order simulations of fluid-structure interaction (Jaiman et al. 2016b). The motion of structure is driven by the traction forces exerted by the fluid flow, whereby the structural motion predicts the new interface position and the geometry changes for the moving fluid domain at each time step. The movements of the internal ALE fluid nodes are updated such that the mesh quality does not deteriorate as the motion of solid structure becomes large. To extract the transient flow characteristics, we solve the Navier-Stokes equations at discrete timesteps which lead to a sequence of linear systems of equations via Newton-Raphson type iterations. We employ the Conjugate Gradient (CG), with a diagonal preconditioner for the symmetric matrix arising from the pressure projection and the standard Generalized Minimal Residual (GMRES) solver based on the modified Gram-Schmidt orthogonalization for the non-symmetric velocity-pressure matrix. The above coupled variational formulation completes the presentation of the full-order model for the fluid-structure interaction.
From a model reduction viewpoint, the coupled system of the nonlinear differential equations for the fluid-body interaction can be written in the following form: where y is the column state vector describing the unknown degrees of freedom and F is a vector-valued function describing the spatially discretized governing equations.
In the present fluid-body system, the state vector comprises of the fluid velocity and the pressure as y = {u f , p} and the structural velocity involves the three translational degrees-of-freedom. For a discretized domain of m elements and n timesteps, the fullorder simulation outputs a high-fidelity data set y ∈ R m×n×q , where q is the number of variables in y. This dataset is extremely valuable to determine the instantaneous physics of the fluid-body system and to construct a low-order representation that preserves the behavior of the original system.

Low-order models
We now turn to the data-driven model reduction technique whose goal is to decompose the aforementioned high-dimensional data set into a set of low-dimensional modes. For that purpose, we can consider the decomposition of the nonlinear mapping F of Eq. (2.7) as F(y) = f + Ay + F (y), (2.8) where f denotes a constant column vector with m rows, A and F are the linear and nonlinear terms. For the ease of explanation, consider the solution vector y(x, t) ∈ R m×n comprising a single quantity of interest which have been determined at discretized locations x of the spatial domain and for a particular time t and the matrix operator A is an m × m matrix which captures the linear dynamics while F (y) is a nonlinear function of y. Using the projection-based model reduction, we can represent the state vector y by an element in a low-rank vector subspace spanned by the column vectors of The state vector y can be approximated by Vŷ, whereŷ is a reduced column vector with k entries. Since the columns of V are orthonormal (i.e., V T V = I), via the Galerkin projection onto the basis V, we get the following reduced dynamics: Next, we have to choose a suitable subspace for the mode decomposition. Using the reduced singular value decomposition (SVD), the above state vector y can be expressed as where the vectors v j are the POD modes of the matrix Y with rank k, W is an orthonormal matrix with n × k , and Σ is a k × k diagonal matrix with diagonal entries σ 1 σ 2 ... σ k 0. For any r k, the subspace spanned by {v 1 , ..., v r } provides an optimal representation of y in the subspace of dimension r using the SVD process. The total energy contained in each POD mode v j can be computed by the singular value σ 2 j . Note that V and W are the orthonormal eigenvectors of YY T and Y T Y, respectively.

Proper orthogonal decomposition
The POD method provides an algorithm to decompose a set of data into a minimal number of modes. We give a brief outline of this projection-based model reduction for the dynamical analysis of wake-body interaction. The general POD algorithm can be expressed as follows. Here, the eigenvectors of YY T are determined instead of performing the SVD. The algorithm adopted from Taira et al. (2017) is summarized in Algorithm 1.
ALGORITHM 1: Snapshot POD Input: Snapshots of spatial field expressed as Y(x, t) where Y ∈ R m×k (m -number of spatial points, k -number of snapshots.) Output: Significant r POD modes V = [v 1 , v 2 , ..., v r ] (i) Develop the fluctuation matrix by subtracting the mean: iii) Find the eigenvalues and eigenvectors of R by RV = ΛV (iv) Determine the number of required POD modes (r) using r j=1 λ j / m j=1 λ j ≈ 1.0, where λ j are the eigenvalues given by Λ.
The standard linear POD is almost in the same order expensive as the full-order analysis since it is using theỸỸ T matrix which has the size of m × m. In a typical timedependent flow analysis, it is unnecessary to generate m POD modes for comparison as the POD mode energy decays exponentially. Hence an alternative method the so called snapshot POD (Sirovich 1987) is applied to extract the most significant modes. In the snapshot method, the eigenvalue decomposition is performed onỸ TỸ ∈ R k×k which is significantly smaller thanỸỸ T as k m. Let the eigenvalues and eigenvectors ofỸ TỸ be given byỸ TỸ W = ΛW, (2.11) then using the relationship between the eigenvectors ofỸỸ T andỸ TỸ , a maximum of k significant POD modes can be extracted by Throughout the study, every POD decomposition will be performed via the snapshot POD method due to its low computational cost and memory usage. After extracting the significant POD modes, the constant and linear components of the instantaneous state vector can be recovered as a linear combination of the identified significant modes where r is the number of significant POD modes. The temporal coefficients of the linear combination are determined by the L 2 inner product . , . between the fluctuation matrix and the modes as followsŶ (2.14) This summarizes the process of POD by performing the SVD on the snapshots of the sampled solutions at certain timesteps. While the above POD-Galerkin process can reconstruct the linear term to the expected error threshold, the nonlinear term will not be reconstructed properly in the context of nonlinear incompressible flow which involves quadratic nonlinearity. The linear POD reconstruction requires a higher number of modes and/or a smaller sampling interval for snapshots to obtain the required local domain accuracy. In other words, the spatial and/or temporal discretizations of the POD method have to be so small such that the nonlinearities behave almost linearly. Subsequently, the POD reconstruction may result in a similar order of computational expense as full-order simulation. This issue can be handled by employing discrete empirical interpolation method (DEIM). The DEIM introduces the nonlinearity by supplementing an additional basis for a low-order representation of nonlinear terms. This gives rise to the reduction in the requirement of POD modes, hence decreasing the computational cost while capturing the nonlinear regions properly.

Discrete empirical interpolation method
To overcome the difficulty in the linear POD, Chaturantabut & Sorensen (2009) proposed the discrete empirical interpolation method to reconstruct the full-order variable as a nonlinear combination of the POD modes. The aim of DEIM is to design a low-order representation for the nonlinear terms by introducing an additional basis. Consider U as a basis generated from the leading l modes of the POD, which is attracted to a low-dimensional subspace. We can approximate the nonlinear term in Eq. (2.7) by the sequence of nonlinear snapshots as F (Vŷ(t)) ≈ Uĉ. The coefficientsĉ can be selected based on Algorithm 2, which relies on a greedy approximation of nonlinear function. In Algorithm 2,ρ and ℘ 1 denote the assigned value and the assigned index of max{|v 1 |}, and e ℘i = [0, ..., 0, 1, 0, ...0] T ∈ R m is the ℘ i th column of the identity matrix of size m × m. The accuracy of DEIM approximation depends on the error induced by the POD projection and the estimation of (P T U ) −1 . Further details of the DEIM process can be found in Chaturantabut & Sorensen (2009).
Here, a set of entries ℘ ⊂ {1, 2, ..., l} often called optimal (best) points are selected to determineĉ by the following relation ( 2.15) Assuming that F is a component-wise function P T F (Vŷ(t)) = F (P T Vŷ(t)), we can rewrite Eq. (2.7) as In the present work, we perform the nonlinear POD on the same fluctuation matrixỹ m×k without separating the linear and nonlinear components. Consider the approximation of y as a nonlinear combination of the POD modes: (2.17) The coefficients θ(t) are calculated by the conditions imposed by the POD-DEIM. While the POD modes are linearly independent, we can obtain a unique number of DEIM points if P T U matrix is invertible. By using just the ℘ rows of V and U , we can establish the following relationship: If the number of points used is higher than the number of significant modes, i.e. l > r, which is often the case, V ℘ becomes a rectangular matrix. This makes the coefficients θ(t) given by the gappy POD reconstruction (2.20) The solution to the least square problem (Eq. 2.20) gives the result The POD-DEIM provides a way to introduce nonlinearity to the POD reconstructions, however, due to this nonlinear behavior, it is not guaranteed to converge to the full-order results. In other words, the use of more POD modes or DEIM points does not assure an improvement in the result. Therefore, determining the optimal sizing of the lowdimensional representation is critical when using POD-DEIM for reconstruction. In the next section, we present the full-order model for generating high-dimensional data.

Problem set-up
In this section, we give an overview of full-order simulations for a freely vibrating structure immersed in a viscous incompressible fluid flow. Specifically, the focus of this section is to present numerical results on the flow past an elastically mounted square cylinder, whereby the cylinder is free to oscillate in the streamwise (X) and the transverse (Y ) directions. The mass and natural frequencies are identical in both X-and Ydirections. The translational flow-induced vibration of a cylinder is strongly influenced by the four key non-dimensional parameters, namely mass-ratio (m * ), Reynolds number (Re), reduced velocity (U r ), and critical damping ratio (ζ) defined as where M is the mass per unit length of the body, C and K are the damping and stiffness coefficients, respectively for an equivalent spring-mass-damper system of a vibrating structure, U ∞ and D denote the free-stream speed and the diameter of cylinder, respectively. The natural frequency of the body is given by f n = (1/2π) K/M and the mass of displaced fluid by the structure is m f = ρ f D 2 L c for a square crosssection, and L c denotes the span of the cylinder. In the above definitions, we make the isotropic assumption for the translational motion of the rigid body, i.e., the mass vector M = (m x , m y ) with m x = m y = M , the damping vector C = (c x , c y ) with c x = c y = C, the stiffness vector K = (k x , k y ) with k x = k y = K. The fluid loading is computed by integrating the surface traction considering the first layer of elements located on the cylinder surface. The instantaneous lift and drag force coefficients are evaluated as Here n x and n y are the Cartesian components of the unit outward normal n. Figure  1a illustrates a schematic of the two-dimensional simulation domain used for the fluidbody interaction problem. The center of the square column is located at the origin of the Cartesian coordinate system. The side length of the square column is denoted as D. The distances to the upstream and the downstream boundaries are 20D and 40D, respectively. The distance between the side walls is 40D, which corresponds to a blockage of 2.5%. The flow velocity U ∞ is set to unity at the inlet and a no-slip wall is implemented at the surface of the square column. While the top and bottom boundaries are defined as slip walls, the computational domain is assumed to be periodic in the spanwise direction for the 3D simulations.

Mesh convergence study
For the high-dimensional approximation of full-order model, the computational domain is discretized using an unstructured finite-element mesh, wherein a boundary layer mesh surrounding the body and three-node triangle (2D) and six-node wedge (3D) elements outside the boundary layer region. Three more grids are generated where the mesh elements are successively increased by approximately a factor of 2, designated as M2, M3 and M4. The discretized domain, along with a close-up view of the corners of the square column is illustrated in figure 1b. Results of grid convergence study are recorded in Table (1) for the lock-in region. All cases for the mesh convergence are simulated at Re = 100, m * = 3 and U r = 5.0. The mesh convergence error is computed by considering the finest mesh M4 as the reference case. The force coefficients, the shedding frequency and the root mean square (rms) of the transverse amplitude are analyzed. It can be seen that values recorded for mesh M3 and M4 differ by less than 1%. Therefore, the mesh M3 is adequate for the present study. Furthermore, the adopted full-order solver and the numerical discretizations have been extensively validated in several earlier studies for both low Re (Miyanawala et al. 2016;Jaiman et al. 2016a) and moderate Re (Jaiman et al. 2016a;Miyanawala & Jaiman 2018) flows. In the next section, the modal decomposition of the pressure field is presented for a representative reduced velocity of U r = 6.0 in the lock-in region at (Re, m * , ζ) = (100, 3.0, 0). The snapshots of the FOM performed for the flow past a vibrating square cylinder is utilized to recover the POD modes and the DEIM points. The accuracy of the linear POD and POD-DEIM are systematically assessed with regard to their effectiveness to extract the flow features.
4. Assessment of low-order model for wake decomposition As described earlier, we incorporate the snapshot POD method described to obtain the low-dimensional decomposition of the wake dynamics. As found in Miyanawala & Jaiman (2018), the laminar bluff body flow involves simply a few significant features. It will be ineffective to generate the entire set of POD modes, e.g., the order of the mesh points of 87,120 for this particular problem. Hence, we use the snapshot POD technique and obtain just the most significant POD modes, which are a few order of magnitude smaller. We reconstruct the pressure field using the linear and nonlinear techniques and compare their effectiveness to capture the organized wake features. In the present analysis, the unsteady pressure field values for all the mesh points, are collected to a m × k matrix P where m (mesh count) = 87, 120 and k (number of snapshots) = 320. The fluctuation matrixỹ m×k is then generated by subtracting the mean value (P) of each point over the snapshotsỸ = P − P. The POD modes are extracted using the eigenvalues Λ k×k = diag[λ 1 , λ 2 , ..., λ k ] and eigenvectors W = [w 1 w 2 ... w k ] of the covariance matrixỸ TỸ ∈ R k×k given bỹ Y TỸ W = ΛW. As presented earlier, the POD modes V = [v 1 v 2 ... v k ] are related to Λ and W by V =ỸWΛ −1/2 . Each eigenvalue represents the energy/strength of the POD mode. Since the mean pressure distribution is initially removed from the pressure field, the relative strength of the mode directly expresses the contribution from each mode for the pressure fluctuations. Figure 2a displays the energy of these modes normalized by the total energy of the 320 modes obtained. It is clear that this energy decays exponentially and the most energetic mode has 56% of the total energy. In fact, the first 9 most significant modes contain 99% of the total energy of the modes, as shown in figure 2b. Initially, these 9 significant modes are used to recover the pressure field in the linear POD reconstruction. We refer to these modes as mode 1, mode 2, etc. and they are in the descending order of mode energy (λ i ). We first incorporate the linear reconstruction method wherein we assume the final flow field is a linear combination of the flow features captured by the POD modes.

Linear POD reconstruction
In the linear POD reconstruction method, the instantaneous pressure field is recovered by the mean and a linear combination of the identified significant modes. In this analysis, r is set to 9, which represents the most energetic modes containing ∼ 99% of the total contribution to the pressure fluctuations. The temporal coefficientsŷ j are determined by the L 2 inner product between the fluctuation matrix and the modes as expressed in Eq. (2.14). The mean pressure distribution and the first 9 POD modes are displayed in figure 3. The mean field is symmetric around the X-axis along the wake centerline. This is expected as the time-averaged distribution of the flow past a symmetrical bluff body should be symmetrical. Furthermore, the modes 2, 4, 5, 6, 7 and 8 are symmetric around the wake centerline while modes 1, 3 and 9 are anti-symmetric with equal values and opposite signs about the wake centerline. It is evident that the first, third and ninth modes correspond to the Karman vortex street with alternating positive and negative pressure regions about the X-axis and the pressure contours resulting from a staggered vortex street. The POD modes 2 and 8 have a high gradient behavior in the near-wake (0.5D − 5D) region almost parallel to the top and bottom edges of the square cylinder suggesting that this mode represents the influence from the shear layer. The modes 4, 5, 6 and 7 originate from the near-wake region and diffuse symmetrically towards the far wake. We can attribute these contributions to the near-wake bubble and its local dynamical property. For the ease of explanation, we refer to these modes as the vortex shedding, the shear layer and the near-wake. Figure 4 quantifies the time-invariant contributions (F i j ) from each mode to the drag and lift forces. For the definition of F i j , j = (x, y) is the direction of the force and i is the mode number. These values are calculated based on the fluid-solid boundary values of the mode fields displayed in figure 3. It is clear that the vortex shedding modes (modes 1, 3 and 9) contribute entirely to the lift force, while the shear layer and near-wake modes contribute entirely to the drag force. Further details on the force decomposition procedure using the modal contributions are presented in Appendix A. Due to this directionally independent contribution of the bluff body features for the forces, the time coefficients (ŷ j (t)) of these modes should display the same frequencies of the lift and drag forces.
The time histories and the FFT spectra of the first 5 POD modes are shown in figure  5a. The time histories have sharp variations, because of the snapshot sampling of the data. The first and third mode coefficients have a low-frequency sinusoidal variation with the natural frequency (f n ). The second mode coefficient has a non-zero mean with a ≈ 2 × f n frequency. The fourth and fifth modes have a nearly zero mean variation with ≈ 2 × f n frequency. Interestingly, as presented in figure 5b, f n and 2f n coincide with the frequencies of lift and drag, respectively. Hence, we can further confirm that the modes 1 and 3 make their sole contribution to the fluctuating lift while the modes 2, 4 and contribute to the drag force. Furthermore, the non-zero mean drag can be attributed to the non-zero mean of mode 2, which contains the largest component of the drag. From these observations, we can further confirm that the vortex shedding process contributes exclusively to the lift force and the near-wake and the shear layer phenomena influence the drag force.
Using the POD procedure, we successfully decompose the flow field into physically significant features. We reconstruct the same field combining these modes in the linear POD technique: such that, utilizing Eq. (2.13) for the pressure field, gives P (t) ≈ P + r j=1ŷ j (t)v j . Figure 6 illustrates the pressure distribution at tU ∞ /D = 100 using the linear POD reconstruction. The recovered POD mode is compared with the result obtained from the full-order model. A good match with a maximum relative local error < 2% can be seen in figure 6. To quantify the accuracy of the entire flow field recovery, the normalized root mean square (rms) error of the entire distribution is considered. The Mode Number (i) rms error rms is given by: where P F OM and P P OD are the pressure values of the mesh nodes extracted from the fullorder model and the POD reconstruction, respectively and |P 100 | is the mean pressure of the field. When 9 modes are used, this error is rms = 3.78%. In this linear reconstruction, the highest error is observed at the regions known to exhibit a nonlinear variation, such as the near-wake region, the shear layer and the vortex cores. Next, we analyze the POD-DEIM technique to improve the accuracy in these nonlinear flow features using the snapshot sequence and their respective DEIM points.

Nonlinear POD-DEIM reconstruction
The linear POD reconstruction has the highest error in the nonlinear regions. To reduce this error, more POD modes should be added to the reconstruction which makes the POD-based reconstruction very expensive. Instead, when the DEIM technique is used, it reduces the calculation load while properly capturing the nonlinearity of the field variable. The DEIM utilizes two POD bases using the snapshot method, namely a first POD basis V from the snapshot sequence, and a second basis U from the nonlinear snapshots via the DEIM points. However, unlike the linear POD reconstruction, the accuracy does not necessarily improve with the number of DEIM points and the number of POD modes employed. Using many DEIM points result in adding contributions from some non-significant indices. In figure 7a, it is clearly seen that for 100 DEIM points there are few mesh points which lie away from the significant nonlinear region are taken into the calculation. Further in Table (2), we quantify the number of points in the nonlinearwake region as we increase the number of DEIM points. It is clear that the percentage of points in the critical region decrease as we include more points for the DEIM calculation. Owing to the nonlinear combination of the POD modes, including additional insignificant modes can increase the total error. As shown in figure 7b, the lowest rms = 4.05% can be obtained when 70 DEIM points are used with 7 POD modes. It further establishes that the linear POD reconstruction is generally accurate in a global sense, i.e. the entire flow field reconstruction, in contrast to the nonlinear POD-DEIM. However, figure 8 demonstrates the reconstructed pressure distribution at tU/D = 100 using the optimum number of points and POD modes. There is a significant reduction in the local error as it allows to capture the nonlinear regions more accurately.
Apart from the global and local accuracy, we further assess the computational time consumed by the linear POD and nonlinear POD-DEIM reconstructions. The detailed analysis is presented in Appendix B. Theoretically, the DEIM reconstruction process should be ≈ 64 times faster than linear POD reconstruction and the total DEIM process should be ≈ 3.14 times faster. In the actual computations, when just the reconstructions are considered, the DEIM is 9.28 times faster than the linear POD. When the total processes are compared, DEIM has a speedup of 3.98. In terms of accuracy, DEIM is more accurate in a local sense since it captures the nonlinearities better than the linear POD reconstruction. However, when the entire fluid domain is considered, the linear POD method is more accurate than the DEIM. It is likely that the DEIM introduces unnecessary nonlinearities to the potential regions slightly changing the reconstructed field values. When decomposing and reconstructing the laminar flow fields, both linear and nonlinear methods perform to a satisfactory level. Both methods are capable of reaching the required threshold in a reasonable computational time while accurately capturing the flow features of the wake. Here onwards, we employ the POD-DEIM since it has an improved accuracy when capturing the nonlinearities in the flow field at a lower computational cost.

Drag and lift modes
In this section, we analyze the behavior of different modes in the near cylinder region and explain the exclusive nature of their contributions to the drag and lift forces exerting on the oscillating cylinder in a uniform flow. Herein, the vortex shedding modes are referred to as the lift modes, while the shear layer and the near wake represent the drag modes. Figure 9 displays the combined variation of the lift modes during a single cycle of lift. Note that the motion of the cylinder is not shown for this reconstruction as the POD modes are time invariant. The lift modes vary in an alternating manner in four quadrant. The variation is anti-symmetric about the streamwise centerline. In the maximum lift case (Point B in Figure 9a), the positive pressure force difference (i.e. +Y direction) in the downstream quadrants dominates the small negative difference in the upstream quadrants and vice versa for the minimum lift (Point D). In the zero lift cases (Point A and C), the upstream and downstream pressure force differences tend to become equally strong and they cancel each other. The lift modes vary in such a way that the force on the top 2 quadrants is equal in magnitude and opposite in direction to the force on the bottom 2 quadrants. Due to this force cancellation, the vortex shedding (lift) modes have no contribution to the drag force. Hence, the FFT of the drag force does not contain the corresponding harmonic of the natural frequency (f n ). Figure 10 describes the variation of drag modes with the fluctuation of the drag force. The drag fluctuation is defined as C D = C D (t)−C D . The drag modes vary symmetrically around the wake centerline hence offer no contribution to the lift. Similar to the lift modes, this explains the absence of a 2f n harmonic in the lift. The maximum drag fluctuation (Point C) is higher than the minimum drag fluctuation (Point A). Further, at the zero drag fluctuation points (B and D), the magnitude of the drag fluctuation remains a positive value. This further confirms that the drag modes exert a non-zero mean drag on the bluff body apart from the drag force due to the base flow.
With the aforementioned observations, the decomposition of fluid force on a moving bluff body based on the contributions from different POD modes can be expressed as where F j is the force in a particular direction (j = x for in-line and j = y for transverse). F 0 j is the time independent contribution from the mean field and F i j is the time independent pressure fluctuation contribution calculated for the i th mode. While Shear layer and near-wake modes is the time dependent coefficient of the i th mode for the force in direction j, n rj is the number of POD modes with a significant contribution to the particular force. Using the snapshot data, we can determine b i j (t) for the streamwise and transverse forces as The complete description of this force decomposition and its usage is provided in Appendix A. Table (3) summarizes the qualitative analysis of the contributions from the mean field and the modes to the drag and lift forces. The mean-field has a symmetric pressure distribution about the wake centerline, hence contributes solely to the time-independent component of the drag force. The vortex shedding modes have an anti-symmetric pressure distribution throughout the time history, hence they have no drag force contribution. We observe that these lift force contributions have a near zero mean (similar to the lift variation) as well. The shear layer and near-wake modes have the same qualitative properties of the mean field, however their contributions to the drag force are timedependent. The POD modes provide a deep insight into important flow features and their contribution to the wake dynamics. It is important to investigate their variation with different flow conditions, i.e. the parameters mentioned in Eq. (3.1). The variation of POD modes and their contribution to wake dynamics with the reduced velocity is examined in the next section, with the goal to explain the role of the wake features in sustaining the synchronized wake-body motion.

Wake feature interaction and sustenance of VIV lock-in
In this section, we investigate the relative contributions from different features to the pressure fluctuations and eventually the fluid forces on the freely oscillating bluff body. When a bluff body is free to oscillate in a current flow it undergoes the lock-in phenomenon: the oscillation amplitude significantly increase when the natural frequency of the bluff body approaches the vortex shedding frequency. In Miyanawala et al. (2016), the lock-in phenomenon for a square cylinder immersed in a laminar flow at Re = 100 is systematically studied. Figures 11a and 11b summarize the bluff body dynamics of a freely vibrating square cylinder. The cylinder undergoes wake-body synchronized lock-in in the range U r ∈ [4.5, 7] and the peak oscillation occurs at U r = 5.0. Figures 11c-f elucidate the variation of relative contributions from different modes as a function of the reduced velocity (U r ). It is interesting to note that the three most energetic modes correspond to the same flow features throughout the U r range namely the first and third modes (vortex shedding) and the second mode (shear layer). However, the modes 4-10 vary in this regard, where most of these modes correspond to the near-wake phenomena. We quantify the relative energy contribution from each wake feature by summing the mode energy of the corresponding modes, i.e.
where E j is the relative energy contribution from the wake feature (j = vs for the vortex shedding, j = sl for the shear layer and j = nw for the near-wake), n j is the number of significant modes corresponding to a particular flow feature and k is the total number of modes. As shown in figure 11c, the total contribution from the vortex shedding increases in the lock-in region. However, the first mode becomes more energetic while the third mode is relatively less energetic in this region. Figure 11d exhibits that the contribution from the shear layer modes reduces significantly in the lock-in region. All the individual shear layer modes also follow a similar trend. The near-wake modes depicted in figure 11e become more energetic in the lock-in region relative to the pre-and post-lock-in regions.
Unlike the vortex shedding and shear layer, the primary and secondary near-wake modes have remarkably similar contributions. A summary of the contributions from the 10 most energetic POD modes corresponding to different physical phenomena is illustrated in figure 11f. Note that these 10 modes capture ≈ 99% of the total mode energy. It is clear that the shear layer contributions decrease while the vortex shedding and the nearwake contributions increase in the lock-in region. In the post-lock-in region, the relative contributions from the flow features remain almost constant. However, our intuition poses a question: it is evident that the extreme transverse motion in the lock-in is sustained by the resonance of the bluff body spring-mass system and the vortex shedding since the transverse motion frequency matches the vortex shedding frequency. Why does not the energy transfer from the fluid to the bluff body make the vortex shedding modes less energetic instead of what we observe here? In that relation, we propose a cycle to explain this counter-intuitive behavior of the decomposed wake features. Figure 12a elaborates the interaction between the wake features and the bluff body motion. When the vortex shedding synchronizes with the bluff body motion, it causes the bluff body to undergo a relatively high-amplitude motion. This widens the wake and eventually the shear layer, decreasing the velocity gradients. This causes the shear layer to give away vorticity flux to the vortex shedding process, intensifying the vortices and the near-wake bubble. The strengthening of vortices increases the in-phase forces with the motion, i.e., the surrounding fluid flow tends to supply higher energy to the structure. As illustrated in Jauvtis & Williamson (2004), the force F v and the energy transfer ratė e v due to the principal vortices can be analyzed using the following simple analytical relations: where Γ is the vortex strength, U v is the streamwise velocity of the predominant vortex relative to the bluff body andẏ is the transverse velocity of the bluff body. It is clear that the increase in the vortex strength will increase the forces and energy transfer to the bluff body. The widening of the high gradient shear layer region in the lock-in regime can be seen in figures 12b-e, which demonstrate the primary shear layer mode for the different U r cases. In the pre-lock-in regime, the near-wake region is positive compared to the shear layer region. When U r = 5.0, the maximum amplitude case, it is clear that the high gradient region has shrunk in the streamwise direction and expanded in the transverse direction. Consequently, the near-wake region and the shear layer region interchanges the distribution when U r = 6.0, i.e. the near-wake region is negative compared to the shear region. This sign change continues to the post-lock-in regime, where the high gradient region extends to the streamwise direction and becomes narrower in the transverse direction. We further generalize this variation of the wake feature contribution for Re > Re cr (= 44.7 (Yao & Jaiman 2017)). Figure 13 demonstrates the bluff body motion response and the modal energy contribution from the large-scale features of the wake. The cylinder motion follows a similar trend for Re = 100, 125 and 150 where the lock-in region is detected as U r ∈ [4.5, 7]. This region is slightly shifted to U r ∈ [5.5, 8] for Re = 70. In all cases, we observe a maximum of A rms y ≈ 0.2D. Regardless of Re, the wake features exhibit a similar trend in terms of modal energy. As displayed in figures 13(b) and (d), the vortex shedding and the near-wake modes become more energetic during the lock-in and the shear layer modes become less energetic. This further confirms the proposed interaction cycle for the coupling of the wake features and the bluff body motion.
Using the modal decomposition, we have quantitatively explained the interaction dynamics of the flow features which have been conjectured by many previous studies. For example, many successful VIV suppression techniques are proposed by passive (Law & Jaiman 2017) and active (Guan et al. 2017;Narendran et al. 2018) methods with the experience based presumption that preventing the interaction between the shear layer, the vortex street and the near-wake will suppress the synchronized wake-body lock-in  phenomena. The cycle proposed above provides a proper physical mechanism for the success of those methods: they prevent the vorticity transfer between the shear layer and the vortex shedding and/or near-wake bubble, which breaks the self-sustenance of the wake interaction cycle. This understanding of the wake features and their interactions will be vital for the development of effective suppression methods and devices for flowinduced vibrations. In this analysis, we observe that the synchronization of the wake and bluff body weakens the shear layer and intensifies the vortices and the near-wake bubble. In the aforementioned analysis, there exists a periodic vortex shedding for the stationary and pre-/post-lock-in regimes. In what follows, we investigate whether the perturbation of the near-wake bubble via flexibility can sustain the synchronized wake-body interaction at below critical Re flow wherein the well-defined periodic vortex shedding does not exist for the stationary counterpart.

Synchronized wake-body interaction at below critical Re
At very low Re, the flow past a bluff body is two-dimensional, steady and symmetric with respect to the wake centerline. The near-wake bubble attached to its surface is the essential feature below the critical Reynolds number Re cr , which is formed by the steady separation from the sharp corners of a square cylinder. Two symmetric and counterrotating recirculation zones are present in the wake bubble. As Re increases above the critical value, a Hopf bifurcation sets in and the flow becomes periodic via vortex shedding process. For circular and square cylinders Park & Yang (2016) demonstrated these values to be Re cr = 46.8 and 44.7 respectively, which were further confirmed by Yao & Jaiman (2017). Interestingly, when the bluff body is free to vibrate, Meliga & Chomaz (2011) predicted for a circular cylinder that this unstable boundary will hold when m * > 1000 and the Hopf bifurcation will occur at much lower Re for low m * values. The authors further conjectured that for a circular cylinder there will be a limiting Re ∼ 32, below which the wake flow will be 2D and steady regardless of the mass ratio m * .
Herein, we observe that for some Re < Re cr the spring-mounted square cylinder undergoes significant synchronized wake-body motion for a specific range of U r . Figure  14a illustrates a variation of high amplitude motion for Re = 30, 35 and 40. When Re becomes closer to Re cr , the synchronization regime widens and the highest amplitude U r shifts from U r = 8 to U r = 7. In contrast to Re > Re cr cases, we observe no motion of the cylinder in pre-and post-synchronization regimes. We further examine the conjecture of Meliga & Chomaz (2011) and demonstrate that for a square cylinder this synchronized motion is present when Re 26 (figure 14b). Additional analysis on these synchronized motion cases revealed that the wake-body system synchronizes to a frequency slightly less than the natural frequency f n of the bluff body similar to the Re > Re cr lock-in regime (figure 14c). We observe that, for all synchronized motion cases, the wake is unsteady with some vortex shedding patterns. For example, we can compare the representative Z-vorticity contours for the zero motion and the synchronized motion cases displayed in figure 14d and (14e), respectively. The zero motion case is almost identical to the stationary cylinder counterpart, while the synchronized motion case is similar to the lockin scenario for Re > Re cr . Notwithstanding, the vortex formation length is considerably high for this below Re cr configuration. We further decompose the unsteady wake of the synchronized motion case and examine similar features as Re > Re cr cases, i.e. the vortex shedding ( figure 14f and h), the shear layer (figure 14g) and the near-wake bubble (figure 14i).
These observations constitute the basic requirement for the wake-bluff body synchronized motion: the bluff body should have an optimal amount of flexibility (i.e., not too rigid nor too flexible) and the flow needs to have sufficiently large inertia (i.e., higher Re) to trigger the unsteadiness in the near-wake bubble. This particular Re is lower than the Re cr for a fixed bluff body. This means that the flexibility of the solid body provides an avenue for the wake and the spring-mounted body to synchronize eventually causing the wake to be unsteady. From this numerical experiment, we can deduce that the flexibility of the bluff body is the primary factor driving the synchronized wakebody motion, neither the vortex shedding nor the shear layer roll-up. Hence the most critical wake feature for the onset of wake-body synchronization is the near-wake bubble. When the Re is very low (< 26) this bubble remains steady and the counter-rotating recirculation zones behind the bluff body are stable. For 26 Re 44.7, it remains same if the bluff body is either too rigid or flexible. However, in this Re regime, when the bluff body is appropriately flexible, slight perturbations cause distortions in the steady wake. These distortions become periodic and begin to synchronize with the bluff body. This synchronization leads to a relatively higher amplitude oscillations at a frequency slightly less than the natural frequency of the solid body in a vacuum. At the same time, due to the vorticity generation by the unsteadiness developed in the wake, the vortices are shed from the downstream end of the wake bubble.
With the aforementioned findings, we summarize the bluff body wake behavior of 20 Re 50 regime in figure 15. We deduce that the root cause of wake-body synchronization is the frequency lock-in between the natural frequency of the bluff body and the nearwake bubble. This further demonstrates that the unsteadiness of the wake can be induced by the flexibility of the bluff body. Moreover, the unsteady wake alone cannot induce high amplitude bluff body oscillations (e.g., pre-and post-lock-in in Re > Re cr ). Hence, we can further infer that the wake-body synchronization is induced by the synchronization of the bluff body motion with the near-wake bubble, not with the vortex street. In the next section, we generalize our findings to three-dimensional turbulent flows at moderately high Re.

Effect of turbulence
In this section, we investigate the dynamic decomposition of the wake behind a threedimensional oscillating square cylinder at Re = 22, 000, wherein the wake is fully turbulent. Our aim is to understand the role of turbulence when we extend the wake feature interaction cycle to turbulent flow. To retrieve the high-fidelity data at this Re, we employ a well-established dynamic subgrid-scale turbulence model in our finite-element formulation. The filtered Navier-Stokes formulation and the determination of the subgrid stress term via the dynamic subgrid-scale model is provided in Jaiman et al. (2016a). We incorporate this full-order model to generate 3D snapshots of the flow fields and the POD-Galerkin projection is applied on this high-fidelity data set. At high-Re turbulent wake flow, the aforementioned large-scale organized flow features are fragmented into smaller scales until the scales are fine enough to dissipate by the fluid viscosity. Therefore, smallscale and high-frequency modes can have a significant impact on the overall dynamics for the high-Re turbulent condition, in contrast to the low-Re study. Figure 16 demonstrates the mode energy distribution for Re = 22, 000. Compared to the low cases, the modal energy is much more distributed among the modes. For instance, the most energetic mode of the Re = 100 case contains 56% of the total energy while it is 32.88% for the high-Re case at Re = 22, 000. Due to this broadening of the mode energy, the energy decay is less steep. For the low Re cases, the first 5 modes contain 95% of the total mode energy and the first 9 modes contain 99%. For the high-Re case, a total of 123 modes is required to capture 95% of the energy and 211 modes are required for 99%. We further investigate this distribution of the mode energy with the presumption that the presence of broadband turbulence is the key factor. Figure 17 displays the mean pressure field and 10 most energetic POD modes obtained using the same snapshot technique for moderate Re case. All the modes exhibit distorted and scattered patterns compared to the low Re cases. However, the POD modes further portray general large-scale features. For example, the modes 1, 3, 4, 5, 6, 7 and 9 exhibit the flow patterns related to vortex shedding while the mode 2 is related to the shear layer and the modes 8 and 10 are of the near-wake phenomena. The distortions occur in each mode throughout the spatial domain due to the broadband and multiscale phenomenon of turbulence, which are not decomposed by the singular value decomposition. Turbulence distributes the modal energy across the modes which make it require significantly more modes to reconstruct the flow field and the underlying wake dynamics. Hence, the POD based reconstruction becomes computationally more expensive in turbulent flows due to the broadband and multiscale character.
Similarly, figure 18 demonstrates the broadband nature of turbulence in the temporal domain. Even with the multiscale spatial distortions, the first mode of the moderate Re has a similar temporal contribution as the first mode of the low Re case with a single dominant frequency close to the natural frequency of the system. However, the temporal coefficients of the other modes have multiple harmonics. Some of the modes exhibit predominant frequencies among the broadband FFTs. For example, the mode 2 has a dominant 2f n frequency behavior, the mode 3 has f n , 2f n and 3f n harmonics and mode 4 has f n and 2f n harmonics. These multiple harmonics occur due to the bombardment of turbulence on the corresponding flow features of the POD modes. Figure 19 displays the pressure field reconstruction using the POD-DEIM technique. The actual instantaneous field contains some distortions and fine near-wake variations which are not completely captured by the reconstruction. However, the general large-scale variations are properly reconstructed with a maximum local error of ≈ 2%. The broad- band energy distributing nature of turbulence has reduced the contribution of significant POD modes to the wake dynamics. Due to this behavior, many flow features are not captured when few of the most energetic modes are considered for the reconstruction. Hence, the inclusion of many POD modes is required for an accurate reconstruction which makes the POD reconstruction computationally expensive and time-consuming. However, this can be mitigated by selective reconstruction of few required timesteps instead of the entire time history. Using this to our advantage, we investigate the validity of the wake-body interaction cycle at this moderate Re. Figure 20 illustrates the response characteristics and the wake feature contributions for the mode energy at Re = 22, 000. In contrast to the low Re cases, the bluff body undergoes galloping at this Reynolds number. We observe the same wake-body synchronization reported in Miyanawala & Jaiman (2018), i.e. 1:1 frequency synchronization of the bluff body motion and force at U r ∈ [5, 6] and 1:3 synchronization at U r = 10.0. Our analysis here is focused on the 1:1 frequency lock-in where the vertical fluid forcing and the bluff body motion are in synchronization. According to figure 20(b), it is clear that the vortex shedding and near-wake modes become more energetic in the synchronized regime while the shear layer modes become less energetic. This is the same behavior observed in the low Re analysis which proves that the proposed wake-body interaction cycle is valid even with the presence of turbulence.
In summary, the presence of turbulence distorts the spatially symmetric/antisymmetric nature of POD modes and distributes the mode energy throughout many POD modes. Hence the reconstruction requires many modes and the classification of features is more complex at moderately high Re. However, even with this complexity, the fundamental wake-body interaction process proposed using low Re analysis is observed to be valid for the three-dimensional turbulent flows. Hence we can conclude that the wake-body interaction cycle proposed in this study is a general cycle for fluid-structure interaction systems.

Concluding remarks
Despite the prevalence of SVD-based modal reduction techniques, there are very few studies on their application to fluid-structure interaction systems. In this paper, we considered the 2-DOF free vibration of a square cylinder under laminar and turbulent flows. We explored the capability of POD decomposition to interpret the most significant wake features and their contributions to the forces on the vibrating body interacting with fluid flow. When the linear and nonlinear POD-DEIM reconstructions are contrasted, we found that the DEIM method is faster and has a higher local accuracy since it captures the nonlinearity of the principle vortices and the near-wake region. For the low Re cases, every POD mode clearly represents one of the large-scale flow features: vortex shedding, shear layer or near-wake bubble. In these cases, we further observed that the nine most energetic modes contain ≈ 99% of the energy. Further, we identified that the vortex shedding modes solely contribute to the transverse (lift) force while the shear layer and the near-wake modes solely contribute to the drag force. Based on these observations, we proposed a novel force decomposition for the drag and lift forces which is different from the conventional force decomposition based on the added mass and the viscous force contributions.
We examined the POD decomposition for a range of U r values and we proposed the mechanism of the sustenance of synchronized wake-body lock-in. This further provided the explanation to the counter-intuitive observation: in the lock-in region, even though the kinetic energy is transferred from the fluid to the bluff body, the principle vortices are much more energetic than the pre-and post-lock-in regimes. It is seen that the bluff body motion widens the wake and it causes a vorticity transfer from the shear layer to the near-wake and vortices. We proposed the wake feature interaction cycle based on these observations. We further confirmed that this mechanism is valid for laminar Re > Re cr range. For below critical Re flows, we observed that the bluff body and the wake still synchronize and undergo large amplitude motion at some U r values when Re 26.
Decomposition of these wakes further exhibited a similar behavior as the synchronized large amplitude motion cases at Re > Re cr . This revealed that the flexibility of the bluff body induced the unsteadiness in the near-prwake bubble causing it to break and generate the vortices. With this observation, we can conclude that the fundamental requirements for the wake-body synchronized motion are, large enough flow inertia and appropriate flexibility of the structure. When the moderate Re turbulent bluff body flow is decomposed, we observe that all the dominant wake modes are bombarded with different scales of turbulence. The broadband nature of turbulence resulted in a wide mode energy distribution, which required up to 123 modes to reach the 95% mode energy threshold. Further analysis of the time coefficients of the modes confirmed that the large-scale features are battered by the multiple frequency turbulence. However, they generally correspond to a large-scale wake feature similar to the laminar cases. The wake decomposition of turbulent flows for U r ∈ [3, 10] confirmed that the wake interaction cycle proposed for laminar cases is valid for turbulent flows as well.

Appendix A. Force decomposition based on modal contribution
The aim of this appendix is to present a general decomposition of the fluid force exerted on a moving body in an incompressible viscous flow. To begin, we provide some background on existing force decomposition techniques that characterize the fluid inertial and the viscous forces on a moving body. In one of the pioneering works, Morison et al. (1950) proposed a force decomposition for the in-line force acting on a cylindrical object which is widely used in many engineering applications. This semi-empirical decomposition can be written as a linear sum of a velocity squared-dependent drag force and an acceleration-dependent inertial force: where C d and C m represent the averaged drag and inertia coefficients, which can be determined by experiments or numerical computations. Owing to the nonlinear dependency of these coefficients on the evolution of vorticity field, Sarpkaya (2001) argued that "It does not perform uniformly well in all ranges of K, β and k/D", where K denotes the Keulegan Carpenter number, β = Re/K and k/D is the relative roughness. In Lighthill (1986), a different approach is taken by the assertion that the viscous drag and the inviscid inertia force operate independently, by re-writing Eq. (A.1): and for a flow defined by U (t) = −U m cos ωt it reduces to where A p and V b denote the projected area and the volume of the body, respectively and C * m is the ideal value of the inertia coefficient. Notwithstanding, many studies demonstrate that it is difficult to represent the actual force with this relation as long as a constant value C * m is considered. In particular, Sarpkaya (2001) clearly demonstrated that the viscous drag force and the inviscid inertia force are not completely independent and it is impossible to decompose the unsteady drag force to an inviscid and a vorticity-drag component. The decomposition of the total force into inviscid and viscous components by Lighthill's relation (Eq. A.2) can be considered as an effort to lump the effects of the complex generation and evolution of vorticity field into mutually independent forces related to the inviscid inertia and the viscous effects. In such force decomposition techniques, the characteristics vorticity patterns and their dynamics generated during the motion of a body are not included.
In what follows, we propose an alternative force decomposition for the in-line (drag) and transverse (lift) pressure forces applied on a bluff body which extends the above decompositions to incorporate significant features of unsteady separated flow. In particular, the unsteady force is decomposed to include the nonlinear generation and evolution of vorticity field around a moving body in a fluid flow. This decomposition is based on the contributions from different POD modes to the forces and can be written as where F j is the force on a particular direction (j = x for the in-line and j = y for transverse). While F 0 j is the time-independent contribution from the mean field, F i j is the unsteady pressure fluctuation contribution associated with i th mode. Here, b i j (t) is the time-dependent coefficient of the i th mode for the force in direction j and n rj is the number of POD modes with a significant contribution for the particular force. Using the snapshot data, we can exactly determine b i j (t) for the in-line and transverse forces by the following relation The magnitude of the modes, the time coefficients and the force contributions introduced in this decomposition slightly fluctuate when flow parameters and the bluff body geometry are changed. Similar to the above methods, we can create databases of F i j for different bluff bodies. These databases can then be used to determine the total forces as well as the contribution from each flow feature to the bluff body dynamics. To further generalize the force decomposition, deep learning techniques  for parametric predictions can be employed. Figure A.1b and A.1c present the reconstructed values of the time dependent coefficient for drag and lift modes. Note that the relevant six out of the first 9 modes are considered for the drag (i.e. n rx = 6 and n ry = 3) based on the values of figure A.1a. In a nutshell, the force component represented by the modal decomposition implicitly characterizes the three constituent components involving an inviscid inertial force, the dynamics of vorticity field, and a skin friction force.

Appendix B. Performance comparison of POD reconstruction methods
Herein, we briefly compare the number of floating point operations (FLOPs) required for the linear and DEIM based POD reconstructions in Table (B.1). The generation of the POD modes which is essential for both reconstructions requires ∼ O(4k 2 m) FLOPs. In this study, we estimate this value to be 3.571×10 10 . The linear POD reconstruction needs ∼ O(rk 2 m) computations where the equal main contributions are from the multiplication between the time coefficient and the POD modes, and the summation of the multiplied POD contributions. With the use of 9 POD modes, the estimated FLOP count is 8.029 × 10 10 . The DEIM technique consumes fewer computations than POD as it requires ∼ O(2n 2 p m + 2mkn d ) FLOPs. Generating the matrix solution in DEIM first step is the most expensive operation as it needs more FLOPs for last DEIM points. With the use of 7 POD modes and 70 DEIM points, the POD-DEIM reconstruction needs 1.244 × 10 9 FLOPs. According to this calculation, the linear POD reconstruction demands ≈ 64.54 times more FLOPs than the DEIM reconstruction. However, the overall linear POD process needs ≈ 3.14 times more FLOPs than the overall DEIM process.
The actual performance of the linear POD and POD-DEIM reconstruction techniques is compared in Table (B.2). All the calculations are performed using the same computer. When just the reconstructions are considered, DEIM is 9.28 times faster than linear POD. When the total processes are compared, DEIM has a speed gain of 3.98. DEIM is more accurate in the nonlinear flow regions. However, linear POD method is more accurate when the cumulative error of the fluid domain is considered. Table B.1: Number of floating points operations (FLOPs) required for linear and nonlinear POD reconstruction. Mesh count (m) = 87,120, number of snapshots (k) = 320, number of significant modes for linear POD (r) = 9, number of significant modes for DEIM (n d ) = 7 and number of DEIM points (n p ) = 70.

Operation
FLOPs Order