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 Reference Ho and Patera1990, Reference Ho and Patera1991) 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 Reference Cossu and Morino2000), $\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. (Reference Fanion, Fernández and Le Tallec2000) and Fernández & Le Tallec (Reference Fernández and Le Tallec2003a,Reference Fernández and Le Tallecb) 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. (Reference Fanion, Fernández and Le Tallec2000) and Fernández & Le Tallec (Reference Fernández and Le Tallec2003a,Reference Fernández and Le Tallecb) 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 (Reference Benjamin1959, Reference Benjamin1960) and Landahl (Reference Landahl1962) 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.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 Reference Fischer, Lottes and Kerkemeier2008). 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 Reference Maday and Patera1989). 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 Reference Negi2017). 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 Reference Schlatter, Stolz and Kleiser2004, Reference Schlatter, Stolz and Kleiser2006) and has been validated with channel flows and flow over wings in Negi et al. (Reference Negi, Vinuesa, Hanifi, Schlatter and Henningson2018). Moving boundaries are treated using the ALE formulation (Ho & Patera Reference Ho and Patera1990, Reference Ho and Patera1991) 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 Reference Fischer, Schmitt and Tomboulides2017). In addition, we have also implemented a fully implicit fixed-point nonlinear iteration method (Küttler & Wall Reference Küttler and Wall2008) 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. (Reference Fischer, Schmitt and Tomboulides2017) 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. Reference Citro, Luchini, Giannetti and Auteri2017). Figure 2 shows the calculated base flow state (streamwise velocity). Sen, Mittal & Biswas (Reference Sen, Mittal and Biswas2009) 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. (Reference Sen, Mittal and Biswas2009).

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 (Reference Eriksson and Rizzi1985) and has been used in several previous works by Barkley, Gomes & Henderson (Reference Barkley, Gomes and Henderson2002), Bagheri et al. (Reference Bagheri, Schlatter, Schmid and Henningson2009b) and Brynjell-Rahkola et al. (Reference Brynjell-Rahkola, Shahriari, Schlatter, Hanifi and Henningson2017) etc. The eigenvalue estimation is done using the implicitly restarted Arnoldi method (Sorensen Reference Sorensen1992) implemented in the open source software package ARPACK (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998). 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 (Reference Cossu and Morino2000) 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 (Reference Cossu and Morino2000) 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 (Reference Cossu and Morino2000). 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 (Reference Cossu and Morino2000), 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 (Reference Cossu and Morino2000) 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 (Reference Navrose and Mittal2016). 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 (Reference Navrose and Mittal2016) 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 (Reference Navrose and Mittal2016) 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 (Reference Navrose and Mittal2016), 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. (Reference Semin, Decoene, Hulin, François and Auradou2012), 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. (Reference Semin, Decoene, Hulin, François and Auradou2012), 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. (Reference Semin, Decoene, Hulin, François and Auradou2012) 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. (Reference Semin, Decoene, Hulin, François and Auradou2012) 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. (Reference Fanion, Fernández and Le Tallec2000) and Fernández & Le Tallec (Reference Fernández and Le Tallec2003a,Reference Fernández and Le Tallecb) 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 Reference Xu, Sen and Gad-el Hak1990, Reference Xu, Sen and Gad-el Hak1993; Bagheri, Mazzino & Bottaro Reference Bagheri, Mazzino and Bottaro2012; Lācis et al. Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014). 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. Reference Bagheri, Mazzino and Bottaro2012). The symmetry breaking phenomenon is known to occur for both subcritical and supercritical Reynolds numbers (Lācis et al. Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014) and in both two- and three-dimensional configurations (Lācis et al. Reference Lācis, Olivieri, Mazzino and Bagheri2017). We investigate the phenomenon through our linear FSI framework. According to the results of Lācis et al. (Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014), 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. (Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014) 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. (Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014) 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. (Reference Pfister, Marquet and Carini2019) 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. (Reference Pfister, Marquet and Carini2019) 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. (Reference Pfister, Marquet and Carini2019) 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 (Reference Pfister2019) 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.