2.1 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:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{f}\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{f}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}^{f}(\boldsymbol{u}^{f}-\boldsymbol{w})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{f}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{f}+\boldsymbol{b}^{f}\quad \text{on }\unicode[STIX]{x1D6FA}^{f}(t), & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{f}=0\quad \text{on }\unicode[STIX]{x1D6FA}^{f}(t), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{s}}{\unicode[STIX]{x2202}t}+\boldsymbol{C}\boldsymbol{u}^{s}+\boldsymbol{K}(\unicode[STIX]{x1D753}^{s}(\boldsymbol{z}_{0},t)-\boldsymbol{z}_{0})=\boldsymbol{F}^{s}\quad \text{on }\unicode[STIX]{x1D6FA}^{s}, & \displaystyle\end{eqnarray}$$

where superscripts $f$ and $s$ denotes the fluid and structural domains, and $\unicode[STIX]{x1D6FA}^{f}(t)$ and $\unicode[STIX]{x1D6FA}^{s}$ represent the fluid and solid domains, respectively. Here $\unicode[STIX]{x1D70C}^{f}$ is the fluid density, $\boldsymbol{u}^{f}$ and $\boldsymbol{w}$ are the fluid and mesh velocities at a spatial point $\boldsymbol{x}\in \unicode[STIX]{x1D6FA}^{f}(t)$ , and $\boldsymbol{b}^{f}$ denotes the body force in the fluid domain. For the structural system, $\boldsymbol{M}$ , $\boldsymbol{C}$ and $\boldsymbol{K}$ are the mass, damping and stiffness matrices of the bluff body and $\boldsymbol{F}^{s}$ is the external force acting on the body. The function $\unicode[STIX]{x1D753}^{s}(\boldsymbol{z}_{0},t)$ maps the initial position vector of the centre of mass ( $\boldsymbol{z}_{0}$ ) to its position at time $t$ , and $\unicode[STIX]{x1D748}^{f}$ is the Cauchy stress tensor for a Newtonian fluid given by

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D748}^{f}=-p\boldsymbol{I}+\unicode[STIX]{x1D707}^{f}(\unicode[STIX]{x1D735}\boldsymbol{u}^{f}+(\unicode[STIX]{x1D735}\boldsymbol{u}^{f})^{\text{T}}),\end{eqnarray}$$

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 $\unicode[STIX]{x1D6E4}$ as follows:

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{f}(t)=\boldsymbol{u}^{s}(t), & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{\unicode[STIX]{x1D6E4}(t)}\unicode[STIX]{x1D748}^{f}(\boldsymbol{x},t)\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\unicode[STIX]{x1D6E4}+\boldsymbol{F}^{s}=\mathbf{0}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{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 modelling of the boundary layer and the vorticity generation over a moving body. While (2.1)–(2.3) of the coupled fluid–body system are directly solved for low- $Re$ flows, we consider the well-established dynamic subgrid-scale model for high- $Re$ turbulent flow. The spatially filtered Navier–Stokes and continuity equations are solved in the variational form. Details of the dynamic subgrid-scale model are provided in Jaiman, Guan & Miyanawala (Reference Jaiman, Guan and Miyanawala2016a ).

The weak variational form of (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 FSI (Jaiman, Pillalamarri & Guan Reference Jaiman, Pillalamarri and Guan2016b ). The motion of the 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 the solid structure becomes large. To extract the transient flow characteristics, we solve the Navier–Stokes equations at discrete time steps, which leads to a sequence of linear systems of equations via Newton–Raphson-type iterations. We employ the conjugate gradient, with a diagonal preconditioner for the symmetric matrix arising from the pressure projection and the standard generalized minimal residual 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 FOM for the FSI.

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:

(2.7) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{y}}{\text{d}t}=\boldsymbol{F}(\boldsymbol{y}),\end{eqnarray}$$

where $\boldsymbol{y}$ is the column state vector describing the unknown degrees of freedom and $\boldsymbol{F}$ is a vector-valued function describing the spatially discretized governing equations. In the present fluid–body system, the state vector comprises the fluid velocity and the pressure as $\boldsymbol{y}=\{\boldsymbol{u}^{f},p\}$ and the structural velocity involves the three translational degrees of freedom. For a discretized domain of $m$ elements and $n$ time steps, the full-order simulation outputs a high-fidelity data set $\boldsymbol{y}\in \mathbb{R}^{m\times n\times q}$ , where $q$ is the number of variables in $\boldsymbol{y}$ . This data set is extremely valuable for determining the instantaneous physics of the fluid–body system and to construct a low-order representation that preserves the behaviour of the original system.

2.2 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 $\boldsymbol{F}$ of (2.7) as

(2.8) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{y})=\boldsymbol{f}+\boldsymbol{A}\boldsymbol{y}+\boldsymbol{F}^{\prime }(\boldsymbol{y}),\end{eqnarray}$$

where $\boldsymbol{f}$ denotes a constant column vector with $m$ rows, and $\boldsymbol{A}$ and $\boldsymbol{F}^{\prime }$ are the linear and nonlinear terms. For ease of explanation, consider the solution vector $\boldsymbol{y}(\boldsymbol{x},t)\in \mathbb{R}^{m\times 1}$ comprising a single quantity of interest which has been determined at discretized locations $\boldsymbol{x}$ of the spatial domain and for a particular time $t$ . The matrix operator $\boldsymbol{A}$ is an $m\times m$ matrix which captures the linear dynamics, while $\boldsymbol{F}^{\prime }(\boldsymbol{y})$ is a nonlinear function of $\boldsymbol{y}$ . Using the projection-based model reduction, we can represent the state vector $\boldsymbol{y}$ by an element in a low-rank vector subspace spanned by the column vectors of an $m\times k$ matrix $\boldsymbol{{\mathcal{V}}}=[v_{1}v_{2}\cdots v_{k}]$ , where $k\ll m$ . The state vector $\boldsymbol{y}$ can be approximated by $\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}$ , where $\hat{\boldsymbol{y}}$ is a reduced column vector with $k$ entries. Since the columns of $\boldsymbol{{\mathcal{V}}}$ are orthonormal (i.e.  $\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{{\mathcal{V}}}=\boldsymbol{I}$ ), via the Galerkin projection onto the basis $\boldsymbol{{\mathcal{V}}}$ , we get the following reduced dynamics:

(2.9) $$\begin{eqnarray}\frac{\text{d}\hat{\boldsymbol{y}}}{\text{d}t}=\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{f}+\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{A}\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}+\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}).\end{eqnarray}$$

Next, we have to choose a suitable subspace for the mode decomposition. Using the reduced singular value decomposition (SVD), the above state vector $\boldsymbol{y}$ can be expressed as

(2.10) $$\begin{eqnarray}\boldsymbol{Y}=\boldsymbol{{\mathcal{V}}}\unicode[STIX]{x1D72E}\boldsymbol{{\mathcal{W}}}^{\text{T}}=\mathop{\sum }_{j=1}^{k}\unicode[STIX]{x1D70E}_{j}\boldsymbol{v}_{j}\boldsymbol{w}_{j}^{\text{T}},\end{eqnarray}$$

where the vectors $\boldsymbol{v}_{j}$ are the POD modes of the matrix $\boldsymbol{Y}$ with rank $k$ , $\boldsymbol{{\mathcal{W}}}$ is an orthonormal matrix with $n\times k$ , and $\unicode[STIX]{x1D72E}$ is a $k\times k$ diagonal matrix with diagonal entries $\unicode[STIX]{x1D70E}_{1}\geqslant \unicode[STIX]{x1D70E}_{2}\geqslant \cdots \geqslant \unicode[STIX]{x1D70E}_{k}\geqslant 0$ . For any $r\leqslant k$ , the subspace spanned by $\{\boldsymbol{v}_{1},\ldots ,\boldsymbol{v}_{r}\}$ provides an optimal representation of $\boldsymbol{y}$ in the subspace of dimension $r$ using the SVD process. The total energy contained in each POD mode $\boldsymbol{v}_{j}$ can be computed by the singular value $\unicode[STIX]{x1D70E}_{j}^{2}$ . Note that $\boldsymbol{{\mathcal{V}}}$ and $\boldsymbol{{\mathcal{W}}}$ are the orthonormal eigenvectors of $\boldsymbol{Y}\boldsymbol{Y}^{\text{T}}$ and $\boldsymbol{Y}^{\text{T}}\boldsymbol{Y}$ , respectively.

2.2.1 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 $\boldsymbol{Y}\boldsymbol{Y}^{\text{T}}$ are determined instead of performing the SVD. The algorithm adopted from Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, Mckeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017) is summarized in Algorithm 1.

The standard linear POD incurs almost the same order of cost as the full-order analysis since it is using the $\tilde{\boldsymbol{Y}}\tilde{\boldsymbol{Y}}^{\text{T}}$ matrix which has size $m\times m$ . In a typical time-dependent 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 Reference Sirovich1987) is applied to extract the most significant modes. In the snapshot method, the eigenvalue decomposition is performed on $\tilde{\boldsymbol{Y}}^{\text{T}}\tilde{\boldsymbol{Y}}\in \mathbb{R}^{k\times k}$ , which is significantly smaller than $\tilde{\boldsymbol{Y}}\tilde{\boldsymbol{Y}}^{\text{T}}$ as $k\ll m$ . Let the eigenvalues and eigenvectors of $\tilde{\boldsymbol{Y}}^{\text{T}}\tilde{\boldsymbol{Y}}$ be given by

(2.11) $$\begin{eqnarray}\tilde{\boldsymbol{Y}}^{\text{T}}\tilde{\boldsymbol{Y}}\boldsymbol{{\mathcal{W}}}=\unicode[STIX]{x1D726}\boldsymbol{{\mathcal{W}}},\end{eqnarray}$$

then, using the relationship between the eigenvectors of $\tilde{\boldsymbol{Y}}\tilde{\boldsymbol{Y}}^{\text{T}}$ and $\tilde{\boldsymbol{Y}}^{\text{T}}\tilde{\boldsymbol{Y}}$ , a maximum of $k$ significant POD modes can be extracted by

(2.12) $$\begin{eqnarray}\boldsymbol{{\mathcal{V}}}=\tilde{\boldsymbol{Y}}\boldsymbol{{\mathcal{W}}}\unicode[STIX]{x1D726}^{-1/2}.\end{eqnarray}$$

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 significant modes identified,

(2.13) $$\begin{eqnarray}\boldsymbol{Y}(\boldsymbol{x},t)\approx \overline{\boldsymbol{Y}}(\boldsymbol{x})+\mathop{\sum }_{j=1}^{r}{\hat{y}}_{j}(t)\boldsymbol{v}_{j},\end{eqnarray}$$

where $r$ is the number of significant POD modes. The temporal coefficients of the linear combination are determined by the $L^{2}$ inner product $\langle \hspace{2.22198pt}.\hspace{2.22198pt},\hspace{2.22198pt}.\hspace{2.22198pt}\rangle$ between the fluctuation matrix and the modes as follows:

(2.14) $$\begin{eqnarray}\hat{\boldsymbol{Y}}(t)=\langle \boldsymbol{Y}-\overline{\boldsymbol{Y}},\boldsymbol{{\mathcal{V}}}\rangle .\end{eqnarray}$$

This summarizes the process of POD by performing the SVD on the snapshots of the sampled solutions at certain time steps.

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 that the nonlinearities behave almost linearly. Subsequently, the POD reconstruction may result in a similar order of computational expense to full-order simulation. This issue can be handled by employing the 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.

2.2.2 Discrete empirical interpolation method

To overcome the difficulty in the linear POD, Chaturantabut & Sorensen (Reference Chaturantabut and Sorensen2009) proposed the DEIM to reconstruct the full-order variable as a nonlinear combination of the POD modes. The aim of the DEIM is to design a low-order representation for the nonlinear terms by introducing an additional basis. Consider $\boldsymbol{{\mathcal{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 (2.8) by the sequence of nonlinear snapshots as $\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}(t))\approx \boldsymbol{{\mathcal{U}}}\hat{\boldsymbol{c}}$ . The coefficients $\hat{\boldsymbol{c}}$ can be selected based on Algorithm 2, which relies on a greedy nonlinear function approximation. In Algorithm 2, $\hat{\unicode[STIX]{x1D70C}}$ and ${\wp}_{1}$ denote the assigned value and the assigned index of $\max \{|\boldsymbol{v}_{1}|\}$ , and $\boldsymbol{e}_{{\wp}_{i}}=[0,\ldots ,0,1,0,\ldots ,0]^{\text{T}}\in \mathbb{R}^{m}$ is the ${\wp}_{i}$ th column of the identity matrix of size $m\times m$ . The accuracy of the DEIM approximation depends on the error induced by the POD projection and the estimation of $\Vert (\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{U}}})^{-1}\Vert$ . Further details of the DEIM process can be found in Chaturantabut & Sorensen (Reference Chaturantabut and Sorensen2009).

Here, a set of entries $\boldsymbol{{\wp}}\subset \{1,2,\ldots ,l\}$ often called optimal (best) points are selected to determine $\hat{\boldsymbol{c}}$ by the following relation:

(2.15) $$\begin{eqnarray}\hat{\boldsymbol{c}}=(\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{U}}})^{-1}\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}(t)).\end{eqnarray}$$

Assuming that $\boldsymbol{F}^{\prime }$ is a componentwise function $\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}(t))=\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}(t))$ , we can rewrite (2.7) as

(2.16) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\hat{\boldsymbol{y}}(t)=(\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{A}\boldsymbol{{\mathcal{V}}})\hat{\boldsymbol{y}}(t)+\boldsymbol{{\mathcal{V}}}^{\text{T}}\boldsymbol{{\mathcal{U}}}(\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{U}}})^{-1}\boldsymbol{F}^{\prime }(\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{V}}}\hat{\boldsymbol{y}}(t)).\end{eqnarray}$$

In the present work, we perform the nonlinear POD on the same fluctuation matrix $\tilde{\boldsymbol{y}}_{m\times k}$ without separating the linear and nonlinear components. Consider the approximation of $\tilde{\boldsymbol{y}}$ as a nonlinear combination of the POD modes:

(2.17) $$\begin{eqnarray}\tilde{\boldsymbol{y}}(t)\approx \boldsymbol{{\mathcal{V}}}\unicode[STIX]{x1D73D}(t).\end{eqnarray}$$

The coefficients $\unicode[STIX]{x1D73D}(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 the $\boldsymbol{{\mathcal{P}}}^{\text{T}}\boldsymbol{{\mathcal{U}}}$ matrix is invertible. By using just the $\boldsymbol{{\wp}}$ rows of $\boldsymbol{{\mathcal{V}}}$ and $\boldsymbol{{\mathcal{U}}}$ , we can establish the relationship

(2.18) $$\begin{eqnarray}\boldsymbol{{\mathcal{V}}}_{{\wp}}\unicode[STIX]{x1D73D}(t)=\tilde{\boldsymbol{y}}_{{\wp}}(t),\end{eqnarray}$$

which further gives

(2.19) $$\begin{eqnarray}\tilde{\boldsymbol{y}}(t)\approx \boldsymbol{{\mathcal{V}}}\boldsymbol{{\mathcal{V}}}_{{\wp}}^{-1}\tilde{\boldsymbol{y}}_{{\wp}}.\end{eqnarray}$$

If the number of points used is greater than the number of significant modes, i.e.  $l>r$ , which is often the case, $\boldsymbol{{\mathcal{V}}}_{{\wp}}$ becomes a rectangular matrix. This makes the coefficients $\unicode[STIX]{x1D73D}(t)$ given by the gappy POD reconstruction

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D73D}(t)=\text{arg}\;\text{min}\Vert \tilde{\boldsymbol{y}}_{{\wp}}(t)-\boldsymbol{{\mathcal{V}}}_{{\wp}}\hat{\boldsymbol{a}}\Vert _{2},\quad \hat{\boldsymbol{a}}\in \mathbb{R}^{r}.\end{eqnarray}$$

The solution to the least squares problem (2.20) gives the result

(2.21) $$\begin{eqnarray}\tilde{\boldsymbol{y}}(t)\approx \boldsymbol{{\mathcal{V}}}\boldsymbol{{\mathcal{V}}}_{{\wp}}^{+}\tilde{\boldsymbol{y}}_{{\wp}},\end{eqnarray}$$

where $\boldsymbol{{\mathcal{V}}}_{{\wp}}^{+}$ is the Moore–Penrose pseudo-inverse of $\boldsymbol{{\mathcal{V}}}_{{\wp}}$ .

The POD-DEIM provides a way to introduce nonlinearity into the POD reconstructions; however, due to this nonlinear behaviour, 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 ensure an improvement in the result. Therefore, determining the optimal sizing of the low-dimensional representation is critical when using POD-DEIM for reconstruction. In the next section we present the FOM for generating high-dimensional data.

3.1 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 $Y$ -directions. The translational flow-induced vibration of a cylinder is strongly influenced by the four key non-dimensional parameters, namely mass ratio $(m^{\ast })$ , Reynolds number $(Re)$ , reduced velocity $(U_{r})$ , and critical damping ratio $(\unicode[STIX]{x1D701})$ , defined as

(3.1a-d ) $$\begin{eqnarray}m^{\ast }=\frac{M}{m_{f}},\quad Re=\frac{\unicode[STIX]{x1D70C}^{f}U_{\infty }D}{\unicode[STIX]{x1D707}^{f}},\quad U_{r}=\frac{U_{\infty }}{f_{n}D},\quad \unicode[STIX]{x1D701}=\frac{C}{2\sqrt{KM}},\end{eqnarray}$$

where $M$ is the mass of the body, $C$ and $K$ are the damping and stiffness coefficients respectively for an equivalent spring–mass–damper system of a vibrating structure, and  $U_{\infty }$ 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\unicode[STIX]{x03C0})\sqrt{K/M}$ and the mass of fluid displaced by the structure is $m_{f}=\unicode[STIX]{x1D70C}^{f}D^{2}L_{c}$ for a square cross-section, where $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 $\boldsymbol{M}=(m_{x},m_{y})$ with $m_{x}=m_{y}=M$ , the damping vector $\boldsymbol{C}=(c_{x},c_{y})$ with $c_{x}=c_{y}=C$ , and the stiffness vector $\boldsymbol{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

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle C_{L}=\frac{1}{\frac{1}{2}\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}DL_{c}}\int _{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D748}^{f}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{n}_{y}\,\text{d}\unicode[STIX]{x1D6E4}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle C_{D}=\frac{1}{\frac{1}{2}\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}DL_{c}}\int _{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D748}^{f}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{n}_{x}\,\text{d}\unicode[STIX]{x1D6E4}. & \displaystyle\end{eqnarray}$$

Here $\boldsymbol{n}_{x}$ and $\boldsymbol{n}_{y}$ are the Cartesian components of the unit outward normal $\boldsymbol{n}$ . In this study, we focus on the lift and drag forces due to the pressure field. Hence, we evaluate the pressure-induced drag ( $C_{Dp}$ ) and lift ( $C_{Lp}$ ) forces given by

(3.4a,b ) $$\begin{eqnarray}C_{Dp}=\frac{1}{\frac{1}{2}\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}D}\int _{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D748}_{p}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{n}_{x}\,\text{d}\unicode[STIX]{x1D6E4},\quad C_{Lp}=\frac{1}{\frac{1}{2}\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}D}\int _{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D748}_{p}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{n}_{y}\,\text{d}\unicode[STIX]{x1D6E4}.\end{eqnarray}$$

Figure 1(a) illustrates a schematic of the two-dimensional simulation domain used for the fluid–body interaction problem. The centre of the square column is located at the origin of the Cartesian coordinate system. The side length of the square column is denoted by $D$ . The distances to the upstream and 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_{\infty }$ 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 three-dimensional simulations.

Figure 1.

Full-order problem set-up for fluid–structure interaction: (a) schematic diagram of the computational domain and boundary conditions; (b) representative $Z$ -plane slice of the unstructured mesh. The top right inset displays the near-cylinder mesh.

3.2 Mesh convergence study

For the high-dimensional approximation of the FOM, the computational domain is discretized using an unstructured finite-element mesh, with a boundary layer mesh surrounding the body and three-node triangle (2-D) and six-node wedge (3-D) 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 1(b). Grid convergence study results are recorded in table 1 for the lock-in region. All cases for the mesh convergence are simulated at $Re=100$ , $m^{\ast }=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 (r.m.s.) of the transverse amplitude are analysed. It can be seen that values recorded for meshes M3 and M4 differ by less than 1 %. Therefore, 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, Guan & Jaiman Reference Miyanawala, Guan and Jaiman2016; Jaiman et al. Reference Jaiman, Guan and Miyanawala2016a ) and moderate- $Re$ (Jaiman et al. Reference Jaiman, Guan and Miyanawala2016a ; Miyanawala & Jaiman Reference Miyanawala and Jaiman2018) flows.

Table 1.

Grid convergence study at $Re=100$ , $m^{\ast }=3$ and $U_{r}=5.0$ .

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^{\ast },\unicode[STIX]{x1D701})=(100,3.0,0)$ . The snapshots of the FOM performed for the flow past a vibrating square cylinder are utilized to recover the POD modes and the DEIM points. The accuracy of the linear POD and POD-DEIM is systematically assessed with regard to their effectiveness in extracting the flow features.

4.1 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 ${\sim}99\,\%$ of the total contribution to the pressure fluctuations. The temporal coefficients ${\hat{y}}_{j}$ are determined by the $L^{2}$ inner product between the fluctuation matrix and the modes as expressed in (2.14). The mean pressure distribution and the first nine POD modes are displayed in figure 3. Note that, throughout this paper, the time-independent fields such as the mean pressure field and POD modes correspond to the initial flow field with zero bluff body motion. (We use a coordinate system fitted to the rigid body: the origin, $(X,Y,Z)=(0,0,0)$ , is always at the centre of the square cylinder.) The mean field is symmetric around the $X$ -axis along the wake centreline. 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 centreline, while modes 1, 3 and 9 are anti-symmetric with equal values and opposite signs about the wake centreline. As shown in figure 5, the POD time coefficients of these modes have the same frequency as the lift coefficient. It is evident that the first, third and ninth modes correspond predominantly to the Kármán vortex street with alternating positive and negative pressure regions about the $X$ -axis and the pressure contours resulting from a staggered vortex street. By examining figure 5, the time coefficients of the symmetric modes have the same frequency as the drag force. Of these modes, modes 2 and 8 have a high transverse gradient $(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y\gg \unicode[STIX]{x2202}/\unicode[STIX]{x2202}x)$ behaviour in the near-wake ( $0.5D$ – $5D$ ) region almost parallel to the top and bottom edges of the square cylinder, suggesting that these modes represent the influence from the shear layer. Modes 4, 5, 6 and 7 originate from the near-wake region and diffuse symmetrically towards the far wake. These modes have a dominant streamwise gradient $(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x\gg \unicode[STIX]{x2202}/\unicode[STIX]{x2202}y)$ compared to the transverse gradient. We can attribute these contributions to the near-wake bubble and its local dynamical property. For ease of explanation, we refer to these modes as the vortex shedding, the shear layer and the near-wake.

Figure 4.

Absolute value of the time-invariant contribution from each mode to the in-line ( $|F_{x}^{i}|$ ) and cross-flow ( $|F_{y}^{i}|$ ) forces. Modes 1, 3 and 9 capture the vortex shedding, modes 2 and 8 correspond to the effect of the shear layer and modes 4, 5, 6 and 7 capture the near-wake bubble effects.

Figure 4 quantifies the time-invariant contributions ( $F_{j}^{i}$ ) from each mode to the drag and lift forces. For the definition of $F_{j}^{i}$ , $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 the near-wake modes contribute entirely to the drag force. Further details of the force decomposition procedure using the modal contributions are presented in appendix B. Due to this directionally independent contribution of the bluff body features for the forces, the time coefficients ( ${\hat{y}}_{j}(t)$ ) of these modes should display the same frequencies of the lift and drag forces.

Figure 5.

Time history and FFT comparisons: (a) temporal coefficients of first five POD modes; (b) force coefficients. Note that modes 2, 4 and 5 have the same frequency ( $2f_{n}$ ) as the drag and modes 1 and 3 have the same frequency ( $f_{n}$ ) as the lift.

The time histories and the fast Fourier transform (FFT) spectra of the first five POD modes are shown in figure 5(a). The first and third mode coefficients have a low-frequency sinusoidal variation with the natural frequency ( $f_{n}$ ). The second, fourth and fifth mode coefficients have a non-zero mean with a frequency ${\approx}2f_{n}$ . Interestingly, as presented in figure 5(b), $f_{n}$ and $2f_{n}$ coincide with the frequencies of lift and drag, respectively. Hence, we can further confirm that modes 1 and 3 make their sole contribution to the fluctuating lift while modes 2, 4 and 5 contribute to the drag force. 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.

Figure 6.

Comparison between POD reconstruction and FOM result: pressure distribution obtained by (a) full-order model, (b) linear POD reconstruction, and (c) the relative error (%). Maximum error percentage is less than 2 %. The highly nonlinear near-wake region and the vortex cores have the highest error. The pressure values are normalized by $(\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2})/2$ . The flow is from left to right.

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 (2.13) for the pressure field gives $P(t)\approx \overline{P}+\sum _{j=1}^{r}{\hat{y}}_{j}(t)\boldsymbol{v}_{j}$ . Figure 6 illustrates the pressure distribution at $tU_{\infty }/D=100$ using the linear POD reconstruction. The recovered POD mode is compared with the result obtained from the FOM. A good match with a maximum relative local error of less than $2\,\%$ can be seen in figure 6. To quantify the accuracy of the entire flow field recovery, the normalized r.m.s. error of the entire distribution is considered. The r.m.s. error $\unicode[STIX]{x1D716}^{rms}$ is given by

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{rms}=\frac{\displaystyle \sqrt{\sum (P_{FOM}-P_{POD})^{2}/n_{c}}}{|\overline{P_{100}}|}\times 100,\end{eqnarray}$$

where $P_{FOM}$ and $P_{POD}$ are the pressure values of the mesh nodes extracted from the FOM and the POD reconstruction, respectively, $n_{c}$ is the node count of the mesh and $|\overline{P_{100}}|$ is the mean pressure of the field. When nine modes are used, this error is $\unicode[STIX]{x1D716}^{rms}=3.78\,\%$ . In this linear reconstruction, the highest error is observed at the regions known to exhibit nonlinear variation, such as the near-wake region, the shear layer and the vortex cores. Next, we analyse the POD-DEIM technique to improve the accuracy in these nonlinear flow features using the snapshot sequence and their respective DEIM points.

4.2 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 $\boldsymbol{{\mathcal{V}}}$ from the snapshot sequence, and a second basis $\boldsymbol{{\mathcal{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 results in adding contributions from some non-significant indices. In figure 7(a) it is clearly seen that for 100 DEIM points there are few mesh points lying away from the significant nonlinear region which are taken into the calculation. Further, in table 2 we quantify the number of points in the nonlinear wake region as we increase the number of DEIM points. It is clear that the percentage of points in the critical region decreases 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 7(b), the lowest $\unicode[STIX]{x1D716}^{rms}=4.05\,\%$ can be obtained when 70 DEIM points are used with seven 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 the nonlinear regions to be captured more accurately.

Figure 7.

(a) DEIM best points; the top right inset illustrates the near-wake DEIM points. (b) Performance of POD-DEIM compared with linear POD; solid lines denote the linear POD error. The least error is observed when 70 DEIM points are used with seven POD modes.

Table 2.

Distribution of DEIM points in the near cylinder wake.

Apart from the global and local accuracy, we further assess the computation time taken by the linear POD and nonlinear POD-DEIM reconstructions. The detailed analysis is presented in appendix C. Theoretically, the DEIM reconstruction process should be ${\approx}64$ times faster than linear POD reconstruction and the total DEIM process should be ${\approx}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, the DEIM has a speed-up of $3.98$ . In terms of accuracy, the 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 computation time while accurately capturing the flow features of the wake. Henceforth, we employ the POD-DEIM since it has improved accuracy when capturing the nonlinearities in the flow field at a lower computational cost.

Figure 8.

Comparison of POD-DEIM and FOM results: pressure distribution obtained by (a) FOM, (b) POD-DEIM reconstruction, and (c) the relative error of the reconstruction (%). Maximum error percentage is less than 1.5 %. The error in the highly nonlinear regions has reduced compared to the linear POD reconstruction. The pressure values are normalized by $1/2\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}$ . The flow is from left to right.

4.3 Drag and lift modes

In this section we analyse the behaviour of different modes in the near-cylinder region and explain the exclusive nature of their contributions to the pressure-induced drag and lift forces exerted on the oscillating cylinder in a uniform flow. Here, 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 the four quadrants. The variation is anti-symmetric about the streamwise centreline. In the maximum lift case (point B in figure 9 a), 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 (points A and C), the upstream and downstream pressure force differences tend to become equally strong and cancel each other out. The lift modes vary in such a way that the force on the top two quadrants is equal in magnitude and opposite in direction to the force on the bottom two quadrants. Due to this force cancellation, the vortex shedding (lift) modes make 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 9.

Reconstruction of the lift modes in the near-cylinder region for $U_{r}=6.0$ : (a) lift variation; (b) zero lift (increasing) at $tU_{\infty }/D=203.75$ ; (c) maximum lift at $tU_{\infty }/D=205.25$ ; (d) zero lift (decreasing) at $tU_{\infty }/D=206.85$ ; (e) minimum lift at $tU_{\infty }/D=208.50$ . The pressure values are normalized by $1/2\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}$ . The contours levels are from $-0.4$ to $0.4$ in increments of $0.1$ . The flow is from left to right.

Figure 10.

Reconstruction of the drag modes in the near-cylinder region for $U_{r}=6.0$ : (a) fluctuation of drag; (b) minimum drag fluctuation at $tU_{\infty }/D=205.05$ ; (c) zero drag fluctuation (increasing) at $tU_{\infty }/D=205.85$ ; (d) maximum drag fluctuation at $tU_{\infty }/D=206.70$ ; (e) zero drag fluctuation (decreasing) at $tU_{\infty }/D=207.40$ . The pressure values are normalized by $1/2\unicode[STIX]{x1D70C}^{f}U_{\infty }^{2}$ . The contours levels are from $-0.34$ to $0.16$ in increments of $0.025$ . The flow is from left to right.

Figure 10 describes the variation of drag modes with the fluctuation of the pressure-induced drag force. The drag fluctuation is defined as $C_{Dp}^{\prime }=C_{Dp}(t)-\overline{C_{Dp}}$ . The drag modes vary symmetrically around the wake centreline, hence make 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 positive. 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 force due to the pressure field on a moving bluff body based on the contributions from different POD modes can be expressed as

(4.2) $$\begin{eqnarray}F_{j}(t)=F_{j}^{0}+\mathop{\sum }_{i=1}^{n_{rj}}b_{j}^{i}(t)F_{j}^{i},\end{eqnarray}$$

where $F_{j}$ is the force in a particular direction ( $j=x$ for in-line and $j=y$ for transverse). $F_{j}^{0}$ is the time-independent contribution from the mean field and $F_{j}^{i}$ is the time-independent pressure fluctuation contribution calculated for the $i$ th mode. While $b_{j}^{i}(t)$ 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_{j}^{i}(t)$ for the streamwise and transverse forces as

(4.3) $$\begin{eqnarray}b_{j}^{i}(t)=\frac{F_{j}(t)-F_{j}^{0}}{F_{j}^{i}n_{rj}}.\end{eqnarray}$$

The complete description of this force decomposition and its usage is provided in appendix B.

Table 3.

Force contributions from the mean field and POD modes. $\overline{b_{j}^{i}}$ denotes the time-averaged force coefficient.

Table 3 summarizes the qualitative analysis of the contributions from the mean field and the modes to the pressure drag and lift forces. The mean field has a symmetric pressure distribution about the wake centreline, 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 make 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, but their contributions to the drag force are time-dependent. 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 (3.1). The variation of POD modes and their contribution to wake dynamics with reduced velocity is examined in the next section, with the aim of explaining the role of the wake features in sustaining the synchronized wake–body motion.