On the linear global stability analysis of rigid-body motion fluid–structure-interaction problems

Abstract

A rigorous derivation and validation for linear fluid–structure-interaction (FSI) equations for a rigid-body motion problem is performed in an Eulerian framework. We show that the ‘added stiffness’ terms arising in the formulation of Fanion et al. (Revue Européenne des Éléments Finis, vol. 9, issue 6–7, 2000, pp. 681–708) vanish at the FSI interface in a first-order approximation and can be neglected when considering the growth of infinitesimal disturbances. Several numerical tests with rigid-body motion are performed to show the validity of the derived formulation by comparing the time evolution between the linear and nonlinear equations when the base flow is perturbed by identical small-amplitude perturbations. In all cases both the growth rate and angular frequency of the instability matches within $0.1\,\%$ accuracy. The derived formulation is used to investigate the phenomenon of symmetry breaking for a rotating cylinder with an attached splitter plate. The results show that the onset of symmetry breaking can be explained by the existence of a zero frequency linearly unstable mode of the coupled FSI system. Finally, the structural sensitivity of the least stable eigenvalue is studied for an oscillating cylinder, which is found to change significantly when the fluid and structural frequencies are close to resonance.

1. Introduction

Fluid–structure-interaction (FSI) studies span a vast and diverse range of applications – from natural phenomenon such as the fluttering of flags (Shelley & Zhang 2011), phonation (Heil & Hazel 2011), blood flow in arteries (Freund 2014), path of rising bubbles (Ern et al. 2012), to the more engineering applications of aircraft stability (Dowell & Hall 2001), vortex induced vibrations (Williamson & Govardhan 2004), compliant surfaces (Riley, Gad-el Hak & Metcalfe 1988; Kumaran 2003) etc. The phenomena that emerge out of an FSI problem often exhibit a highly nonlinear, dynamically rich and complex behaviour with different flow regimes such as fluttering and tumbling of plates falling under gravity (Mittal, Seshadri & Udaykumar 2004), or the unsteady path of rising bubbles (Mougin & Magnaudet 2001; Ern et al. 2012). Despite its nonlinear nature, the initial transitions from one state to another may often be governed by a linear instability mechanism. The simplified case of a linear instability of a parallel boundary layer over a compliant surface was first derived and investigated by Benjamin (1959, 1960) and Landahl (1962). For boundary layers over Krammer-type compliant surfaces, the linear hydrodynamic instability study was performed by Carpenter & Garrad (1985, 1986) using this formulation and subsequently several studies have used the same formulation for studying flow instabilities over flexible surfaces, for example, Carpenter & Morris (1990) and Davies & Carpenter (1997).

When the parallel flow assumption is no longer valid a global approach needs to be taken. From the perspective of global instabilities in FSI problems, the first such study was reported by Cossu & Morino (2000) for the instability of a spring-mounted cylinder. The authors considered a rigid two-dimensional cylinder that was free to oscillate in the cross-stream direction, subject to the action of a spring-mass-damper system. The linearized equations were solved in a non-inertial frame of reference attached to the cylinder. For low (solid to fluid) mass density ratios, the authors report that the critical Reynolds number for vortex shedding drops to $Re_{c}\approx 23$. Mittal & Singh (2005) also reported a similarly low critical $Re$ for a cylinder oscillating in both streamwise and transverse directions. A non-inertial reference is also used by Navrose & Mittal (2016) to study the lock-in phenomenon of cylinders oscillating in cross-stream through the linear stability analysis. More recently, Magnaudet and co-workers approached the problem of rising and falling bodies through the perspective of linear instability. Fabre, Assemat & Magnaudet (2011) developed a quasi-static model to describe the stability of heavy bodies falling in a viscous fluid. Assemat, Fabre & Magnaudet (2012) performed a linear study of thin and thick two-dimensional plates falling under gravity. The problem is again formulated in a non-inertial frame of reference of the moving body undergoing both translation and rotation. The authors showed a quantitative agreement between the quasi-static model of Fabre et al. (2011) and the linear stability results for high mass-ratio cases, although the agreement systematically deteriorated as the mass ratio was reduced. Tchoufag, Fabre & Magnaudet (2014a) performed similar studies for three-dimensional disks and thin cylinders, Tchoufag, Magnaudet & Fabre (2014b) applied the linear stability analysis to spheroidal bubbles and Cano-Lozano et al. (2016) to oblate bubbles.

In general, linear FSI investigations have largely followed one of two methodologies. Either through the parallel flow approach or through frames of reference attached to the rigid body in motion. There are a few exceptions to this. Lesoinne et al. (2001) formulated the linear FSI problem in the arbitrary-Eulerian–Lagrangian (ALE) form of the inviscid Navier–Stokes where they treated the grid velocity as a pseudo-variable. Fanion, Fernández & Le Tallec (2000) derived a more general formulation starting from the ALE equations in the weak form, which has been used by Fernández & Le Tallec (2003a,b). The authors derive a formulation independent of the fluid grid velocity, but both the boundary conditions and the fluid stresses at the FSI interface are modified. The velocity continuity boundary condition transforms to a transpiration boundary condition while additional stress terms at the interface are obtained comprising of higher-order derivatives of the base flow. These additional stresses have been termed as ‘added stiffness’ terms (Fernández & Le Tallec 2003a,b). Goza, Colonius & Sader (2018) used an immersed boundary method for the global stability analysis of inverted flag flapping. In a very recent development, Pfister, Marquet & Carini (2019) also derived the linear FSI formulation with the ALE framework. In their formulation the fluid and structural quantities are solved on the moving material points. The standard Navier–Stokes and divergence equations are therefore modified to account for the motion of the material points. Using this formulation on moving material points, the authors have investigated several different FSI problems including ones with nonlinear structural models, elastic flutter instabilities and finite aspect ratio structures (Pfister et al. 2019; Pfister & Marquet 2020).

In the current work we follow the methodology of Fanion et al. (2000) and Fernández & Le Tallec (2003a) which, in spirit, is similar to the methodology used by Benjamin (1960) and Landahl (1962) for compliant wall cases. However, unlike Fanion et al. (2000) we proceed with the linearization of the equations in their strong form. The final expressions for the linearized equations are all evaluated on a stationary grid (Eulerian formulation). The resulting expressions for the linearized FSI problem are very similar to the ones obtained by Fanion et al. (2000); however, some crucial differences arise for the stresses at the interface. The ‘added stiffness’ terms arising in the work of Fanion et al. (2000) are shown to vanish at first order. While we invoke no special assumptions for the linearization of the fluid equations (other than small-amplitude perturbations), the current capabilities of our numerical solver limits the class of FSI problems that can be validated. Therefore, we confine the focus of the current work to FSI problems undergoing rigid-body motion and defer a more general formulation and validation to future work. The derived formulation for rigid-body linear FSI is numerically validated by comparing linear and nonlinear evolution of different cases which are started from the same base flow state, perturbed by identical small-amplitude disturbances. It is shown that the linear and nonlinear simulations evolve nearly identically through several orders of magnitude of growth of the perturbations. The nonlinear cases eventually reach a saturated state while the linear simulations continue with the exponential growth. The derived equations are then used to investigate the instability of an oscillating cylinder at subcritical Reynolds numbers, for an ellipse in both pure translation and pure rotation, and for the case of symmetry breaking in a cylinder with an attached splitter plate. Finally, in § 4 we investigate the structural sensitivity of the least stable eigenvalue in the coupled FSI problem of an oscillating cylinder at $Re=50$, with varying structural parameters.

The remainder of the paper is organized as follows. In § 2 we describe the problem in a general setting and derive the linearized equations for a FSI problem for rigid-body motion. Numerical validation and results for the derived formulation are presented in § 3. In § 4 we introduce the adjoint problem for the linear FSI system of a cylinder oscillating in cross-flow and show the changes in structural sensitivity of the unstable eigenvalue. Section 5 concludes the paper.

2. Linearization of FSI

2.1. General problem description

We primarily follow the index notation with the implied Einstein summation. Bold face notation is used to represent vectors where necessary. Consider a FSI problem as illustrated in figure 1. The bounded region in white marked by $\varOmega ^{f}$ represents the fluid part of the domain and the grey region marked by $\varOmega ^{s}$ represents the structural domain. The combined fluid and structural regions are bounded by the boundaries represented by $\partial \varOmega ^{f}_{v}$ and $\partial \varOmega ^{f}_{o}$. Here $\partial \varOmega ^{f}_{v}$ represents the far-field boundary conditions which define the problem and $\partial \varOmega ^{f}_{o}$ represents the open boundary condition applied for computations performed on a finite sized domain. The structural domain is bounded by the time varying FSI interface $\varGamma (t)$ on which the fluid forces act. The Navier–Stokes equations, along with the incompressibility constraint govern the evolution of the fluid in $\varOmega ^{f}$. For a problem with moving interfaces, the Navier–Stokes is usually formulated in the ALE formulation (Ho & Patera 1990, 1991) defined on moving material points as

(2.1a)\begin{gather} \left.\frac{\partial U_{i}}{\partial t}\right\vert_{\boldsymbol{W^{g}}} + (U_{j} -W^{g}_{j}) \frac{\partial U_{i}}{\partial x_{j}} = -\frac{\partial P}{\partial x_{i}} + \frac{1}{Re} \frac{\partial^{2} U_{i}}{\partial x_{j}\partial x_{j}}, \end{gather}

(2.1b)\begin{gather}\frac{\partial U^{0}_{i}}{\partial x_{i}} = 0 . \end{gather}

Here $\boldsymbol {U}$ and $P$ represents the fluid velocity and pressure, respectively, while $\boldsymbol {W^{g}}$ represents the velocity of the material points. We refer to the time derivative in this formulation $(\partial U_{i}/\partial t\vert _{\boldsymbol {W^{g}}})$ as the ALE time derivative evaluated when following a material point with velocity $\boldsymbol {W^{g}}$. Note that $\boldsymbol {W^{g}}$ is not uniquely defined. The only restriction on material point velocity is for the points on the interface $\varGamma (t)$, where the fluid velocity is equal to the velocity of the interface points. Typically $\boldsymbol {W^{g}}$ is defined by some function which smoothly extends the velocity at the interface to the rest of the fluid domain.

Figure 1. Domain for a general FSI problem. The equilibrium position of the interface is marked by $\varGamma ^{0}$ and the perturbed position is marked by $\varGamma (t)$.

The equations governing the structural motion depend on the type of modelling or degrees of freedom of the structure being considered. Confining the discussion to rigid-body motion problems, we represent the structural equations in a general form as

(2.2)\begin{equation} \mathcal{M}\frac{\partial^{2} \eta_{i}}{\partial t^{2}} + \mathcal{D} \frac{\partial \eta_{i}}{\partial t} + \mathcal{K} \eta_{i} = \mathcal{F}_{i} . \end{equation}

Here $\mathcal {M}$ is a generalized inertia, $\mathcal {D}$ is a generalized damping, $\mathcal {K}$ is a generalized stiffness and $\mathcal {F}_{i}$ represents the fluid forces acting on $\varGamma (t)$. The equation represents a typical linear spring-mass-damper system for rigid-body motion. The structural degrees of freedom are represented by $\eta _{i}$, whose definition depends on the type of modelling used for the structure. For example, for a cylinder free to oscillate in one direction subject to a spring-damper action (as in Cossu & Morino 2000), $\eta _{i}$ represents the position of the centre of mass in that direction, $\mathcal {M}$ is the cylinder mass, $\mathcal {D}$ is the damping coefficient, $\mathcal {K}$ is the spring stiffness constant and $\mathcal {F}_{i}$ is the integral of the fluid forces acting on the cylinder in the $i$th direction.

The fluid and structural domains are coupled at the FSI interface through a no-slip and no penetration condition. Defining $\boldsymbol {x}^{\varGamma }$ as the (time-dependent) position of points on the interface $\varGamma (t)$ and $\boldsymbol {v}^{\varGamma }$ as their instantaneous velocity, we may write the boundary condition as

(2.3)\begin{equation} W^{g}_{i} = U_{i} = \frac{\textrm{d} x_{i}^{\varGamma}}{\textrm{d} t}= v^{\varGamma}_{i}\quad \text{on } \varGamma(t). \end{equation}

The coupled problem is now defined by (2.1a) and (2.1b) in the fluid domain, (2.2) in the structural domain, coupled with the no-slip condition at the interface (2.3). The problem is completed by the application of appropriate Dirichlet boundary conditions on $\partial \varOmega ^{f}_{v}$ and open boundary conditions on $\partial \varOmega ^{f}_{o}$. We assume that the outer domain boundaries always remain stationary.

2.2. Steady state

For the problem described in § 2.1, consider a full system state given by $(\boldsymbol {U}^{0},P^{0},\boldsymbol {{\eta }}^{0})$ (referred to as the base flow state), defined on material points $\boldsymbol {x}^{0}$, which satisfies the time-independent incompressible Navier–Stokes and structural equations (along with the respective boundary conditions), i.e.

(2.4a)\begin{gather} U^{0}_{j}\frac{\partial U^{0}_{i}}{\partial x_{j}} +\frac{\partial P^{0}}{\partial x_{i}} - \frac{1}{Re}\frac{\partial^{2} U^{0}_{i}}{\partial x_{j}\partial x_{j}} = 0, \end{gather}

(2.4b)\begin{gather}\frac{\partial U^{0}_{i}}{\partial x_{i}} = 0 , \end{gather}

(2.4c)\begin{gather}\mathcal{K} \eta^{0}_{i} - \mathcal{F}^{0}_{i} = 0, \end{gather}

(2.4d)\begin{gather}U_{i} = v^{\varGamma}_{i} = 0;\quad \text{on } \varGamma^{0}. \end{gather}

Here $\mathcal {F}^{0}_{i}$ are the fluid forces acting on the structure at equilibrium and $\varGamma ^{0}$ is the equilibrium position of the FSI interface.

2.3. Linearization of the structure

We begin with the linearization of the structural equations. In defining the generalized structural equation (2.2) we make an implicit assumption that the structural equations are inherently linear, which we consider appropriate for the class of rigid-body motion problems that are being considered. In principle, nonlinear spring-mass-damper systems may also be considered but for the moment we restrict the discussion to linear spring-mass-damper systems. Given an initial stationary solution $\eta ^{0}_{i}$ of the structural equation, we may introduce a perturbed state given by the superposition of the stationary state and small-amplitude perturbations, i.e. $\eta _{i} = \eta ^{0}_{i} + \eta '_{i}$. Similarly, we decompose the fluid forces acting on the structure as $\mathcal {F}_{i} = \mathcal {F}^{0}_{i} + \mathcal {F}'_{i}$, where $\mathcal {F}'_{i}$ are the fluid forces acting on the structure due to the (yet unknown) linearized fluid perturbations. Introducing the decomposition into (2.2), we obtain the equations for the structural perturbations

(2.5)\begin{equation} \mathcal{M}\frac{\partial^{2} \eta'_{i}}{\partial t^{2}} + \mathcal{D} \frac{\partial\eta'_{i}}{\partial t} + \mathcal{K} \eta^{0}_{i} + \mathcal{K} \eta'_{i} - \mathcal{F}^{0}_{i} - \mathcal{F}'_{i} = 0. \end{equation}

Equation (2.5) will be completed if we can evaluate the expression for the term $\mathcal {F}'_{i}$ arising due to the linearized fluid perturbations. The two systems (fluid and structure) are coupled through the term $\mathcal {F}'_{i}$ and the velocity boundary conditions at the interface.

2.4. Taylor expansion based linearization of the fluid equations

To begin with, we define an affine mapping between the material points in the reference (equilibrium) configuration and the perturbed points as

(2.6a)\begin{gather} \boldsymbol{x} = [\mathcal{I} + \mathcal{R'}]\boldsymbol{x}^{0} + b, \end{gather}

(2.6b)\begin{gather}\therefore \boldsymbol{{\rm \Delta} x} = \mathcal{R'}\boldsymbol{x}^{0} + b. \end{gather}

Here $b$ is a vector representing a constant translation of the material points, $\mathcal {I}$ is the identity matrix and $\mathcal {R}'$ can be considered to be a ‘perturbation matrix’ for the material points. Using such a decomposition allows us to conveniently express the inverse mapping

(2.7)\begin{equation} \frac{\partial \boldsymbol{x}^{0}}{\partial \boldsymbol{x}} = \left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{x}^{0}}\right)^{-1} = \left[\mathcal{I} + \mathcal{R'}\right]^{-1}. \end{equation}

For a small deformation, we may assume that $||\mathcal {R'}||\ll 1$ and the matrix inverse may then be written as

(2.8)\begin{equation} \left[\mathcal{I} + \mathcal{R'}\right]^{-1} = \mathcal{I} - \mathcal{R'} + \mathcal{R'}^{2} - \mathcal{R'}^{3} \cdots. \end{equation}

Retaining only the first-order term of the expansion, we obtain the approximate inverse as

(2.9)\begin{equation} \frac{\partial \boldsymbol{x}^{0}}{\partial \boldsymbol{x}} = \left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{x}^{0}}\right)^{-1} \approx \left[\mathcal{I} - \mathcal{R'}\right]. \end{equation}

We refer to (2.9) as the geometric linearization of the problem. This allows us to conveniently reformulate all derivative terms evaluated on the perturbed grid in terms of the derivatives in the reference configuration. In addition, since $\mathcal {R}'$ is a perturbation matrix, this allows us to identify higher-order terms arising due to this geometric nonlinearity and later discard them when we retain only the first-order terms. In what follows, all quantities of the form $\partial /\partial x_{j}^{0}$ represent the evaluation of derivatives in the reference configuration.

Next we define a perturbation field for the fluid components $(\boldsymbol {u}',p')$. The perturbation field is also defined on the same equilibrium grid $\boldsymbol {x}^{0}$ on which the base flow $(\boldsymbol {U}^{0},P^{0})$ has been defined. The total velocity and pressure fields on the perturbed locations are then evaluated using a superposition of the first-order Taylor expansions for the base flow and the perturbation field, i.e.

(2.10a)\begin{gather} U_{i}(x,y,z) = U^{0}_{i} + \left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right) + u'_{i} + \left(\frac{\partial u'_{i}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right), \end{gather}

(2.10b)\begin{gather}P(x,y,z) = P^{0} + \left(\frac{\partial P^{0}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right) + p' + \left(\frac{\partial p'}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right) . \end{gather}

The last term in (2.10a) and (2.10b) is dropped since it represents a second-order quantity arising due to the interaction of fluid and geometric perturbation terms. This results in the expressions for the total velocity and pressure as

(2.11a,b)\begin{equation} u^{\xi}_{i}(x,y,z) = \left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right),\quad p^{\xi}(x,y,z) = \left(\frac{\partial P^{0}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right), \end{equation}

(2.12a)\begin{gather} U_{i}(x,y,z) = U^{0}_{i} + u^{\xi}_{i} + u'_{i} , \end{gather}

(2.12b)\begin{gather}P(x,y,z) = P^{0} + p^{\xi} + p' . \end{gather}

Where (2.12a) and (2.12b) are derived from (2.10a) and (2.10b), respectively, after dropping the second-order terms, and using the compact notations $u^{\xi }_{i}$ and $p^{\xi }$ to represent the first-order quantities of the Taylor expansion of the base flow velocity and pressure, respectively. Physically, this means that, for a point moving in space, we ignore changes in the perturbation field experienced due to the motion of the point, however, we include the first-order changes in the base flow field. The expressions in (2.12a) amount to a triple decomposition of the total velocity field, where the three terms $U^{0}_{i},u^{\xi }_{i}$ and $u'_{i}$ signify the base flow field, the perturbation of the base field observed by a point due to its motion and the perturbation velocity field, respectively. Additionally, associated with the motion of the material points, a velocity of the material points can also be defined such that

(2.13)\begin{equation} w_{i} = \frac{\partial x_{i}}{\partial t} = \frac{\partial {\rm \Delta} x_{i}}{\partial t}. \end{equation}

Taking the time derivative of (2.12a), we obtain the ALE time derivative of the total velocity along the trajectory of the material point motion

(2.14)\begin{align} \left.\frac{\partial U_{i}}{\partial t}\right\vert_{\boldsymbol{w}} &= \frac{\partial U^{0}_{i}}{\partial t} + \frac{\partial }{\partial t} \left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}{\rm \Delta} x_{j} \right) +\frac{\partial u'_{i}}{\partial t} \nonumber\\ &\implies\left.\frac{\partial U_{i}}{\partial t}\right\vert_{\boldsymbol{w}} = \left(w_{j}\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}\right) +\frac{\partial u'_{i}}{\partial t}. \end{align}

Next we substitute the triple decomposition of the velocity and pressure fields into the governing equations for the fluid motion and perform the geometric linearization so that all derivative quantities are consistently evaluated based on the reference configuration. Starting with the divergence-free constraint at the perturbed configuration leads to the following set of expressions:

(2.15)\begin{equation} \frac{\partial U_{i}}{\partial x_{i}} =\left\lbrace \begin{aligned} & \underbrace{\frac{\partial U^{0}_{i}}{\partial x^{0}_{i}}}_{I} + \underbrace{{\rm \Delta} x_{j}\frac{\partial }{\partial x^{0}_{j}} \left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{i}}\right)}_{II} +\underbrace{\frac{\partial u'_{i}}{\partial x^{0}_{i}}}_{III} \\ & \quad -\underbrace{\frac{\partial u'_{i}}{\partial x^{0}_{k}}\mathcal{R'}_{ki} }_{IV} -\underbrace{\frac{\partial }{\partial x^{0}_{k}}\left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}\right) {\rm \Delta} x_{j}\mathcal{R'}_{ki}}_{V}-\underbrace{\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} \frac{\partial {\rm \Delta} x_{j}}{\partial x^{0}_{k}}\mathcal{R'}_{ki}}_{VI}. \end{aligned}\right. \end{equation}

Term $I$ vanishes since it represents the divergence-free constraint of the base flow in the reference configuration. Similarly, term $II$ contains the same divergence-free constraint inside the brackets and hence also vanishes. Terms $IV,V$ and $VI$ all represent terms that are second order in the perturbations quantities $\boldsymbol {u}',\boldsymbol {{\rm \Delta} x}$ and $\mathcal {R'}$. Therefore, to a first-order approximation these quantities can be dropped. This leaves the final divergence constraint as

(2.16)\begin{equation} \frac{\partial U_{i}}{\partial x_{i}} = \frac{\partial u'_{i}}{\partial x^{0}_{i}} = 0. \end{equation}

Note that the derivative is evaluated in the reference configuration (but satisfies the divergence-free constraint on the perturbed locations up to a first-order approximation).

In a similar manner, we may introduce the triple decomposition of the velocity and pressure field into the ALE form of the Navier–Stokes evaluated at the perturbed locations as

(2.17)\begin{align} & \left(U^{0}_{j}\frac{\partial U^{0}_{i}}{\partial x_{j}} + U^{0}_{j}\frac{\partial u^{\xi}_{i}}{\partial x_{j}} + u^{\xi}_{j}\frac{\partial U^{0}_{i}}{\partial x_{j}}\right) + \left( U^{0}_{j}\frac{\partial u'_{i}}{\partial x_{j}} + u'_{j}\frac{\partial U^{0}_{i}}{\partial x_{j}}\right) +\left( \left.\frac{\partial U_{i}}{\partial t}\right\vert_{\boldsymbol{w}} - w_{j}\frac{\partial U_{i}}{\partial x_{j}} \right) \nonumber\\ &\quad =\left(-\frac{\partial P^{0}}{\partial x_{i}} + \frac{1}{Re}\frac{\partial^{2}U^{0}_{i}}{\partial x_{j}\partial x_{j}}\right) + \left(-\frac{\partial p^{\xi}}{\partial x_{i}} + \frac{1}{Re}\frac{\partial^{2}u^{\xi}_{i}}{\partial x_{j}\partial x_{j}}\right) +\left(-\frac{\partial p'}{\partial x_{i}} + \frac{1}{Re}\frac{\partial^{2}u'_{i}}{\partial x_{j}\partial x_{j}}\right). \end{align}

The nonlinear transport terms have already been dropped. We now substitute the ALE time derivative from (2.14) and introduce the geometric linearization to consistently evaluate derivatives based on the original configuration. A full expansion of all the terms can be found in appendix A. We write the final form of the Navier–Stokes obtained after the expansion and simplification of the terms and dropping all terms higher than first order:

(2.18)\begin{align} & \left[U^{0}_{j}\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} + \frac{\partial P^{0}}{\partial x^{0}_{i}} -\frac{1}{Re}\frac{\partial^{2} U^{0}_{i}}{\partial x^{0}_{j}\partial x^{0}_{j}}\right] + {\rm \Delta} x_{l}\frac{\partial }{\partial x^{0}_{l}}\left[U^{0}_{j}\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} + \frac{\partial P^{0}}{\partial x^{0}_{i}} -\frac{1}{Re}\frac{\partial^{2} U^{0}_{i}}{\partial x^{0}_{j}\partial x^{0}_{j}}\right] \nonumber\\ &+\quad \left[\frac{\partial u'_{i}}{\partial t} +U^{0}_{j}\frac{\partial u'_{i}}{\partial x^{0}_{j}} + u'_{j}\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} +\frac{\partial p'}{\partial x^{0}_{i}} - \frac{1}{Re}\frac{\partial^{2} u'_{i}}{\partial x^{0}_{j}\partial x^{0}_{j}}\right] = 0. \end{align}

The terms in the first square bracket may be identified as the steady-state equation for the base flow in the reference configuration. The terms inside the second square bracket are also the steady-state equations in the reference configuration. In fact, the set of terms from the first two brackets together represent the first-order Taylor expansion of the steady-state equations for the base flow at the perturbed points. Thus, all the terms in the first two brackets vanish. In hindsight, the final set of expressions obtained seems rather obvious. In fact one may observe the first three terms in the expansion of the divergence-free constraint and realize that they amount to a similar expression. Terms $I$ and $II$ in (2.15) together form the first-order Taylor expansion of the base flow divergence-free constraint at the perturbed locations. Since the solution of the steady-state equations and the base field divergence is identically zero everywhere, their Taylor expansions also vanish at the perturbed locations. Thus, we are left with the final linear equations for the perturbed quantities which is also independent of the arbitrarily defined velocity of the grid motion $\boldsymbol {w}$:

(2.19)\begin{equation} \frac{\partial u'_{i}}{\partial t}+U^{0}_{j}\frac{\partial u'_{i}}{\partial x^{0}_{j}} +u'_{j}\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}}+\frac{\partial p'}{\partial x^{0}_{i}} -\frac{1}{Re}\frac{\partial^{2} u'_{i}}{\partial x^{0}_{j}\partial x^{0}_{j}}= 0.\end{equation}

2.5. Linearized boundary conditions

Before proceeding to evaluate the total linearized forces on the perturbed points, we consider the global balance of fluxes for the base flow at the equilibrium and the perturbed positions. This allows us to evaluate the forces arising due to the variation of the boundary in a steady base flow. We use the conservative form of the equations, starting with the base flow state in the equilibrium configuration and evaluate the integral over the whole domain:

(2.20)\begin{align} & \frac{\partial}{\partial x^{0}_{j}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij}) = 0 \nonumber\\ &\quad \implies \int_{\varOmega^{0}}\frac{\partial}{\partial x^{0}_{j}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})\,\textrm{d}\varOmega = 0 \nonumber\\ &\quad \implies \int_{\partial \varOmega^{f}_{v}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n^{0}_{j}\,\textrm{d} S + \int_{\partial \varOmega^{f}_{o}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n^{0}_{j}\,\textrm{d} S - \int_{\varGamma^{0}}\sigma^{0}_{ij}n^{0}_{j}\,\textrm{d} S = 0. \end{align}

Here $n^{0}_{j}\,\textrm {d} S$ represents the local surface vector in the equilibrium configuration. For compact notation, we have denoted the stresses due to the fluid at equilibrium as $\sigma ^{0}_{ij}$, defined as

(2.21)\begin{equation} \sigma^{0}_{ij} = -P^{0}\delta_{ij} +\frac{1}{Re}\left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} + \frac{\partial U^{0}_{j}}{\partial x^{0}_{i}} \right) . \end{equation}

The three integrals in (2.20) represent the momentum flux through the far field ($\partial \varOmega ^{f}_{v}$), outflow ($\partial \varOmega ^{f}_{o}$) and FSI interface ($\varGamma ^{0}$), respectively. Of course it is not necessary to integrate over the whole domain, and the integral may be performed over any arbitrarily chosen part of the domain. The use of divergence theorem will then provide us with the flux balance across the surfaces of this arbitrarily chosen volume. Accordingly, we may choose an arbitrary volume enclosed by three surfaces. Two of the surfaces of this arbitrary volume are chosen such that they coincide with $\partial \varOmega _{v}$ and $\partial \varOmega _{o}$ and a third surface may be chosen arbitrarily, which we denote as $C$. Performing the volume integral over this arbitrary volume and using the divergence theorem gives the relation between the fluxes across the boundaries,

(2.22)\begin{equation} \int_{\partial \varOmega^{f}_{v}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n^{0}_{j}\,\textrm{d} S + \int_{\partial \varOmega^{f}_{o}} (U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n^{0}_{j}\,\textrm{d} S + \int_{C}(U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n_{j}\,\textrm{d} S = 0. \end{equation}

We highlight here that at this point we are not calculating the balance of total momentum in a flow field that has been perturbed. But rather, we are evaluating the expressions for the momentum flux, across different enclosing surfaces, for the same steady base flow. Comparing (2.20) and (2.22) it is easy to see that the fluxes through the arbitrary surface $C$ and the FSI interface $\varGamma ^{0}$ must be equal, i.e.

(2.23)\begin{equation} \int_{C}(U^{0}_{i}U^{0}_{j} - \sigma^{0}_{ij})n_{j}\,\textrm{d} S = \int_{\varGamma^{0}}(-\sigma^{0}_{ij})n^{0}_{j}\,\textrm{d} S. \end{equation}

Up to this point no assumptions have been invoked and (2.23) is simply a result of momentum flux conservation. The surface $C$ is arbitrary and can be deformed in any desired fashion. Therefore, the surface can be assumed to be a perturbation of the equilibrium FSI interface $\varGamma ^{0}$ (for example, the perturbed interface obtained during an unsteady FSI problem). For the exact evaluation of the base-flow momentum flux across the perturbed surfaces, the exact values at the perturbed points would be required. However, following the Taylor expansion based procedure, these may be evaluated using the truncated Taylor series. Noting that $\boldsymbol {U^{0}}$ is identically zero on $\varGamma ^{0}$, the resulting linearized expression for the surface integral over this perturbed interface may be evaluated as

(2.24)\begin{align} & \int_{\varGamma^{0}}\sigma^{0}_{ij}n^{0}_{j}\,\textrm{d} S + \int_{\varGamma^{0}}\sigma^{0}_{ij}n'_{j}\,\textrm{d} S + \int_{\varGamma^{0}}\sigma^{\xi}_{ij}n^{0}_{j}\,\textrm{d} S = \int_{\varGamma^{0}}\sigma^{0}_{ij}n^{0}_{j}\,\textrm{d} S \nonumber\\ &\quad \implies \int_{\varGamma^{0}}\big[\sigma^{0}_{ij}n'_{j} + \sigma^{\xi}_{ij}n^{0}_{j}\big]\,\textrm{d} S = 0. \end{align}

Here $n'_{j}\,\textrm {d} S$ represents the first-order change in the surface vector due to the perturbation of the boundary and we have dropped the second-order terms for convection $u^{\xi }_{i}u^{\xi }_{j}$ and the surface force $\sigma ^{\xi }_{ij}n'_{j}$. The expression $\sigma ^{\xi }_{ij}$ is simply the first-order term of the Taylor expansion of the base flow stresses, defined as

(2.25)\begin{equation} \sigma^{\xi}_{ij} = {\rm \Delta} x_{k}\frac{\partial }{\partial x^{0}_{k}} \left[-P^{0}\delta_{ij} +\frac{1}{Re}\left(\frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} + \frac{\partial U^{0}_{j}}{\partial x^{0}_{i}} \right)\right] . \end{equation}

The first term in (2.24) represents the variation of the base flow forces due to the change in surface normal and the second term represents the variation due to the change in boundary position. To a first-order approximation, these terms balance to zero. Note that no assumption has been made on the type of deformation at the boundary and the condition holds for any arbitrary deformation of the boundary. We note that these are precisely the terms that arise in Fanion et al. (2000) and Fernández & Le Tallec (2003a,b) that have been termed as ‘added stiffness’. To a first-order approximation, they simply sum up to zero and play no role in the linear dynamics. However, we note that it is the integral of the added stiffness terms that vanishes and not necessarily the point wise values. One may now evaluate the total linearized forces arising from (2.18) on the perturbed FSI interface and obtain the linearized boundary conditions for the stress balance (2.5) as

(2.26)\begin{align} & \left.\begin{aligned} & \mathcal{M}\frac{\partial^{2} \eta'_{i}}{\partial t^{2}} + \mathcal{D}\frac{\partial\eta'_{i}}{\partial t} + \mathcal{K} (\eta^{0}_{i} + \eta'_{i}) + \int_{\varGamma^{0}}\sigma^{0}_{ij}n^{0}_{j}\,\textrm{d} S \nonumber\\ & \quad + \int_{\varGamma^{0}}\sigma^{0}_{ij}n'_{j}\,\textrm{d} S + \int_{\varGamma^{0}}\sigma^{\xi}_{ij}n^{0}_{j}\,\textrm{d} S + \int_{\varGamma^{0}}\sigma'_{ij}n^{0}_{j}\,\textrm{d} S \end{aligned}\right\rbrace = 0, \nonumber\\ &\quad \implies \mathcal{M}\frac{\partial^{2} \eta'_{i}}{\partial t^{2}} + \mathcal{D}\frac{\partial\eta'_{i}}{\partial t} + \mathcal{K}\eta'_{i} + \int_{\varGamma^{0}}\sigma'_{ij}n^{0}_{j}\,\textrm{d} S = 0, \end{align}

where $\sigma '_{ij}$ is the fluid stress due to the perturbation field $(\boldsymbol {u}',p')$ defined as

(2.27)\begin{equation} \sigma'_{ij} = -p'\delta_{ij} +\frac{1}{Re}\left(\frac{\partial u'_{i}}{\partial x^{0}_{j}} + \frac{\partial u'_{j}}{\partial x^{0}_{i}} \right) . \end{equation}

The velocity continuity boundary conditions on the moving interface $\varGamma (t)$ simply become

(2.28)\begin{align} & U_{i} = v^{\varGamma}_{i}\quad \textrm{on}\ \varGamma(t) \nonumber\\ &\quad \implies {\rm \Delta} x^{\varGamma}_{k}\frac{\partial U^{0}_{i}} {\partial x_{k}} + u'_{i} = v^{\varGamma}_{i}. \end{align}

Thus, the perturbation velocities must account for the Taylor expansion term of the base flow at the perturbed boundary. Fanion et al. (2000) and Fernández & Le Tallec (2003a,b) refer to this as the transpiration boundary condition.

The linear FSI problem for the perturbations is turned into a standard linear perturbation problem with the exception of the modified boundary conditions (transpiration instead of no-slip). The additional forces at the perturbed boundary are simply due to the perturbation field $(\boldsymbol {u}',p')$ and no ‘added stiffness’ terms arise in a linear approximation. One may immediately notice that the final form is a generalization of the form derived by Benjamin (1959, 1960) and Landahl (1962) for parallel flow cases. Finally, we highlight an assumption that is inherent throughout the mathematical derivation. When evaluating the base flow terms at the perturbed points, a first-order Taylor expansion has been used. This is perfectly valid when the perturbed points are inside the equilibrium fluid domain however, when the perturbed points move outside this domain, the validity of the Taylor expansion is less clear. The implicit validity of the Taylor expansion beyond the equilibrium fluid domain underpins the mathematical results of both the unsteady problem in § 2.4 as well as for the steady problem in the current subsection. To the best of our knowledge,we are unaware of any previous work which may resolve this conundrum. We provide a heuristic answer in the next section where several examples are considered for rigid-body motion in symmetric and asymmetric configurations to show the validity of the current approach.

3. Linear instability results

3.1. Numerical method

Numerical tests are performed to validate the linear equations for FSI derived in the previous section. The cases considered are an oscillating cylinder with a spring-damper action, oscillating and rotating ellipse initially held at an angle to the flow, and rotating cylinder-splitter body. All computations were performed using a high-order spectral-element method code (Fischer, Lottes & Kerkemeier 2008). The code uses $n$th-order Lagrange interpolants at Gauss–Lobatto–Legendre (GLL) points for the representation of the velocity and $(n-2)$th-order interpolants at Gauss–Legendre points for the representation of the pressure in a $\mathbb {P}_{n}-\mathbb {P}_{n-2}$ formulation (Maday & Patera 1989). A third-order backward difference is used for the time integration of the equations. The viscous terms are evaluated implicitly while extrapolation is used for the nonlinear terms. Over integration is used for a consistent evaluation of the nonlinear terms and a relaxation term based on the high-pass filtered velocity field is used to stabilize the method (Negi 2017). The stabilization method is based on the approximate deconvolution model with relaxation term (ADM-RT) method used for the large-eddy simulation of transitional flows (Schlatter, Stolz & Kleiser 2004, 2006) and has been validated with channel flows and flow over wings in Negi et al. (2018). Moving boundaries are treated using the ALE formulation (Ho & Patera 1990, 1991) and the fluid and structural equations are coupled using the Green's function decomposition approach whereby the geometrical nonlinearity is evaluated explicitly via extrapolation, while the added-mass effects are treated implicitly (Fischer, Schmitt & Tomboulides 2017). In addition, we have also implemented a fully implicit fixed-point nonlinear iteration method (Küttler & Wall 2008) for coupling of the fluid and structural equations. Both methods give nearly identical results in all tested cases. For the nonlinear and linear stability results reported in this work, the semi-implicit method of Fischer et al. (2017) has been used. For the case of adjoint analysis, the fully implicit fixed-point iterations have been used. All quantities are non-dimensionalized using the fluid density $\rho$, free stream velocity $U_{\infty }$ and the cylinder/ellipse diameter $D$.

3.2. Oscillating cylinder at subcritical Reynolds numbers

For the first case, we investigate a two-dimensional circular cylinder free to oscillate in the cross-stream direction subject to the action of a spring-damper system. The Reynolds number of the flow based on the cylinder diameter is $Re=23.512$. The inlet of the computational domain is $25$ diameters upstream of the cylinder while the outflow boundary is $60$ diameters downstream of the cylinder. The lateral boundaries are 50 diameters away on either side. A uniform inflow Dirichlet boundary condition is applied on the inflow and the lateral boundaries while the stress-free boundary condition is applied on the outflow boundary. The computational domain is discretized using $2284$ spectral elements which are further discretized into $10\times 10$ GLL points for a ninth-order polynomial representation for the velocity. This amounts to a total of 228 400 degrees of freedom in the domain. The base flow for all cases was calculated by keeping the structure fixed at its initial position. At $Re=23.512$, the flow with a fixed cylinder is linearly stable. The convergence to steady state was accelerated by using BoostConv (Citro et al. 2017). Figure 2 shows the calculated base flow state (streamwise velocity). Sen, Mittal & Biswas (2009) report a linear empirical relation for the length of the reverse flow region behind the cylinder which predicts $L/D=1.146$ for the given Reynolds number. The length of the reverse flow region was found to be $L/D=1.150$ for the calculated steady flow, matching well with the results of Sen et al. (2009).

Figure 2. Streamwise velocity of the base flow state for a cylinder in cross-flow at $Re=23.512$.

We denote $\eta$ as the vertical position of the cylinder and the structural equation is modelled using a spring-mass-damper system for the vertical displacement of the cylinder. Following the earlier notation, we denote the fluid domain as $\varOmega ^{f}$, the structural domain as $\varOmega ^{s}$ and the FSI interface as $\varGamma$. In all our cases, the structural equation contains a second-order derivative in time. The order of the time derivative is reduced by introducing a variable substitution for the velocity $\varphi$ of the structure. The combined nonlinear equations for the FSI system can be written as

(3.1a)\begin{gather} \left.\frac{\partial U_{i}}{\partial t}\right\vert_{\boldsymbol{W^{g}}} + (U_{j} -W^{g}_{j}) \frac{\partial U_{i}}{\partial x_{j}} +\frac{\partial P}{\partial x_{i}} - \frac{1}{Re} \frac{\partial^{2} U_{i}}{\partial x_{j}\partial x_{j}} = 0\quad \text{in } \varOmega^{f}, \end{gather}

(3.1b)\begin{gather}\frac{\partial U_{i}}{\partial x_{i}} = 0\quad \text{in } \varOmega^{f}, \end{gather}

(3.1c)\begin{gather}\mathcal{M}\frac{\textrm{d} \varphi}{\textrm{d} t} + \mathcal{D}\varphi + \mathcal{K}\eta - \mathcal{F}_{2} = 0\quad \text{for } \varOmega^{s}, \end{gather}

(3.1d)\begin{gather}\frac{\partial \eta}{\partial t} - \varphi = 0\quad \text{for } \varOmega^{s}, \end{gather}

(3.1e)\begin{gather}U_{1} = 0\quad \text{on } \varGamma, \end{gather}

(3.1f)\begin{gather}U_{2} -\varphi = 0\quad \text{on } \varGamma, \end{gather}

(3.1g)\begin{gather}\mathcal{F}_{i} - \oint_{\varGamma} \left[p\delta_{ij} -\frac{1}{Re} \left(\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}} \right) \right]n_{j}\partial\varOmega = 0\quad \text{fluid forces on } \varGamma. \end{gather}

The parameters for the structural system are specified in table 1. The mass $\mathcal {M}$ corresponds to a density ratio of $7$ (solid to fluid), and the damping constant $\mathcal {D}$ is set to $0.754\,\%$ of the critical damping. Here $\omega _{n} = \sqrt {\mathcal {K}/\mathcal {M}}$ represents the undamped natural frequency of the spring-mass system, $\lambda _{s}=-\omega _{n}(\zeta \pm \textrm {i}\sqrt {1-\zeta ^{2}})$ represents the damped structural eigenvalue of the spring-mass-damper system, with $\zeta =\mathcal {D}/(2\sqrt {\mathcal {K}\mathcal {M} })$ being the damping ratio.

Table 1. Structural parameters for a two-dimensional cylinder oscillating in the cross-stream direction.

Following the Taylor expansion based triple decomposition of the total fluid velocity and pressure, the linearized system of equations governing the fluid and structural perturbations $(\boldsymbol {u}',p',\eta ',\varphi ')$ can be written as

(3.2a)\begin{gather} \frac{\partial u'_{i}}{\partial t} + U^{0}_{j}\frac{\partial u'_{i}} {\partial x^{0}_{j}} + u'_{j} \frac{\partial U^{0}_{i}}{\partial x^{0}_{j}} +\frac{\partial p'}{\partial x^{0}_{i}} - \frac{1}{Re}\frac{\partial^{2} u'_{i}}{\partial x^{0}_{j}\partial x^{0}_{j}} = 0\quad \text{in } \varOmega^{f} , \end{gather}

(3.2b)\begin{gather}\frac{\partial u'_{i}}{\partial x^{0}_{i}} = 0\quad \text{in } \varOmega^{f}, \end{gather}

(3.2c)\begin{gather}u'_{i} = 0,\quad \text{on } \partial\varOmega_{v} , \end{gather}

(3.2d)\begin{gather}\sigma'_{ij}n^{0}_{j} = 0\quad \text{on } \partial\varOmega_{o} , \end{gather}

(3.2e)\begin{gather}\mathcal{M}\frac{\textrm{d} \varphi'}{\textrm{d} t} + \mathcal{D}\varphi' + \mathcal{K}\eta' - \mathcal{F}'_{2} = 0\quad \text{for } \varOmega^{s} , \end{gather}

(3.2f)\begin{gather}\frac{\partial \eta'}{\partial t} - \varphi' = 0\quad \text{for } \varOmega^{s} , \end{gather}

(3.2g)\begin{gather}u'_{1} + \eta'\frac{\partial U^{0}_{1}}{\partial x^{0}_{2}} = 0\quad \text{on } \varGamma, \end{gather}

(3.2h)\begin{gather}u'_{2} + \eta'\frac{\partial U^{0}_{2}}{\partial x^{0}_{2}} - \varphi' = 0\quad \text{on } \varGamma, \end{gather}

(3.2i)\begin{gather}\mathcal{F}'_{i} - \oint_{\varGamma^{0}}\left[p'\delta_{ij} -\frac{1}{Re}\left(\frac{\partial u'_{i}}{\partial x^{0}_{j}} + \frac{\partial u'_{j}}{\partial x^{0}_{i}} \right)\right]n^{0}_{j} \partial\varOmega = 0\quad \text{fluid forces on } \varGamma. \end{gather}

Note that even though the cylinder is free to move only in the vertical direction, the streamwise perturbation velocity $u'_{1}$ is non-zero at the cylinder surface due to the Taylor expansion term of the base flow. Additionally, when evaluating the fluid forces in (3.2i), the direction of the normal $n^{0}_{j}$ is based on the equilibrium position of the FSI interface.

Before calculating the spectra of the linearized problem, we compare the evolution of the linear and nonlinear flow cases when the base flow is perturbed by identical small-amplitude perturbations. If the perturbation amplitude is small then nonlinear effects will be negligible and one can expect the linear and nonlinear evolution to be the same. The solutions would eventually diverge due to amplitude saturation in the nonlinear case. Performing this comparison has another advantage in this particular case. Given the low Reynolds number of the case, we expect the dynamics of be governed by a single least damped global mode. Therefore, by tracking the growth in perturbation amplitude through the linear regime, one can determine both the frequency and the growth rate of the unstable mode from the nonlinear simulations. This allows us to validate both the frequency and the growth rate obtained from the linearized equations.

The base flow is disturbed by pseudo-random perturbations of order $O(10^{-6})$ and the flow evolution is tracked. Figure 3(a) shows the variation of $\eta$ with time during the initial stages of the evolution. Both the linear and nonlinear results fall on top of each other providing the first evidence of the validity of the derived linear equations. Note that the nonlinear simulations have been performed using the ALE framework including the mesh movement, while there is no mesh motion in the linear equations which have been derived to be independent of the grid velocity $\boldsymbol {W^{g}}$. The time evolution of $\eta$ clearly indicates a single growing mode in the flow. By tracking the peak amplitudes of the oscillations, denoted as $\eta ^{pks}$, one can determine the growth rate and the frequency of the unstable mode. Figure 3(b) shows the time evolution of $\eta ^{pks}$ in a semi-log plot. After an initial transient phase both the growth rate and angular frequency of the disturbances stabilize to a constant value and the $\eta ^{pks}$ plot traces a straight line in the semi-log plot, signifying exponential growth in time.

Figure 3. Comparison of linear and nonlinear evolution of cylinder position $\eta$ for identical small-amplitude disturbance. (a) The time evolution of $\eta$ in the first few cycles of the oscillation. (b) The evolution of peak amplitudes of the oscillation in a semi-log scale.

Finally, we introduce the ansatz $(\boldsymbol {u}',p',\eta ',\varphi ') = (\hat {\boldsymbol {u}},\hat {p},\hat {\eta },\hat {\varphi })\,\textrm {e}^{\lambda t}$ to reduce the linear equations into an eigenvalue problem for the angular frequency $\lambda$. The system is solved by estimating the eigenvalues of the time-stepping operator. The method was first used by Eriksson & Rizzi (1985) and has been used in several previous works by Barkley, Gomes & Henderson (2002), Bagheri et al. (2009b) and Brynjell-Rahkola et al. (2017) etc. The eigenvalue estimation is done using the implicitly restarted Arnoldi method (Sorensen 1992) implemented in the open source software package ARPACK (Lehoucq, Sorensen & Yang 1998). Figure 4 shows the one-sided eigenspectra obtained with a single unstable mode. Here $\lambda _{r}$ represents the real part of the eigenvalue (growth rate) and $\lambda _{i}$ represents the imaginary part of the eigenvalue (angular frequency). In table 2 we report the comparison of estimates obtained through the nonlinear simulations, linear simulations and the Arnoldi method. For both the linear and nonlinear simulations, the initial transient period of 10 oscillation cycles was discarded and the estimates are evaluated from the time period when both the growth rate and frequency have stabilized. All three methods have a very good agreement with each other, with the relative difference in the growth rate being less than $0.1\,\%$. The streamwise velocity component of the eigenvector for the unstable frequency is shown in figure 5.

Figure 4. One-sided eigenspectrum for a cylinder in cross-flow at $Re=23.512$.

Figure 5. (a) Real and (b) imaginary parts of the streamwise velocity component of the eigenvector corresponding to the unstable eigenvalue for $Re=23.512$.

Table 2. Unstable eigenvalue estimates obtained from different methods for an oscillating cylinder at $Re=23.512$.

Thus, we find the flow is unstable for an oscillating cylinder at nearly half the critical Reynolds number for the stationary case. Cossu & Morino (2000) also found instability at the same Reynolds number for a cylinder oscillating in the cross-flow direction, albeit for slightly different structural parameters.

The undamped natural frequency of the structural system $\omega _{n}$ is varied parametrically to investigate the changes in instability while keeping the density ratio constant. The damping constant is also varied such that the damping is always equal to $0.754\,\%$ of the critical damping of the spring-mass-damper system. The one-sided spectra of the various cases is shown in figure 6. The instability of the system arises for a narrow range of $\omega _{n}$ with the peak of the instability centred around $\omega _{n}\approx 0.796$. The system rapidly becomes stable again as the structural frequency is varied away from the peak value. In the range where the system is unstable, the unstable frequency of the combined FSI system is close to the undamped frequency $\omega _{n}$. Much of the low frequency spectra remains unaffected by the variation of $\omega _{n}$. We note that some of the results are in contrast with some of the findings of Cossu & Morino (2000) who performed global stability of the oscillating cylinder at the same Reynolds number and density ratio. In their study, the authors find a low frequency unstable mode with eigenvalue $\lambda = 1.371\pm 0.194i$ for a structural eigenvalue of $\lambda _{s}=-0.01 \pm 1.326i$. The flow case marked with diamonds in figure 6 corresponds to the same case investigated by Cossu & Morino (2000). While we find the system to be unstable at the subcritical Reynolds number of $Re=23.512$, we do not find the instability for the same structural parameters. Unlike Cossu & Morino (2000), we also do not find the existence of a low frequency unstable mode within the investigated structural parameters. In all our investigated cases, the effect of the structure on the spectrum remains confined close to the angular frequency of the spring-mass-damper system. (Note that the non-dimensionalization of length in Cossu & Morino (2000) is with respect to the radius of the cylinder, while it is with respect to the diameter in the current study. Hence, the respective angular frequencies are doubled in the current study.)

Figure 6. One-sided eigenspectra for a cylinder in cross-flow at $Re=23.512$ for varying undamped frequency, $\omega _{n}$, of the spring-mass system.

We consider another subcritical case at $Re=40$ which has been investigated in Navrose & Mittal (2016). In this case the authors consider an oscillating cylinder in cross-flow, with a mass ratio of 10 and without any structural damping ($\mathcal {D}=0$). The results are reported for the variation of the reduced velocity $U^{*}$ which is defined as $U^{*}=U_{\infty }/(\,f_{n}D)$, where $f_{n}=\omega _{n}/(2{\rm \pi} )$ is the natural frequency of the spring-mass system. We investigate the case with $U^{*}=8$, again comparing evolution of small-amplitude perturbations in the linear and the nonlinear cases. The initial cycles and the evolution of peak amplitude is shown in figures 7(a) and 7(b). The estimated eigenvalues from the linear, nonlinear and Arnoldi method are reported in table 3. The estimates match to within $0.1\,\%$ accuracy. For the nonlinear simulations, the amplitude of oscillations saturates at $y/D=0.4$. This compares well with the saturation amplitude reported in Navrose & Mittal (2016) for this case.

Figure 7. Comparison of linear and nonlinear evolution of cylinder position $\eta$ at $Re=40$ and $U^{*}=8$. (a) Time evolution of $\eta$ in the first few cycles of the oscillation. (b) Evolution of peak amplitudes in a semi-log scale.

Table 3. Unstable eigenvalue estimates obtained from different methods for an oscillating cylinder at $Re=40$ and $U^{*}=8$.

The variation of stability characteristics with varying reduced velocities is also investigated through the evaluation of the eigenvalue spectra for several different cases. The variation of spectra for different parameters is shown in figure 8(a). Similar to the trend seen for the $Re=23.512$ case, the flow exhibits unstable eigenvalues within a narrow band of frequencies. The variation of growth rate with the reduced velocities is shown in figure 8(b). The growth rates for $U^{*}=6$ and $U^{*}=10$ lie just above the stability threshold. This compares very well with the results of Navrose & Mittal (2016) who report the unstable range as $5.9<U^{*}<10.1$ as the parameter range of instability for the oscillating cylinder at $Re=40$. The peak growth rates occur at $U^{*}\approx 8.0$ in Navrose & Mittal (2016), which is also the case in the current work.

Figure 8. (a) One-sided spectra for an oscillating cylinder at $Re=40$ for varying reduced velocities $U^{*}$. (b) Variation of the growth rate with reduced velocity $U^{*}$.

3.3. Confined oscillating cylinder

In the previous two cases we simulated physical scenarios with effectively unbounded domains. In this subsection we consider the case of an FSI instability in a bounded domain. We replicate the physical set-up studied in Semin et al. (2012), where a cylinder is free to oscillate in the vertical direction and is confined between two parallel walls. The channel height is denoted at $h$ and the ratio of the cylinder diameter to channel height is $D/h=0.66$. The computational domain is set up to match the simulations of Semin et al. (2012), with the inlet located at $-5h$ and the outflow boundary located at $7h$. A parabolic profile is imposed on the inlet for the streamwise velocity and the mean streamwise velocity ($\bar {U}$) at the inlet is used for the normalization of the velocity scales. The length scale is normalized by the channel height $h$, and the Reynolds number is defined using these two scales as $Re=\bar {U}h/\nu$. The solid to fluid density ratio is set to $1.19$, which corresponds to the density ratio studied by Semin et al. (2012) both numerically and experimentally. The cylinder is free to oscillate in the vertical direction without any external restoring forces or damping. Thus, the only forces acting on the cylinder are the fluid forces. We investigate the loss of stability in the FSI system for three different Reynolds numbers of $Re=25, 22$ and $20$.

Following the usual procedure outlined earlier, the base flow for the three different Reynolds numbers is obtained for a stationary cylinder. Small-amplitude perturbations are then seeded into the system and the linear and nonlinear evolution of the system is tracked in time. The base flow state for $Re=25$ is shown in figure 9. Figure 10(a) shows the evolution of the peak amplitudes of the oscillation. Again, evident from the graph is the fact that exponential growth is well captured by the linear formulation, and both linear and nonlinear regimes undergo the same amplification through several orders of magnitude. Figure 10(b) shows the eigenvalue spectra for the three cases. A comparison of the linear, nonlinear and Arnoldi approximations of the unstable eigenvalue is given in table 4 and the unstable eigenvector is shown in figure 11. The saturated limit-cycle amplitudes saturate at $\eta _{max}=0.090, 0.076$ and $0.062$ for $Re=25, 22$ and $20$, respectively. The corresponding Strouhal numbers, defined as $Sr=fh/\bar {U}$, with $f$ being the saturated oscillation frequency, are calculated as $Sr=1.13, 1.17$ and $1.22$. The values correspond well with the results of Semin et al. (2012) for the given density ratio.

Figure 9. Streamwise velocity of the base flow state for a confined cylinder at $Re=25$.

Figure 10. (a) Comparison of linear and nonlinear evolution of peak amplitudes vertical displacement $\eta$ for a confined cylinder at three different Reynolds numbers. (b) One-sided eigenvalue spectra for the three different Reynolds number cases.

Figure 11. Real (a) and imaginary (b) parts of the streamwise velocity component of the unstable eigenvector for a confined cylinder at $Re=25$.

Table 4. Unstable eigenvalue estimates for a confined cylinder at different Reynolds numbers.

3.4. Asymmetric flow case of a rotated ellipse

The previous subsections investigated the flow around an oscillating circular cylinder through the linear stability analysis. The base flow for these cases exhibits symmetry about the horizontal axis passing through the origin. In order to test the linear formulation a case where no such symmetries arise we consider the case of a rotated ellipse in an open flow. The ellipse geometry is generated with the minor axis length of $a=0.25$ aligned with the streamwise direction, the major axis length of $b=0.5$ aligned in the cross-flow direction and the centre of the ellipse located at the origin of the coordinate system $(0,0)$. The ellipse is then rotated by an angle of $30^{\circ }$ clockwise. The stabilized base flow around the rotated ellipse is calculated at $Re=50$, where the diameter along the major axis is used as the length scale for the Reynolds number. The streamwise velocity for the stabilized base flow is shown in figure 12, which clearly shows the lack of symmetry of the base flow close to the ellipse.

Figure 12. Streamwise velocity of the base flow state for a rotated ellipse at $Re=50$.

We consider two cases of FSI: an ellipse free to oscillate in the cross-flow direction and a case with the ellipse free to rotate about the out-of-plane axis passing through its centre. In both cases the ellipse is considered to be held stationary due to a constant external force which balances the fluid forces at equilibrium. The density ratio is set to 10 for the ellipse and no spring or damping forces are considered. Thus, the ellipse motion is governed entirely due to fluid forces and the inertia of the ellipse. We expect any modelling errors in fluid forces to show up strongly for such a case. Proceeding in the usual manner, the comparison of small-amplitude linear and nonlinear evolution for the oscillating case is shown in figures 13(a) and 13(b), while the comparison for the rotational case is shown in figures 13(c) and 13(d). The eigenvalue estimates for both cases are shown in table 5. A good agreement is found for all the cases considered. We note that for the case with vertical oscillations, a zero frequency unstable mode also exists. This was deduced from the time evolution of the simulations since the positive and negative peaks showed a marginally different growth rate. The Arnoldi procedure also showed the existence of both an oscillatory and zero frequency unstable mode. A nonlinear least squares procedure was then used to estimate the growth rate and the unstable frequencies. The one-sided spectra for the two cases is shown in figure 14 and the streamwise components of the eigenvector associated with the (oscillatory) unstable eigenvalues are shown in figure 15. We mention that we have evaluated the ‘added stiffness’ terms of Fanion et al. (2000) and Fernández & Le Tallec (2003a,b) for all of the tested cases and we find that these terms always remain several orders of magnitude smaller than the forces arising due to the fluid perturbations. This further confirms the earlier mathematical result (2.24) that the ‘added stiffness’ terms represent a higher-order correction and do not play a role in the linearized dynamics.

Figure 13. Comparison of linear and nonlinear evolution of the position of a rotated ellipse at $Re=50$. The top two panels represent the vertical translation cases while the bottom two represent the rotational cases. (a) Time evolution of vertical position $\eta$ in the first few cycles of the oscillation. (b) Evolution of peak amplitudes in a semi-log scale for vertical oscillations. (c) Time evolution of rotational angle $\eta$. (d) Evolution of peak amplitudes of the rotational angle in a semi-log scale.

Figure 14. One-sided spectra for a rotated ellipse at $Re=50$. (a) Spectrum for an ellipse free to oscillate in the vertical direction and (b) spectrum for the ellipse free to rotate about its centre.

Figure 15. Streamwise velocity component of the eigenvector for (a,b) vertical oscillation and (c,d) rotational FSI cases. The left panels show the real part of the eigenvector and the right panels show the imaginary part of the eigenvector corresponding to the unstable eigenvalue.

Table 5. Unstable eigenvalue estimates for vertical oscillation and rotational cases for a rotated ellipse at $Re=50$, obtained with three different methods.

3.5. Spontaneous symmetry breaking

We use the derived linear formulation to investigate the case of a circular cylinder with an attached splitter plate set at zero incidence to the oncoming flow. The structural equations take the same form as in (3.1c), with $\eta$ now representing the deviation of the rotational angle from the equilibrium position. The terms $\mathcal {M}, \mathcal {D}, \mathcal {K}$ and $\mathcal {F}$ represent the moment of inertia, rotational damping, rotational stiffness and the out-of-plane moment acting on the body, respectively. The cylinder is free to rotate about its centre and the rotational stiffness $\mathcal {K}$ and damping $\mathcal {D}$ are both set to zero. Thus, the cylinder rotates only due to the action of the fluid forces. The particular set-up, along with a variant where the splitter plate is flexible, has been a subject of investigation by several previous authors (Xu, Sen & Gad-el Hak 1990, 1993; Bagheri, Mazzino & Bottaro 2012; Lācis et al. 2014). For certain lengths of the cylinder-splitter body, the system exhibits an interesting dynamic where the body spontaneously breaks symmetry and the splitter plate settles at a non-zero angle to the oncoming flow. The breaking of symmetry leads to the generation of lift force which could play a role in locomotion through passive mechanisms (Bagheri et al. 2012). The symmetry breaking phenomenon is known to occur for both subcritical and supercritical Reynolds numbers (Lācis et al. 2014) and in both two- and three-dimensional configurations (Lācis et al. 2017). We investigate the phenomenon through our linear FSI framework. According to the results of Lācis et al. (2014), symmetry breaking is expected to occur for splitter-plate lengths of less than $2D$. Accordingly, we set the splitter-plate length to be $1D$ and a thickness of $0.02D$. Following Lācis et al. (2014) we use a solid to fluid density ratio of $1.001$ for $Re=45$ and $1.01$ for $Re=156$. Figures 16(a) and 16(b) show the base flow states for the diameter based Reynolds numbers of $Re=45$ (subcritical) and $Re=156$ (supercritical).

Figure 16. Base flow state for a cylinder with splitter plate at (a) $Re=45.0$ and (b) $Re=156$. Note that the splitter plate is not discernible in the visualization due to its very small thickness ($0.02D$).

Figure 17 shows the spectra for the two different Reynolds numbers. Both cases have an eigenvalue lying on the positive $y$-axis ($\lambda =0.10+0i$ for $Re=45$ and $\lambda =0.19+0i$ for $Re=156$). This represents the symmetry breaking eigenmode of the system since it does not oscillate about a zero mean but rather leads to a monotonic growth in the rotational angle. For the higher Reynolds number, several other modes of instability exist in the flow which represent the von-Kármán modes of the flow. From the perspective of linear analysis, these modes simply oscillate about the mean angle with growing oscillation amplitudes. However, the zero frequency mode leads to a monotonic rise in the angle and, thus, causes the symmetry breaking effect. For $Re=45$, we also simulate the linear and nonlinear cases after adding a random small-amplitude perturbations to both flow cases. Figure 18 shows the time evolution of the angle $\eta$ of the cylinder-splitter body. Both the simulations undergo the same exponential growth of about five orders of magnitude before the nonlinear case saturates, while the linear case continues its exponential growth. The saturation of the nonlinear simulations occur at $\eta \approx 15.5^{\circ }$, which is the same turn angle reported by Lācis et al. (2014) for a splitter plate of length $1D$. Thus, the linear FSI framework predicts the onset of symmetry breaking instability. Of course for the system to finally exhibit symmetry breaking a nonlinear mechanism is required since the flow must equilibrate at the new position. However, the onset can be traced to the zero frequency unstable mode.

Figure 17. One-sided spectra for the cases exhibiting spontaneous symmetry breaking at $Re=45$ and $Re=156$ for a rotatable cylinder with splitter plate.

Figure 18. Time evolution of the angle for a rotatable cylinder with the splitter plate at $Re=45$.

As noted earlier, the work of Pfister et al. (2019) also considered the problem of linear instability in FSI problems and considered more general structural models including nonlinear structural models and finite aspect ratio geometries. The approach towards linearization of the fluid quantities however has been different from the one taken in the current study. It is therefore interesting to contrast the linearization approaches of the two studies. In particular, both studies start with the ALE form of the full Navier–Stokes and clearly realize the importance of treating the motion of material points consistently. The point of departure between the two studies has been on the approach towards the treatment of the material point motion in the fluid domain. Pfister et al. (2019) take the approach of incorporating the changes into the equations itself. On the other hand, in the current work the linearized operator at the deformed material points is further approximated by way of a Taylor series expansion and geometric linearization. The final form of the linearized Navier–Stokes obtained is the obvious conclusion of these different approaches. Pfister et al. (2019) obtain a modified linear operator on a moving grid which incorporates the effects of domain deformation, while the current work obtains the modified linear operator defined on the original stationary grid (the modification arising only at the boundaries). It appears the crucial aspect in the FSI linearization is the consistent treatment of material point motion. Finally, we make a note that the PhD thesis of Pfister (2019) reports some cases with more complex structural models where the added stiffness terms appear to have a non-trivial effect on the instability characteristics of the FSI problems.

4. Structural sensitivity of the eigenvalue

As noted in figures 6 and 8(a), the unstable eigenvalue of the coupled FSI system changes as the structural parameters are varied, while the base flow is held constant. The variation occurs both for the growth rate as well as the frequency of the eigenvalue. One might expect the sensitivity of the eigenvalue to structural perturbations to change as well. We investigate the eigenvalue sensitivity to structural perturbations for different values of the structural parameters for a spring-mounted two-dimensional circular cylinder, which is free to oscillate in the cross-stream direction. The linear FSI problem is defined by (3.2a)–(3.2i). In the following sections we assume all derivatives are evaluated in the reference configuration and we drop the superscript $^{0}$ from the derivative terms. In addition, uppercase letters denote base flow quantities and lowercase letters denote perturbation quantities, and the superscripts $^{0}$ and $'$ are dropped from the base flow and perturbation quantities.

The eigenvalue sensitivity is typically studied through the use of adjoint equations (Giannetti & Luchini 2007; Luchini & Bottaro 2014). For any linear operator $\mathcal {L}$, the adjoint operator $\mathcal {L}^{\dagger }$ is defined such that it satisfies the Lagrange identity

(4.1)\begin{equation} \langle \mathcal{L}^{\dagger}p,q\rangle = \langle p,\mathcal{L}q\rangle \end{equation}

for any arbitrary vectors $q$ and $p$ in the domain of $\mathcal {L}$ and $\mathcal {L}^{\dagger }$, respectively. The symbol $\langle \cdot ,\cdot \rangle$ denotes the inner product under which the above identity holds. In the current context the operator $\mathcal {L}$ represents the linearized Navier–Stokes equations for FSI, also referred to as the direct problem, $\mathcal {L}^{\dagger }$ is the corresponding adjoint operator with the definition of the inner product as the integral over the domain $\varOmega$ and time horizon $T$:

(4.2)\begin{equation} \langle p,q\rangle = \int_{T}\int_{\varOmega} p^{H}q\,\textrm{d}\varOmega\,\textrm{d} t. \end{equation}

Defining the vector $\boldsymbol {q}=(\boldsymbol {u},p,\eta ,\varphi )$, which lies in the domain of $\mathcal {L}$, and an adjoint vector $\boldsymbol {q}^{\dagger }=(\boldsymbol {u}^{\dagger },p^{\dagger },\eta ^{\dagger },\varphi ^{\dagger })$ which lies in the domain of $\mathcal {L^{\dagger }}$, the adjoint equations may be found by taking the inner product between the adjoint variables and the direct equations. The use of integration by parts, divergence theorem and a judicious manipulation of terms leads to a set of expressions for the adjoint equations. The full derivation of the equations along with the boundary conditions may be found in appendix B. The final form of the adjoint equations is expressed as

(4.3a)\begin{gather} -\frac{\partial u^{\dagger}_{i}}{\partial t} + u^{\dagger}_{j} \frac{\partial U_{j}}{\partial x_{i}}-U_{j}\frac{\partial u^{\dagger}_{i}}{\partial x_{j}} + \frac{\partial p^{\dagger}}{\partial x_{i}} - \frac{1}{Re}\frac{\partial }{\partial x_{j}} \left(\frac{\partial u^{\dagger}_{i}}{\partial x_{j}} + \frac{\partial u^{\dagger}_{j}}{\partial x_{i}}\right) =0\quad \text{in } \varOmega^{f}, \end{gather}

(4.3b)\begin{gather}\frac{\partial u^{\dagger}_{i}}{\partial x_{i}} = 0\quad \text{in } \varOmega^{f}, \end{gather}

(4.3c)\begin{gather}-\mathcal{M}\frac{\textrm{d} \varphi^{\dagger}}{\textrm{d} t} +\mathcal{D}\varphi^{\dagger} -\eta^{\dagger} + \oint_{\partial\varOmega} \sigma^{\dagger}_{2j}n_{j} \,\textrm{d}\varOmega = 0\quad \text{in } \varOmega^{s}, \end{gather}

(4.3d)\begin{gather}-\frac{\textrm{d}\eta^{\dagger}}{\textrm{d} t} + \mathcal{K}\varphi^{\dagger} - \oint_{\partial\varOmega} \left(\frac{\partial U_{k}}{\partial x_{2}}\right) \sigma^{\dagger}_{kj}n_{j} \,\textrm{d}\varOmega = 0\quad \text{in } \varOmega^{s}, \end{gather}

(4.3e)\begin{gather}u^{\dagger}_{1} = 0\quad \text{on } \varGamma, \end{gather}

(4.3f)\begin{gather}u^{\dagger}_{2} - \varphi^{\dagger}_{2} = 0\quad \text{on } \varGamma, \end{gather}

(4.3g)\begin{gather}\sigma_{ij}^{\dagger} - \left[-p^{\dagger}\delta_{ij} + \frac{1}{Re}\left(\frac{\partial u^{\dagger}_{i}}{\partial x_{j}} + \frac{\partial u^{\dagger}_{j}}{\partial x_{i}} \right)\right] = 0\quad \text{adjoint forces on } \varGamma, \end{gather}

(4.3h)\begin{gather}u^{\dagger}_{i} = 0\quad \text{on } \partial\varOmega_{v}, \end{gather}

(4.3i)\begin{gather}(\sigma^{\dagger}_{ij}+u_{i}^{\dagger}U_{j})n_{j} = 0\quad \text{on } \partial\varOmega_{o}. \end{gather}

Here $\sigma _{ij}^{\dagger }$ is referred to as the adjoint fluid stresses acting on the cylinder. Unlike in the direct equations, the adjoint fluid velocities are spatially homogeneous on the FSI interface $\varGamma$. The sign of the diffusive term indicates the equations evolve an adjoint field backward in time. A variable substitution $\tau = - t$ can be made to return the equations to forward (in $\tau$) propagation.

A validation of the adjoint equations can be found by comparing the spectra of the direct and adjoint problems. The eigenvalues of the adjoint problem are the complex conjugate of the eigenvalues of the direct problem. For real matrices, the spectra for both problems is the same. Figure 19 shows the direct and adjoint spectra for a spring-mounted two-dimensional cylinder at a diameter based Reynolds number of $Re=50$, mass ratio of $10$, corresponding to $\mathcal {M}=7.854$, $\mathcal {K}=1.1537$ which corresponds to the natural frequency of the spring-mass system of $\omega _{n}=0.3833$. The damping is set to zero ($\mathcal {D}=0$). As observed in the figure, the spectra for both the direct and adjoint problems have a very good agreement with each other.

Figure 19. Comparison of the eigenspectra for the direct and adjoint problems for a density ratio of $10$ and natural frequency $\omega _{n}=0.3833$.

For a general linear eigenvalue problem $\mathcal {L}q = \lambda q$, a first-order approximation to the perturbed eigenvalue problem $(\mathcal {L} +{\rm \Delta} \mathcal {L})(q + {\rm \Delta} q) = (\lambda + {\rm \Delta} \lambda )(q +{\rm \Delta} q)$ can be obtained by making use of the corresponding adjoint eigenvector vector field $q^{\dagger }$ (Giannetti & Luchini 2007), resulting in an expression for the eigenvalue perturbation as

(4.4)\begin{gather} {\rm \Delta}\lambda = \frac{\langle q^{\dagger},{\rm \Delta}\mathcal{L}q\rangle}{\langle q^{\dagger},q\rangle}, \end{gather}

(4.5)\begin{gather}(\mathcal{L}^{\dagger} - \lambda^{*}\mathcal{I})q^{\dagger} = 0, \end{gather}

where $^{*}$ implies complex conjugation and $q^{\dagger }$ represents the right eigenvector of the adjoint operator $\mathcal {L}^{\dagger }$ with the eigenvalue of $\lambda ^{*}$. To estimate the drift in eigenvalue, knowledge of the specific operator perturbation ${\rm \Delta} \mathcal {L}$ is required. However, as shown by Giannetti & Luchini (2007), a spatial sensitivity map of the eigenvalue perturbation can be found by assuming the perturbation operator to be of the form of a spatially localized feedback, with the feedback being proportional to the local values of the field variables (velocities). For such a case, we may estimate the eigenvalue drift along with its upper bound for each spatially localized operator perturbation as

(4.6)\begin{gather} {\rm \Delta}\lambda = \frac{\langle q^{\dagger},\delta(x-x_{0},y-y_{0})C_{0} q\rangle}{\langle q^{\dagger},q\rangle}, \end{gather}

(4.7)\begin{gather}|{\rm \Delta}\lambda| \le ||C_{0}||\boldsymbol{\varTheta}(x_{0},y_{0}), \end{gather}

(4.8)\begin{gather}\boldsymbol{\varTheta}(x,y) = \frac{||q^{\dagger}(x,y)||\cdot ||q(x,y)||}{||\langle q^{\dagger},q\rangle||}, \end{gather}

where $C_{0}$ is a matrix that defines the feedback due to the localized operator perturbation and the local coupling between different velocity components. The quantity $\boldsymbol {\varTheta }(x,y)$ gives an indication of regions where the localized feedback will produce a large drift in the eigenvalue, thus representing the sensitivity of the eigenvalue to structural perturbations.

The structural sensitivity of a spring-mounted cylinder is investigated for varying natural frequencies $\omega _{n}$ of the spring-mass system. In all subsequent cases the structural damping is set to zero and the density ratio is set to $10$. The natural frequency of the spring-mass system is varied and the sensitivity map $\boldsymbol {\varTheta }$ of the least stable eigenvalue is evaluated. Table 6 shows the variation of the natural frequency $\omega _{n}$ and the corresponding least stable eigenvalue. The table also shows the unstable eigenvalue for a stationary cylinder at $Re=50$. Figure 20 shows the corresponding changes in structural sensitivity of the least stable eigenvalue. When the structural frequency is much lower than the unstable frequency (figure 20b), the sensitivity map resembles the sensitivity map of the stationary cylinder with small changes close to the cylinder. A small region of sensitivity is generated close to the cylinder, which is located symmetrically with respect to the horizontal axis. This region grows in intensity as $\omega _{n}$ approaches closer to the unstable frequency (20c). Close to resonance (20d), the dominant region shifts to the cylinder, lying symmetrically on the top and bottom. The near wake region is still sensitive to structural perturbations, however, it is no longer the dominant region. As $\omega _{n}$ is increased to be much larger than the unstable frequency, the sensitivity map resembles the stationary case again. In all cases with FSI, the sensitivity is non-zero at the cylinder surface. However, the values remain small if $\omega _{n}$ is far from resonance. Close to resonance, small regions on the cylinder surface have high sensitivities. A close-up of the sensitivity map for $\omega _{n}=0.7665$ is shown in figure 21, which shows non-zero values of $\boldsymbol {\varTheta }$ at the cylinder surface. The non-zero sensitivity on the cylinder surface has implications for control of vortex induced vibrations. In particular, it is found that control techniques for elimination of vortex streets which apply to stationary cylinders do not always work for vibrating cylinders (Dong, Triantafyllou & Karniadakis 2008).

Figure 20. Comparison of the structural sensitivity map for the least stable eigenvalue for different values of the natural frequency $\omega _{n}$. The top panel (a) is the structural sensitivity for the stationary cylinder case. The bottom four panels (b–e) represent the cases with FSI. (a) Stationary cylinder; (b) $\omega _{n}=0.0767$; (c) $\omega _{n}=0.3833$; (d) $\omega _{n}=0.7665$ and (e) $\omega _{n}=3.8327$.

Figure 21. Close-up of the sensitivity map for $\omega _{n}=0.7665$.

Table 6. Unstable eigenvalue corresponding to the structural natural frequency $\omega _{n}$ for the oscillating cylinder cases.

5. Conclusions

An Eulerian formulation for the linear stability analysis of rigid-body motion FSI problems is rigorously derived and validated. The final form of the linear equations is evaluated on the equilibrium grid on which the base state is defined and the first-order effects of moving interfaces are captured by the modification of the interface boundary conditions. The ‘added stiffness’ terms that arise in the formulation of Fanion et al. (2000) and Fernández & Le Tallec (2003a,b) are shown to vanish to a first-order approximation and play no role in the linearized perturbation dynamics. The formulation is validated for rigid-body motion FSI problems by comparing the evolution of the linear equations with the nonlinear system when both systems are perturbed with identical small-amplitude disturbances. This has been performed for several different cases for both rotation and translation motion, in symmetric as well as asymmetric flow cases, and bounded and unbounded flows. The linear and nonlinear equations evolve in a near identical manner and the extracted frequency and growth rates for the two cases match to within $0.1\,\%$ relative difference.

The FSI linear framework is used to analyse the case of symmetry breaking for a rotating cylinder with a rigid splitter plate. The linear stability analysis predicts the existence of a zero frequency unstable mode. This mode is identified as the cause of symmetry breaking since it causes the system to monotonically deviate from the equilibrium position and thus causing the onset of the symmetry breaking effect. The zero frequency unstable mode can be found for both the subcritical and supercritical Reynolds numbers.

Finally, the eigenvalue sensitivity to structural perturbations is investigated using the adjoint equations for FSI, for a two-dimensional cylinder oscillating in the cross-stream direction at $Re=50$, for varying natural frequencies of the spring-mass system. When the structural frequencies are far from the unstable frequency, the sensitivity is hardly affected. However, close to resonance, the sensitivity map changes completely and the dominant region of sensitivity lies close to the cylinder, located symmetrically above and below. In all cases of FSI, the sensitivity at the surface is non-zero, however, it attains large values only close to the resonance condition between the fluid and the structure.