Global stability of rigid-body-motion fluid-structure-interaction problems

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. (2000) vanish at the FSI interface in a first-order approximation. 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 non-linear 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 fluid-structure-interaction 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.


Introduction
Fluid-structure interaction (FSI) studies span a vast and diverse range of applications -from natural phenomenon such as the fluttering of flags (Shelley and Zhang, 2011), phonation (Heil and Hazel, 2011), blood flow in arteries (Freund, 2014), path of rising bubbles (Ern et al., 2012) etc., to the more engineering applications of aircraft stability (Dowell and Hall, 2001), vortex induced vibrations (Williamson and Govardhan, 2004), compliant surfaces (Riley et al., 1988;Kumaran, 2003) etc. The phenomena that emerge out of an FSI problem often exhibit highly non-linear, dynamically rich and complex behavior with different flow regimes such as fluttering and tumbling of plates falling under gravity (Mittal et al., 2004), or the unsteady path of rising bubbles (Mougin and Magnaudet, 2001;Ern et al., 2012). Despite its non-linear 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 (1959Benjamin ( , 1960 and Landahl (1962). For boundary layers over Krammer-type compliant surfaces, the linear hydrodynamic instability study was performed by 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 and Morris (1990); Davies and Carpenter (1997) etc.
When the parallel flow assumption is no longer valid a global approach needs to be taken. From the perspective of global instabilities in fluid-structure-interaction problems, the first such study was reported by Cossu and Morino (2000) for the instability of a spring mounted cylinder. The authors considered a rigid 2D 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 ≈ 23. Mittal and Singh (2005) also report a similarly low critical Re for a cylinder oscillating in both streamwise and transverse directions. A non-inertial reference is also used by Navrose and Mittal (2016) to study the lock-in phenomenon of cylinders oscillating in crossstream through the linear stability analysis. More recently, Magnaudet and co-workers have approached the problem of rising and falling bodies through the perspective of linear instability. Fabre et al. (2011) developed a quasi-static model to describe the stability of heavy bodies falling in a viscous fluid. Assemat et al. (2012) perform 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 et al. (2014a) performed similar studies for three-dimensional disks and thin cylinders, Tchoufag et al. (2014b) apply the linear stability analysis to spheroidal bubbles and Cano-Lozano et al. (2016) for oblate bubbles.
In general the 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 his. 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 pseudovariable. Fanion et al. (2000) have derived a more general formulation starting from the ALE equations in the weak form, which has been used by Fernández and 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 and Le Tallec, 2003a,b). Goza et al. (2018) have used an immersed boundary method for the global stability analysis of inverted flag-flapping. In a very recent development, Pfister et al. (2019) also derive the 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.
In the current work we follow the methodology of Fanion et al. (2000) and Fernández and 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 obtainted 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 non-linear 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 non-linear simulations evolve nearly identically through several orders of magnitude of growth of the perturbations. The non-linear 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 section 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 section 2 we describe the problem in a general setting and derive the linearized equations for a fluid-structureinteraction problem for rigid-body motion. Numerical validation and results for the derived formulation are presented in section 3. In section 4 we introduce the adjoint problem for the linear FSI system of a cylinder oscillating in crossflow and show the changes in structural sensitivity of the unstable eigenvalue. Section 5 concludes the paper.

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 fluid-structureinteraction problem as illustrated in figure 1. The bounded region in white marked by fluid forces act. The Navier-Stokes equations, along with the incompressibility constraint govern the evolution of the fluid in Ω f . For a problem with moving interfaces, the Navier-Stokes is usually formulated in the Arbitrary-Lagrangian-Eulerian (ALE) formulation Patera, 1990, 1991) defined on moving material points as Here U and P represents the fluid velocity and pressure respectively, while W g represents the velocity of the material points. We refer to the time derivative in this formulation (∂U i /∂t| W g ) as the ALE time derivative evaluated when following a material point with velocity W g . Note that W g is not uniquely defined. The only restriction on material point velocity is for the points on the interface Γ(t), where the fluid velocity is equal to the velocity of the interface points. Typically W g is defined by some function which smoothly extends the velocity at the interface to the rest of the fluid domain. The equations governing the structural motion depend on the type of modeling 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 Here M is a generalized inertia, D is a generalized damping and K is a generalized stiffness. F i represents the fluid forces acting on Γ(t). The equation represents a typical linear spring-mass-damper system for rigid-body motion. The structural degrees of freedom are represented by η i , whose definition depends on the type of modeling used for the structure. For example, for a cylinder free to oscillate in one direction subject to a spring-damper action (as in Cossu and Morino (2000)), η i represents the position of the center of mass in that direction, M is the cylinder mass, D is the damping coefficient and K is the spring stiffness constant. 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 x Γ as the (time-dependent) position of points on the interface Γ(t) and v Γ as their instantaneous velocity, we may write the boundary condition as The coupled problem is now defined by equations 1a and 1b in the fluid domain, equation 2 in the structural domain, coupled with the no-slip condition at the interface (equation 3

Steady-state
For the problem described in Section 2.1, consider a full system state given by (U 0 , P 0 , η 0 ) (referred to as the base flow state), defined on material points x 0 , which satisfies the time-independent incompressible Navier-Stokes and structural equations (along with the respective boundary conditions). i.e.
Where F 0 i are the fluid forces acting on the structure at equilibrium and Γ 0 is the equilibrium position of the FSI interface.

Linearization of the structure
We begin with the linearization of the structural equations. In defining the generalized structural equation 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 non-linear 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 η 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. η i = η 0 i + η i . Similarly, we decompose the fluid forces acting on the structure as F i = F 0 i + F i , where F i are the fluid forces acting on structure due to the (yet unknown) linearized fluid perturbations. Introducing the decomposition into equation 2 we obtain the equations for the structural perturbations Equation 5 will be completed if we can evaluate the expression for the term F i arising due to the linearized fluid perturbations. The two systems (fluid and structure) are coupled through the term F i and the velocity boundary conditions at the interface.

Taylor expansion based linearization of the fluid equations
To begin with, we define a linear invertible mapping between the material points in the reference (equilibrium) configuration and the perturbed points as Here I is the identity matrix and 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 For a small deformation, we may assume ||R || 1 and the matrix inverse may then be written as Retaining only the first order term of the expansion, one obtains the approximate inverse as We refer to equation 7 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 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 ∂/∂x 0 j represent the evaluation of derivatives in the reference configuration.
Next we define a perturbation field for the fluid components (u , p ). The perturbation field is also defined on the same equilibrium grid x 0 on which the base flow (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 The last term in equations 8a and 8b 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 P(x, y, z) = P 0 + p ξ + p .
Where, equations 9a and 9b are derived from equations 8a and 8b respectively after dropping the second-order terms, and using the compact notations u ξ i and p ξ 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 equation 9a amount to a triple decomposition of the total velocity field, where the three terms U 0 i , u ξ 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 Taking the time derivative of equation 9a one obtains the ALE time derivative of the total velocity along the trajectory of the material point motion 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 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 V I all represent terms that are second order in the perturbations quantities u , ∆x and R . Therefore, to a firstorder approximation these quantities can be dropped. This leaves the final divergence constraint as 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 The non-linear transport terms have already been dropped. We now substitute the ALE time derivative from equation 11 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: The terms in 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. On 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 equation 12 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 w:

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 Here n 0 j dS represents the the local surface vector in the equilibrium configuration. For compact notation we have denoted the stresses due the fluid at equilibrium as σ 0 i j , defined as The three integrals in equation 17 represent the momentum flux through the far-field (∂Ω f v ), outflow (∂Ω f o ) and FSI interface (Γ 0 ) respectively. In a similar manner the integral over the perturbed domain may be performed and without loss of generality we may assume that the outer boundaries ∂Ω f v and ∂Ω f o remain fixed, leading to an expression for the momentum flux at the perturbed surfaces as Where n j dS represents the change in the surface vector due to the motion of the boundary and we have dropped the second order terms for convection u ξ i u ξ j and the surface force σ ξ i j n j . Equation 19 represents a first-order approximation of the fluxes at the perturbed boundaries. The expression σ ξ i j is simply the first-order term of the Taylor expansion of the base flow-stresses, defined as Subtracting equation 17 from 19 one obtains The first term 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); Fernández and Le Tallec (2003a) and Fernández and Le Tallec (2003b) 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. One may now evaluate the total linearized forces arising from equation 15 on the perturbed FSI interface and obtain the linearized boundary conditions for the stress balance (equation 5) as where σ i j is the fluid stress due to the perturbation field (u , p ) defined as The velocity continuity boundary conditions on the moving interface Γ(t) simply become Thus the perturbation velocities must account for the Taylor expansion term of the base flow at the perturbed boundary. Fanion et al. (2000); Fernández and 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 (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 (1959Benjamin ( , 1960 and Landahl (1962) for parallel flow cases. We make a note that while in the current work the structural equations comprise of rigidbody motion equations, the linearization of the fluid equations did not invoke any such assumption. For a more general FSI problem, one would only need to focus on the linearization of the structural equations and the fluid linearization remains unchanged. In the next section several examples are considered for rigid-body motion in symmetric and asymmetric configurations to show the validity of the current approach.

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 spectralelement method (SEM) code (Fischer et al., 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 (GL) points for the representation of the pressure in a n − n−2 formulation (Maday and 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 non-linear terms. Over-integration is used for a consistent evaluation of the non-linear terms and a relaxation term based on the high-pass filtered velocity field (HPF-RT) is used to stabilize the method (Negi, 2017). The stabilization method is based on the ADM-RT method used for the large-eddy simulation (LES) of transitional flows (Schlatter et al., 2004(Schlatter et al., , 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 Patera, 1990, 1991) and the fluid and structural equations are coupled using the Green's function decomposition approach whereby the geometrical non-linearity is evaluated explicitly via extrapolation, while the added-mass effects are treated implicitly (Fischer et al., 2017). In addition we have also implemented a fully implicit fixed-point nonlinear iteration method (Küttler and Wall, 2008) for coupling of the fluid and structural equations. Both methods give nearly identical results in all tested cases. All quantities are non-dimensionalized using the fluid density ρ, free-stream velocity U ∞ , and the cylinder/ellipse diameter D.

Oscillating cylinder at sub-critical Reynolds Numbers
For the first case we investigate a 2D circular cylinder free to oscillate in the crossstream 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 × 10 GLL points for a 9 th order polynomial representation for the velocity. This amounts to a total of 228400 degrees of freedom in the domain. The base flow for all cases was calculated by keeping structure fixed at its initial position. At the 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).
We denote η as the vertical position of the cylinder and the structural equation is modeled using a spring-mass-damper system for the vertical displacement of the cylinder. Following the earlier notation we denote the fluid domain as Ω f , the structural domain as Ω s and the FSI interface as Γ. 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 ϕ of the structure. The combined M K D ω n λ s 5.4977 3.4802 6.597 × 10 −2 0.7956 −0.006 ± 0.7956i Table 1: Structural parameters for a 2D cylinder oscillating in the cross-stream direction.
non-linear equations for the FSI system can be written as The parameters for the structural system are specified in table 1. The mass M corresponds to a density ratio of 7 (solid to fluid), and the damping constant D is set to 0.754% of the critical damping. ω n = √ K/M, represents the undamped natural frequency of the spring-mass system, λ s = −ω n (ζ ± i 1 − ζ 2 ), represents the damped structural eigenvalue of the spring-mass-damper system, with ζ = D/(2 √ KM) being the damping ratio. 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 (u , p , η , ϕ ) can be written as 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 equation 26i, the direction of the normal n 0 j , is based on the equilibrium position of the FSI interface.
Before calculating the spectra of linearized problem we compare the evolution of the linear and non-linear flow cases when the base flow is perturbed by identical smallamplitude perturbations. If the perturbation amplitude is small then non-linear effects will be negligible and one can expect the linear and non-linear evolution to be the same. The solutions would eventually diverge due to amplitude saturation in the non-linear 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 non-linear 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 3a shows the variation of η with time during the initial stages of the evolution. Both the linear and non-linear results fall on top of each other providing the first evidence of the validity of the derived linear equations. Note that the non-linear 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 W g . The time evolution of η clearly indicates a single growing mode in the flow. By tracking the peak amplitudes of the oscillations, denoted as η pks , one can determine the the growth rate and the frequency of the unstable mode. Figure 3b shows the time evolution of η 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 η pks plot traces a straight line in the semi-log plot, signifying exponential growth in time.
Finally, we introduce the ansatz (u , p , η , ϕ ) = (û,p,η,φ)e λt to reduce the linear equations into an eigenvalue problem for the angular frequency λ. The system is solved by estimating the eigenvalues of the time-stepping operator. The method was first used by Eriksson and Rizzi (1985) and has been used in several previous works by Barkley et al. (2002); Bagheri et al. (2009);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 et al., 1998). Figure 4 shows the one-sided eigenspectra obtained with a single unstable mode. λ r represents the real part of the eigenvalue (growth rate) and λ i represents the imaginary part of the eigenvalue (angular frequency). In table 2 we report the comparison of estimates obtained through the non-linear simulations, linear simulations and the Arnoldi method. For both the linear and non-linear 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.
Thus we find the flow is unstable for an oscillating cylinder at nearly half the critical Reynolds number for the stationary case. Cossu and 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 ω n is varied parametrically to investigate the changes in instability while keeping the density ratio constant.

Case
Non-linear Linear Arnoldi Oscillation (Re = 23.512) 9.86 × 10 −3 ± 0.704i 9.85 × 10 −3 ± 0.704i 9.86 × 10 −3 ± 0.704i Table 2: Unstable eigenvalue estimates obtained from different methods for an oscillating cylinder at Re = 23.512.  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 ω n with the peak of the instability centered around ω n ≈ 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 ω n . Much of the low frequency spectra remains unaffected by the variation of ω n . We note that some of the results are in contrast with some of the findings of Cossu and 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 λ = 1.371 ± 0.194i for a structural eigenvalue of λ s = −0.01 ± 1.326i. The flow case marked with diamonds in figure 6 corresponds to the same case investigated by Cossu and Morino (2000). While we find the system to be unstable at the sub-critical Reynolds number of Re = 23.512, we do not find the instability for the same structural parameters. Unlike Cossu and 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 nondimensionalization of length in Cossu and 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). We consider another sub-critical case at Re = 40 which has been investigated in Navrose and 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 (D = 0) . The results are reported for the variation of the reduced velocity U * which is defined as U * = U ∞ /( f n D), where f n = ω n /(2π) is the natural frequency of the spring-mass
system. We investigate the case with U * = 8, again comparing evolution of smallamplitude perturbations in the linear and the non-linear cases. The initial cycles and the evolution of peak amplitude is shown in figure 7a and 7b. The estimated eigenvalues from the linear, non-linear and Arnoldi method are reported in table 3. The estimates match to within 0.1% accuracy. For the non-linear simulations the amplitude of oscillations saturates at y/D = 0.4. This compares well with the saturation amplitude reported in Navrose and Mittal (2016) for this case. 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 8a. 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 8b. The growth rates for U * = 6 and U * = 10 lie just above the stability threshold. This compares very well with the results of Navrose and 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 * ≈ 8.0 in Navrose and Mittal (2016), which is also the case in the current work.

Asymmetric flow case of a rotated ellipse
The previous subsection 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 crossflow 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 • 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 9 which clearly shows the lack of symmetry of the base flow close to the ellipse.
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 center. 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 modeling errors in fluid forces to show up strongly for such a case. Proceeding in the usual manner, the comparison of small-amplitude linear and non-linear evolution for the oscillating case is shown in figures 10a and 10b, while the comparison for the rotational case is shown in figures 10c and 10d. The eigenvalue estimates for both cases are shown in table 4. 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 non-linear 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 11 and the streamwise components of the eigenvector associated with the (oscillatory) unstable eigenvalues are shown in figure 12.
We mention that we have evaluated the "added stiffness" terms of Fanion et al. (2000); Fernández and Le Tallec (2003a) and Fernández and Le Tallec (2003b) 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 (equation 21) that the "added stiffness" terms represent a higher-order correction and do not play a role in the linearized dynamics.

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 equation 25c, with η now representing the deviation of the rotational angle from the equilibrium position. The terms M, D, K and 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 center and the rotational stiffness K and damping D are both set to zero. Thus the cylinder rotates only due to the action of the fluid forces. The particular setup, along with a variant where the splitter plate is flexible, has been a subject of investigation  by several previous authors (Xu et al., 1990(Xu et al., , 1993Bagheri et al., 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 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 sub-critical and super critical 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 13a and 13b show the baseflow states for the diameter based Reynolds numbers of Re = 45 (sub-critical) and Re = 156 (super-critical). Figure 14 shows the spectra for the two different Reynolds numbers. Both cases have an eigenvalue lying on the positive y-axis (λ = 0.10 + 0i for Re = 45 and λ = 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 monotonic rise in the angle and thus causes the symmetry breaking effect. For Re = 45,  we also simulate the linear and non-linear cases after adding a random small-amplitude perturbations to both flow cases. Figure 15 shows the time evolution of the angle η of the cylinder-splitter body. Both the simulations undergo the same exponential growth of about five orders of magnitude before the non-linear case saturates, while the linear case continues its exponential growth. The saturation of the non-linear simulations occur at η ≈ 15.5 • , 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 non-linear mechanism is required since the flow must equilibrate at the new position. However the onset can be traced to the zero frequency unstable mode.

Structural sensitivity of the eigenvalue
As noted in figure 6 and 8a, 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 2D circular cylinder, which is free to oscillate in the cross-stream direction. The linear FSI problem is defined by equations 26a-26i. 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 (Gianetti and Luchini, 2007;Luchini and Bottaro, 2014). For any linear operator L, the adjoint operator L † is defined such that it satisfies the Lagrange identity for any arbitrary vectors q and p in the domain of L and L † respectively. The symbol ·, · denotes the inner product under which the above identity holds. In the current context the operator L represents the linearized Navier-Stokes equations for FSI, also referred to as the direct problem, L † is the corresponding adjoint operator with the definition of the inner product as the integral over the domain Ω and time horizon T p, q = T Ω p H q dΩ dt.
Defining the vector q = (u, p, η, ϕ), which lies in the domain of L, and an adjoint vector q † = (u † , p † , η † , ϕ † ) which lies in the domain of L † , 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 below.
Here σ † i j is referred to as the adjoint fluid forces acting on the cylinder. Unlike in the direct equations, the adjoint fluid velocities are homogeneous on the FSI interface Γ. The sign of the diffusive term indicates the equations evolve an adjoint field backward in time. A variable substitution τ = −t can be made to return the equations to forward (in τ) propagation. Figure 16: Comparison of the eigenspectra for the direct and adjoint problems for density ratio of 10 and natural frequency ω n = 0.3833 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 16 shows the direct and adjoint spectra for a springmounted 2D cylinder at a diameter based Reynolds number of Re = 50, mass ratio of 10, corresponding to M = 7.854, K = 1.1537 which corresponds to the natural frequency of the spring-mass system of ω n = 0.3833. The damping is set to zero (D = 0). As observed in the figure, the spectra for both the direct and adjoint problems have a very good agreement with each other.
For a general linear eigenvalue problem Lq = λq, a first-order approximation to the perturbed eigenvalue problem (L + ∆L)(q + ∆q) = (λ + ∆λ)(q + ∆q) can be obtained by making use of the corresponding adjoint eigenvector vector field q † (Gianetti and Luchini, 2007), resulting in an expression for the eigenvalue perturbation as where * implies complex conjugation and q † represents the right eigenvector of the adjoint operator L † with the eigenvalue of λ * . To estimate the drift in eigenvalue, knowledge of the specific operator perturbation ∆L is required. However, as shown by Gianetti and 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 a such a case we may estimate the eigenvalue drift along with its upper bound for each spatially localized operator perturbation as 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 Θ(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 spring-mounted cylinder is investigated for varying natural frequencies ω 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 Θ of the least stable eigenvalue is evaluated. Table 5 shows the variation of the natural frequency ω n and the corresponding least stable eigenvalue. The table also shows the unstable eigenvalue for a stationary cylinder at Re = 50. Figure 17 shows the corresponding changes in structural sensitivity of the least stable eigenvalue. When the structural frequency is much lower than the unstable frequency (figure 17b), 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 ω n approaches closer to the unstable frequency (17c). Close to resonance (17d), 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 ω 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 ω 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 ω n = 0.7665 is shown in figure 18 which shows non-zero values of Θ 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 et al., 2008).

Conclusions
An Eulerian formulation for the linear stability analysis of rigid-body-motion fluidstructure-interaction 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); Fernández and Le Tallec (2003a,b) are shown to vanish to a first order-order approximation and play no role in the linearized perturbation dynamics. The formulation is validated by comparing the evolution of the linear equations with the non-linear 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. 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 analyze 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 sub-critical and super-critical Reynolds numbers.
Finally, the eigenvalue sensitivity to structural perturbations is investigated using the adjoint equations for FSI, for a 2D 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.

A Navier-Stokes of Taylor expanded velocity
We evaluate the Navier-Stokes at the perturbed locations using the Taylor expansion of the velocity and pressure fields and then use the geometric linearization to consistently evaluate derivatives at the perturbed configuration. Non-linear terms arising due to interaction of perturbations terms u , p , R and ∆x are dropped.
∂x 0 k ∂x 0 n (I k j I n j − I k j R n j − R k j I n j + R k j R n j ) For clarity the terms are numbered and the expansion of the individual terms after dropping the terms higher than first order is given below.
For linear FSI, we consider the case of a 2D circular cylinder free to oscillate in the vertical direction subject to the action of a spring-mass-damper system and the fluid forces. The structural system is then defined by the position η and velocity ϕ of its center of mass in the vertical direction. In the following i = 2 represents the vertical direction. The vectors are defined as q = (u, p, η, ϕ), which lies in the domain of L, and q = (u † , p † , η † , ϕ † ), which lies in the domain of L † . The Lagrange identity may be wriiten as Integrating by parts and using the divergence theorem leads to The temporal boundary terms for appropriately defined direct and adjoint operators vanish (Bagheri et al., 2009). For the remaining boundary terms, we denote the direct and adjoint stresses in compact notation as For the outer boundaries of the domain the structural variables play no role and and the boundary terms may be written compactly as For homogeneous Dirichlet and stress-free conditions for the direct problem, these lead to On the interface boundary Γ 0 , the structural equations are part of the boundary terms. The baseflow U is identically zero and using the velocity boundary condition at the interface for the direct variables, one obtains The second and the third brackets form the two adjoint structural equations, while the first and the fourth bracket together forms the boundary conditions for the coupling between the fluid and the structure in the adjoint equations. The final set of adjoint equations for FSI of the 2D oscillating cylinder can be written as

C Convergence tests
The convergence tests for an oscillating cylinder are reported in this section. The convergence test is performed for the case of a spring-mounted 2D cylinder at a diameter based Reynolds number of Re = 23.512. The structural parameters are the same as the ones reported in table 1. Table 6 reports the domain size and polynomial order of the various tests performed to ensure converged results. A comparison of the spectra for case 1 and case 2 is shown in figure 19. The panel on the left shows the spectra comparison for two different polynomial orders. The two cases have nearly identical spectra implying the convergence of results with grid resolution. The panel on the right shows the spectra for three different domain sizes. For all three cases the physical eigenvalue of the system is always converged. The highly damped modes represent the box modes of the system which vary slightly depending on the size of the computational box. A similar behavior of the damped eigenmodes can be found in Assemat et al., 2012. Overall, the results show the convergence of the physical eigenvalue both in terms of grid resolution and numerical domain.