1. Introduction
Fluid-induced vibrations (FIV) are of great interest to many fields of engineering. They are generally classified either as vortex-induced vibrations (VIV), wake-induced vibrations (WIV) or galloping. In the field of engineering, we find two main design philosophies. On the one hand, the design of structures that prevent such vibrations to avoid damage is an obvious example. Griffin & Ramberg (Reference Griffin and Ramberg1982) reviewed studies on VIV of marine risers and listed means to suppress such oscillations. On the other hand, we find the design of oscillating–deformable structures that are conceived to harvest energy. Here, one aims to optimise the motion of the submerged rigid body or deformation of the flexible structure in order to harvest the most energy. Some examples include submerged oscillating–deformable structures that are able to convert energy from marine currents and waves (see review from Bernitsas (Reference Bernitsas2016)). The `wave carpet’ project developed by Alam (Reference Alam2012) aims to extract energy from waves using a deformable carpet. In contrast, the VIVACE (vortex induced vibration aquatic clean energy) concept from Bernitsas et al. (Reference Bernitsas, Raghavan, Ben-Simon and Garcia2008) proposes to convert kinetic energy from marine currents to electricity using VIV of multiple oscillating cylinders. Energy harvesting strategies also find applications in microfluidics: for instance, Lee et al. (Reference Lee, Qi, Zhou and Lua2019) proposed a MEMS (microelectromechanical system) energy harvester based on the oscillation of a cylinder mounted on a piezoelectric chip. In this context, at low Reynolds numbers, they found that the efficiency of the device was greater when placed in a dense field of oscillating cylinders.
The canonical case of a single freely oscillating cylinder has been extensively studied with heavy focus on the lock-in phenomenon (Williamson et al. Reference Williamson and Govardhan2004). It is defined as a synchronisation between the frequency associated with transverse oscillations of the rigid body and that of the vortex shedding in the wake of the cylinder. Outside of the lock-in regime, however, the frequency tends to the vortex shedding frequency of a fixed cylinder. Navrose & Mittal (2016) found that the lock-in phenomenon induces high amplitude vibrations of the cylinder. It has also been shown that a decrease in the reduced mass ratio between the density of the body and the fluid leads to a wider synchronisation regime. Note that most studies with a single cylinder considered a single degree of freedom (1DOF) corresponding to transverse motion. According to Williamson et al. (Reference Williamson and Govardhan2004), in-line motion, if structurally allowed (2DOF), does not change much the dynamics, and mostly turns the transverse oscillation into a lemniscate trajectory where the streamwise motion is induced by a nonlinear effect. On the other hand, pure in-line oscillations seem not to have been observed for a single cylinder. Such motion would be linked to symmetric vortex shedding, which is likely not observed.
A number of configurations involving multiple freely oscillating bodies have been explored. Authors have first focused on WIV. King & Johns (Reference King and Johns1976) first explored WIV of flexible cylinders (2DOF) in tandem traversing a free surface, either rigidly connected or not, for spacings of
$L/D=[1.25 - 7]$
, where
$D$
is the diameter of the cylinders and
$L$
is the distance between the centres, at Reynolds numbers
${\textit{Re}}= {U_\infty D}/{\nu }=[10^3 - 2\times 10^4]$
, where
$U_\infty$
is the inlet velocity and
$\nu$
the fluid’s viscosity. They observed both transverse oscillations and in-line oscillations, but the latter were only reported at much larger values of the Reynolds number, in the range
$[10^3 - 2 \times10^4$
].
Bokaian & Geoola (Reference Bokaian and Geoola1984) focused on the transversal WIV (1DOF) of the rear body by fixing the front one. In the interval
${\textit{Re}} = [2900,5900]$
, they found that the vortex shedding behind the front cylinder is suppressed for spacings of
$L/D \leqslant 2$
. Later studies explored in detail the WIV of a rear oscillating body in the 2DOF (Brika & Laneville Reference Brika and Laneville1999) and 1DOF configurations (Assi et al. Reference Assi, Meneghini, Aranha, Bearman and Casaprima2006, Reference Assi, Bearman and Meneghini2010). Assi, Bearman & Meneghini (Reference Assi, Bearman and Meneghini2010) found that for high separation between bodies, the amplitude of the rear body is decreased and resembles a VIV amplitude. Assi et al. (Reference Assi, Bearman, Carmo, Meneghini, Sherwin and Willden2013) (1DOF) developed the concept of wake stiffness in the galloping of cylinders placed in tandem. The steady lift across the wake is defined as a restoring force towards the centreline, acting as a fluid dynamic spring. The Strouhal number associated with the wake stiffness was found to be constant with the Reynolds number. Mittal & Kumar (Reference Mittal and Kumar2001) numerically studied the tandem and staggered configurations with a 2DOF configuration for low Reynolds number (
${\textit{Re}}=100$
) in the wake interference regime (
$L/D=5.5$
). For this large spacing, the front body behaves like an isolated cylinder with trajectories resembling a lemniscate shape. Soft lock-in was observed and the vortex-shedding frequency of the bodies is detuned from the natural frequency. The rear body displays trajectories in the shape of a lemniscate or a tilted ovoid, whether it is placed in tandem or in staggered configuration. Papaioannou et al. (Reference Papaioannou, Yue, Triantafyllou and Karniadakis2008) used an arbitrary Lagrangian–Eulerian (ALE) method to further explore the effect of spacing on the 2DOF tandem. For a
${\textit{Re}} = 160$
and reduced mass
$m^*=10$
, they explored spacings (
$L/D = [2.5,3.5,5.0]$
) corresponding to different flow regimes in the fixed tandem case (Zdravkovich Reference Zdravkovich1987). Small values of the spacing lead to stronger oscillations of the upstream cylinder over a wider reduced velocity range and shift the response curves to higher reduced velocities. Borazjani & Sotiropoulos (Reference Borazjani and Sotiropoulos2009) directly simulated a tandem of cylinders with 1DOF for a low reduced mass,
$L/D=1.5$
and
${\textit{Re}}=200$
. For low values of the reduced velocity (
$U^*=2 \pi U_{\infty } \sqrt {m_{c} / k}/ D$
,
$k$
and
$m_c$
being the stiffness and cylinder’s mass respectively), they found that the oscillation amplitudes are small and therefore outside of the lock-in region. The front cylinder exhibits larger oscillation amplitudes than the rear one. The effect of an increase of the reduced velocity is to bring the cylinders’ oscillations out of phase, thus increasing their amplitudes of motion. At a critical reduced velocity, the cylinders continue to oscillate out of phase but the rear cylinder’s amplitude becomes greater than the front one. In particular, the authors found a wider lock-in region than for an isolated cylinder. Besides, a structure that would be outside of the lock-in region can be brought into it by placing it in tandem with a similar structure.
Kim et al. (Reference Kim, Alam, Sakamoto and Zhou2009b
) experimentally studied the VIV of the tandem configuration with 1DOF transversal to the fluid flow for several spacings (
$L/D = 1.1 {-} 4.2$
) at
${\textit{Re}} = 4365{-}74\,200$
. Five distinct regimes were identified. Regime I (
$1.1 \lt L/D \lt 1.2$
) features negligible vibrations due to minimal fluctuating lift forces, while Regime II (
$1.2 \lt L/D \lt 1.6$
) exhibits strong vibrations, particularly in the upstream cylinder, for higher reduced velocities. Regime III (
$1.6 \lt L/D \lt 3.0$
) shows significant vibrations of both bodies, with the upstream cylinder’s response being influenced by the downstream cylinder. In Regime IV (
$3.0 \lt L/D \lt 3.7$
), vibrations are again minimal; the downstream cylinder stabilises the wake. Finally, Regime V (
$L/D \gt 3.7$
) displays higher vibrations in the downstream cylinder, attributed to periodic von Kármán vortex shedding. In a subsequent study, Kim et al. (Reference Kim, Alam, Sakamoto and Zhou2009a
) used tripping wires to suppress VIV. They found that placing the wires at an optimal position effectively suppressed vibrations in flow Regimes I–IV by altering the shear layer behaviour and preventing vortex formation.
Prasanth & Mittal (Reference Prasanth and Mittal2009a
,
Reference Prasanth and Mittalb
) numerically examined the free vibrations of two cylinders in the staggered and tandem configurations with 2DOF at
${\textit{Re}}=100$
for
$m^*=10$
,
$L/D =5.5$
and compared the dynamic responses with those of a single cylinder. In the staggered configuration, the upstream cylinder behaves similarly to a single cylinder but with slightly higher oscillation amplitudes, while the downstream cylinder exhibits significantly larger transverse oscillations. Lock-in occurs over a range wider than for a single cylinder, with shared vortex shedding frequencies. The downstream cylinder in the staggered case displays both a lemniscate shape and orbital motions, influenced by complex vortex interactions and asymmetrical flow patterns. For the tandem configuration, the upstream cylinder shows early lock-in and significant influence from the downstream cylinder, despite having a qualitatively similar transverse response to an isolated one. The downstream cylinder experiences much larger oscillations that are twice that of a single cylinder in the laminar regime. Its behaviour mimics high Reynolds number responses, including the presence of an upper branch in vibration response. Both cylinders undergo synchronisation, with frequency and phase shifts tied to vortex shedding and lift forces. Phase differences and hysteresis effects are observed, and the flow regime is divided into different regions based on flow-structure interactions.
Griffith et al. (Reference Griffith, Jacono, Sheridan and Leontini2017) investigated the dynamic response of staggered cylinders at
${\textit{Re}}=200$
with 1DOF, for a fixed streamwise spacing (
$L/D = 1.5$
). They found that gap flow, which reverses direction as the cylinders oscillate, plays a critical role. A regime map was developed, categorising major vortex shedding modes and temporal behaviours. Unlike a single cylinder, matched natural and shedding frequencies do not produce synchronised oscillations; instead, quasiperiodic and chaotic responses emerge. For rigid cylinders, three base modes were observed: no gap flow, gap pair dominated and wake pair dominated – shifting with cross-stream offset. Near the gap–wake pair transition, more complex flow states appear. When cylinders are free to oscillate, low reduced velocities yield minimal motion and rear-cylinder vortex shedding. At intermediate velocities, out-of-phase oscillations enlarge the gap and produce an irregular vortex street. At higher velocities, the rear cylinder chases the front, with joint vortex shedding. As the spacing increases, vortex pairs dominate and the system approaches single-cylinder behaviour. Huera-Huarte & Gharib (Reference Huera-Huarte and Gharib2011) conducted an experimental study of VIV and WIV of a tandem of flexible cylinders (2DOF) in the wake interference regime. They found that both flexible cylinders in a tandem arrangement exhibit classical VIV near lock-in reduced velocities, regardless of the gap distance. At higher reduced velocities, their dynamic responses diverge depending on the spacing between the bodies. The upstream cylinder shows stronger VIV for smaller gaps, while the downstream cylinder may experience WIV at larger gaps, from the presence of vortex shedding in the gap region. Zhang et al. (Reference Zhang, Tang, Lu, Jin, An and Cheng2024a
) investigated FIV (1DOF) of two square cylinders in tandem through simulations and reduced-order modelling. Multiple vibration branches, such as VIV, biased oscillation and galloping, are identified depending on reduced velocity and spacing ratio and their link to wake and structural modes is analysed. Tirri et al. (Reference Tirri, Nitti, Sierra-Ausin, Giannetti and de Tullio2023) and Zhang et al. (Reference Zhang, Lu and Zhang2024b
) conducted a linear stability analysis (LSA) at low Reynolds number of the tandem configuration with 1DOF over a very limited range of structural parameters. They both found two leading unstable eigenmodes, one being associated with the classical vortex shedding behaviour in the wake of a tandem of fixed bodies and the other being of structural nature since they do not have a counterpart in the case of fixed bodies.
Some other authors explored configurations with more bodies and outlined the complexity of such systems. Hosseini, Griffith & Leontini (Reference Hosseini, Griffith and Leontini2021) have investigated the tandem configuration of both two and three fixed cylinders at
${\textit{Re}}=200$
. They found that when the distance between the two bodies is large, the wake exhibits a two-row vortex structure. Adding a third body in that wake has no impact in the majority of cases. However, a region was identified where the placement of an additional body suppresses the vortex shedding in the gap between the two upstream bodies. The tandem of three cylinders oscillating transversely at low Reynolds number was also investigated by Chen et al. (Reference Chen, Ji, Williams, Xu, Yang and Cui2018) and Zhu et al. (Reference Zhu, Zhong, Shao, Zhou and Alam2024).
Note that, in most of the cited bibliography, configurations in which the cylinders are allowed to move in two directions (2DOF) essentially lead to the same dynamics as those in which the cylinder are only allowed to vibrate transversely (1DOF), driven by the transverse force due to antisymmetric vortex shedding. Pure in-line oscillation, which would be linked to symmetric vortex shedding, have not been reported, with the notable exception of King & Johns (Reference King and Johns1976). However, such behaviour was only observed for
${\textit{Re}}\gt 10^3$
, away from the range considered here. Likewise, no known study seems to have reported coupling with the third degree of freedom corresponding to rotational motion of the cylinder. Such motion would be linked to a torque, which is not expected to be significant for a cylinder, as it would be exerted by viscous stress only and not pressure, unlike for more elongated bodies where this kind of motion, known as flutter, is an important subject of research (Chai et al. Reference Chai, Gao, Ankay, Li and Zhang2021).
The consequent number of parameters in a multibody system renders exhaustive parametric studies time-consuming, which justifies the search for new methodologies that enable a systematic scanning of the problem’s parameters. A strategy for analysing fluid–structure interactions has been to project the fluid forces onto the structural degrees of freedom and treat the problem as a harmonically forced structural oscillator, using generalised aerodynamic forces (GAF) to represent the fluid loading. This line of thought was formalised early on by Hassig (Reference Hassig1971), who introduced the
$p{-}k$
method later also referred as Schur complement formulation. In the original formulation, the unsteady aerodynamic forces are available for harmonic oscillations at discrete reduced frequencies. In practical implementations, the use of simplified aerodynamic models (analytical or semiempirical), frequency interpolation for the approximation of the GAFs (Houtman & Timme Reference Houtman and Timme2023), Kriging interpolation (Timme, Marques & Badcock Reference Timme, Marques and Badcock2011) or proper orthogonal decomposition-based (Bekemeyer & Timme Reference Bekemeyer and Timme2019) has historically been a pragmatic way to reduce computational cost.
The idea of exploiting the unsteady forces exerted on a body undergoing prescribed oscillatory motion to predict instability was also pursued by Sabino et al. (Reference Sabino, Fabre, Leontini and Jacono2020). In that work, the unsteady forces were obtained by the exact resolution of the linearised Navier–Stokes equations, rather than from aerodynamic models or frequency interpolation. Subsequently, they defined a mechanical impedance as the ratio between the lift force and the velocity of the cylinders. This allowed them to treat the problem by exploiting an analogy with electric systems (Conciauro & Puglisi Reference Conciauro and Puglisi1981) for which an impedance of the negative real part is a necessary condition for instability. Note that a similar concept of impedance can be applied to acoustic systems, as was done by Fabre et al. (Reference Fabre, Longobardi, Bonnefis and Luchini2019) for an oscillating flow through a thin circular aperture. It was then applied for the stability prediction of the flow through a circular aperture in a thick plate (see Fabre et al. Reference Fabre, Longobardi, Citro and Luchini2020; Sierra-Ausin, Citro & Giannetti Reference Sierra-Ausin, Citro and Giannetti2022).
The objective of this study is to extensively explore the stability of the flow past a tandem of oscillating cylinders, exploiting the impedance method just outlined. As the coupling with the other degrees of freedom is not expected to be significant, we restrict to a 1DOF for each cylinder corresponding to transverse motion. The instability of the tandem is first investigated through the fluid–structure problem in an ALE frame and a LSA is performed for several
$(Re, m^*, U^*)$
. Second, we derive an impedance-based criterion from calculations of the forced case, i.e. where the cylinders are forced to oscillate sinusoidally in the transverse direction of the flow. The results of the coupled case and the impedance-based predictions are compared and we then use the impedance-based method to explore the stability of the tandem for a vast range of parameters.
2. Problem formulation
We consider
$N$
spring-mounted cylinders immersed in a Newtonian, incompressible fluid figure 1. Let
$\tilde {\varOmega }(t)$
denote the (deformable) fluid domain and
$\tilde {\varGamma }_i(t)$
its interface with the cylinders. The flow is physically parametrised by dynamic viscosity
$\nu$
, density
$\rho$
, incoming velocity
$U_{\infty }$
, diameter
$D$
and spacing
$L$
between the cylinders, yielding two dimensionless parameters, the Reynolds number
${\textit{Re}} = {U_\infty D}/{\nu }$
and the spacing ratio
$L/D$
. The fluid flow is described by the velocity and pressure fields
$\tilde {\boldsymbol{u}}, \tilde {p}$
, which are governed, in non-dimensional form, by the following set of equations and boundary conditions at the cylinder walls:
\begin{align} \left .\frac {\partial \tilde {\boldsymbol{u}}}{\partial t}\right |_{\tilde {\boldsymbol{x}}}+(\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}}) \boldsymbol{\cdot }\tilde {\boldsymbol{u}}-\tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }\tilde {\boldsymbol{\sigma }}(\tilde {\boldsymbol{u}}, \tilde {p})&=\boldsymbol{0} \quad \text{in } \quad \tilde {\varOmega }(t), \end{align}
\begin{align} \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }\tilde {\boldsymbol{u}}&=0 \quad \text{ in } \quad \tilde {\varOmega }(t), \end{align}
\begin{align} \tilde {\boldsymbol{u}} & = \dot {Y}_i \boldsymbol{e}_y \quad \text{ on } \quad \tilde {\varGamma }_i(t). \end{align}
Figure 1. Array of
$N$
spring-mounted cylinders with densities, spring stiffness and damping parameters:
$\rho _i$
,
$k_i$
and
$\gamma _i$
. The cylinders are immersed in a fluid domain
$\hat {\varOmega }$
of density
$\rho _f$
. The domain is delimited by inlet and outlet boundaries,
$\varGamma _{in}$
and
$\varGamma _{out}$
as well as lateral boundaries
$\varGamma _l$
.
The symbol
$(\tilde {\hphantom {h}})$
is used for time-dependent quantities as well as time and spatial derivatives evaluated in the time-dependent domain. Here
$(\dot {\hphantom {h}})$
refers to time derivatives. The stress tensor is defined as
$\tilde {\boldsymbol{\sigma }}(\tilde {\boldsymbol{u}}, \tilde {p})=-\tilde {p} \boldsymbol{\text{I}} + 2\mu \tilde {\boldsymbol{D}}(\tilde {\boldsymbol{u}})$
with
$\tilde {\boldsymbol{D}}(\tilde {\boldsymbol{u}}) = (\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}}+\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}}^{ {T}} )$
. The notation
$(|_{\tilde {\boldsymbol{x}}})$
refers to the time derivative in the deformed domain. This term contains the coupling with the motion of the body through the deformation of the domain, handled by the ALE framework as will be explained below. The coupling is also contained in the boundary condition (2.3). In the latter,
$Y_i(t)$
is the instantaneous displacement of the
$i$
th cylinder in the transverse direction.
The cylinders are physically parametrised by their masses
$m_{c_i}$
, spring stiffnesses
$k_i$
and damping parameters
$g_i$
. In a non-dimensional way, this yields three non-dimensional parameters for each cylinder, namely a mass ratio
$m_i^* = ( {4 m_{c_i} }/{\pi D^2 \rho _f})$
, a reduced velocity
$U^*_i = 2 \pi U_{\infty } \sqrt {m_{c_i} / k_i}/ D$
and a damping coefficient
$\gamma _i = g_i/(2\sqrt {m_{c_i}k_i})$
. The equation governing the motion of the
$i$
th cylinder is, in a non-dimensional form,
\begin{align} \ddot {Y}_i+\frac {4 \pi \gamma _i}{U_i^*} \dot {Y}_i+\left (\frac {2 \pi }{U_i^*}\right )^2 Y_i &= \dfrac {4 F_{y_i}(t)}{\pi m_i^*}, \quad \text{ for }i=1,\ldots , N, \end{align}
\begin{align} \mbox{ with } F_{y_i} &= \int _{\tilde {\varGamma }_i} \boldsymbol{e}_{y} \boldsymbol{\cdot }\tilde {\boldsymbol \sigma } \boldsymbol{\cdot }\boldsymbol{n} {\,\text{d}} \tilde {\varGamma }_i. \end{align}
2.1. Arbitrary Lagrangian–Eulerian formulation
2.1.1. General formalism
The ALE method is a conforming method that allows one to treat interfaces in a Lagrangian frame of reference while the fluid is treated in an Eulerian frame of reference. We consider a fixed reference domain
$\hat {\varOmega }$
where unknowns are evaluated in an Eulerian frame of reference. Lagrangian variables on the other hand are evaluated on the actual physical domain
$\tilde {\varOmega }(t)$
, which is time-dependent. Let us define the diffeomorphism
$\mathcal{A}$
that allows us to express the position
$\tilde {\boldsymbol{x}}(t)$
of the actual domain with respect to the position
$\hat {\boldsymbol{x}}$
of the reference domain,
\begin{equation} \begin{aligned} \mathcal{A} : \hat {\varOmega } \times \mathbb{R}^{+} &\longrightarrow \tilde {\varOmega } \times \mathbb{R}^{+}, \\ (\hat {\boldsymbol{x}}, t) &\longmapsto \left (\tilde {\boldsymbol{x}}(t), t\right ). \end{aligned} \end{equation}
This diffeomorphism allows a mapping of the actual domain through the position
\begin{equation} \tilde {\boldsymbol{x}} = \hat {\boldsymbol{x}}+\hat {\boldsymbol{\xi }}_e(\hat {\boldsymbol{x}}, t), \end{equation}
where
$\hat {\boldsymbol{\xi }}_e$
is an extension displacement field that propagates the Lagrangian interface deformation to the fluid domain (as schematised in figure 2). This field is arbitrary and it is determined as a solution of an elliptic equation,
$-\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\mathcal{E}} (\hat {\boldsymbol{\xi }}_e )=\boldsymbol{0}$
, which ensures a smooth distribution over the whole domain. Following the methodology employed by Pfister, Marquet & Carini (Reference Pfister, Marquet and Carini2019), we apply the diffeomorphism to the Lagrangian variables and we substitute them into (2.1) and (2.2) which yields the ALE formulation of the incompressible Navier–Stokes equation in a stress-free configuration,
\begin{align} \hat {J}\left (\hat {\boldsymbol{\xi }}_e\right ) \frac {\partial \hat {\boldsymbol{u}}}{\partial t}+\left ( (\hat {\boldsymbol{\nabla }} \hat {\boldsymbol{u}}) \hat {\boldsymbol{\varPhi }}\left (\hat {\boldsymbol{\xi }}_e\right )\right )\left (\hat {\boldsymbol{u}}-\frac {\partial \hat {\boldsymbol{\xi }}_e}{\partial t}\right )-\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\hat {\boldsymbol{\varSigma }}\left (\hat {\boldsymbol{u}}, \hat {p}, \hat {\boldsymbol{\xi }}_e\right )=\boldsymbol{0} \quad \text{ in } \quad \hat {\varOmega }, \end{align}
\begin{align} -\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\mathcal{E}}\left (\hat {\boldsymbol{\xi }}_e\right )=\boldsymbol{0} \quad \text{ in } \quad \hat {\varOmega }, \end{align}
\begin{align} -\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\left (\hat {\boldsymbol{\varPhi }}\left (\hat {\boldsymbol{\xi }}_e\right ) \hat {\boldsymbol{u}}\right )=0 \quad \text{ in } \quad \hat {\varOmega }. \end{align}
Figure 2. Sketch of the geometrical transformations involved in the ALE approach: (a) actual time-dependent domain and (b) reference domain.
In the previous expressions,
$\hat {\boldsymbol{\nabla }}$
is the gradient in the reference coordinates. Here
$\hat {\boldsymbol{\varPhi }} (\hat {\boldsymbol{\xi }}_e )=\hat {J} (\hat {\boldsymbol{\xi }}_e ) \hat {\boldsymbol{F}} (\hat {\boldsymbol{\xi }}_e )^{-1}$
denotes the deformation operator introduced by the change of variables, with
$\hat {J} (\hat {\boldsymbol{\xi }}_e )=\operatorname {det} (\hat {\boldsymbol{F}} (\hat {\boldsymbol{\xi }}_e ) )$
the Jacobian of the deformation gradient
$\hat {\boldsymbol{F}} (\hat {\boldsymbol{\xi }}_e )=\boldsymbol{\text{I}}+\hat {\boldsymbol{\nabla }} \hat {\boldsymbol{\xi }}_e$
. The ALE fluid stress tensor expressed in the reference configuration writes as
\begin{equation} \hat {\boldsymbol{\varSigma }}\left (\hat {\boldsymbol{u}}, \hat {p}, \hat {\boldsymbol{\xi }}_e\right )=\hat {\boldsymbol{\sigma }}\left (\hat {\boldsymbol{u}}, \hat {p}, \hat {\boldsymbol{\xi }}_e\right ) \hat {\boldsymbol{\varPhi }}\left (\hat {\boldsymbol{\xi }}_e\right )^{ {T}}, \end{equation}
where
$\hat {\boldsymbol{\sigma }} =-\hat {p} \boldsymbol{\text{I}}+\tfrac {1}{Re} \hat {\boldsymbol{D}}$
, with the viscous dissipation tensor defined as
\begin{equation} \hat {\boldsymbol{D}}\left (\hat {\boldsymbol{u}}, \hat {\boldsymbol{\xi }}_e\right )=\frac {1}{2} \frac {1}{\hat {J}\left (\hat {\boldsymbol{\xi }}_e\right )}\left ((\hat {\boldsymbol{\nabla }} \hat {\boldsymbol{u}}) \hat {\boldsymbol{\varPhi }}\left (\hat {\boldsymbol{\xi }}_e\right )+\hat {\boldsymbol{\varPhi }}\left (\hat {\boldsymbol{\xi }}_e\right )^{ {T}}(\hat {\boldsymbol{\nabla }} \hat {\boldsymbol{u}})^{ {T}}\right ). \end{equation}
Hereinafter, we particularise the elliptic operator
$\boldsymbol{\mathcal{E}} = \boldsymbol{\nabla}$
, that is, the extension displacement field is determined by solving a Laplace equation. The complete formulation used to determine the extension field is as follows:
\begin{align} \begin{cases} \Delta \hat {\boldsymbol{\xi }}_e = \boldsymbol{0}, \quad \quad \text{(2.13)}\\ \hat {\boldsymbol{\xi }}_{e} = Y_i \boldsymbol{e}_y \quad \text{on} \quad \hat {\varGamma }_i. \quad \quad \text{(2.14)} \end{cases} \end{align}
2.1.2. The discrete-ALE ansatz
Equation (2.13) is linear, and the number of cylinders is finite; thus, according to the superposition principle, we can look for a solution of the extension field as a function of the cylinders’ vertical displacement,
\begin{equation} \hat {\boldsymbol{\xi }}_e = \sum _{i=1}^{N} Y_i(t) \boldsymbol{\xi }_{e_i}, \end{equation}
where
$\boldsymbol{\xi }_{e_i}$
is an elementary field associated with the displacement of the
$i$
th cylinder, namely the solution of the elementary problem
\begin{align}\smash{\raise-18pt\hbox{$\bigg\{$}} &\Delta \boldsymbol{\xi }_{e_i} = \boldsymbol{0}, \end{align}
\begin{align} &\boldsymbol{\xi }_{e_i} = \boldsymbol{e}_y \quad \text{on} \quad \hat {\varGamma }_i, \end{align}
\begin{align} &\boldsymbol{\xi }_{e_i} = \boldsymbol{0} \quad \text{ on} \quad \hat {\varGamma }_{j, i \neq j}. \end{align}
We call (2.15) the discrete-ALE ansatz. The main advantage of this expression is that it is sufficient to solve
$N$
elementary problems once to reconstruct the deformation field for all possible values of
$Y_i$
, allowing for a reduction of the computational time required to solve the entire system.
2.2. Linearised VIV problem
Following the usual approach, both the fluid and structural variables are decomposed into a steady component and a small-amplitude perturbation which is searched under a modal form. Namely, for the fluid part of the problem, we start with the expansion
\begin{equation} \hat {\boldsymbol{q}}_f(\hat {\boldsymbol{x}}, t) =\boldsymbol{q}_{f,0}(\hat {\boldsymbol{x}}) + \varepsilon {\boldsymbol{q}}_f(\hat {\boldsymbol{x}}) e^{-i \omega t} \end{equation}
were
$\boldsymbol{q}_{f,0}$
is the so-called base flow, corresponding to the steady solution of the Navier–Stokes equations in the reference domain,
$\varepsilon \ll 1$
,
$\boldsymbol{q}_f$
is the fluid part of the eigenmode, and
$\omega$
is an a priori complex eigenvalue. Unlike in other conventions, the perturbation are chosen to be written without a caret symbol. They remain, however, quantities that are evaluated on the reference domain.
Similarly, for the structural part of the problem, we parametrise the displacement of the cylinders by
\begin{equation} Y_i(t) = \varepsilon y_i e^{-i \omega t} ; \qquad \dot {Y}_i(t) = \varepsilon z_i e^{-i \omega t} \quad \mbox{ with } z_i = - i \omega y_i. \end{equation}
The eigenmode of the fluid–structure problem is thus defined as
${\boldsymbol{q}} = ({\boldsymbol{q}}_f, y_1, \ldots ,y_N, z_1, \ldots , z_N )$
, where
${\boldsymbol{q}}_f= ({\boldsymbol{u}},{p} )$
denotes its `fluid’ part and
$[y_1,\ldots ,y_N,z_1,..,z_N]$
its `solid part’.
2.2.1. The ALE fluid–structure coupled formulation
Substituting the ansatz (2.19) into (2.8), (2.10), (2.4) and (2.3) and writing the equations in the steady deformed configuration (see Pfister et al. (Reference Pfister, Marquet and Carini2019), for more details) leads to the following formulations:
\begin{equation} \begin{aligned} &-i\omega \left (\underbrace { {\boldsymbol{u}}}_{\boldsymbol{B}_{\textit{ff}}\, \boldsymbol{q}_f } + \sum _{i=1}^{N} y_i (\underbrace {-\boldsymbol{\xi }_{e_i} \boldsymbol{\cdot }\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0}_{\boldsymbol{B}_{ui} y_i })\right ) = \underbrace {- \left ( {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0 + \boldsymbol{u}_0 \boldsymbol{\cdot }\hat {\boldsymbol{\nabla }} {\boldsymbol{u}} \right ) + 2\mu \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\hat {\boldsymbol{D}}\left ({\boldsymbol{u}}\right )}_{\boldsymbol{A}_{uu} \boldsymbol{q}_f} \underbrace {-\hat {\boldsymbol{\nabla }} p}_{\boldsymbol{A}_{up} \boldsymbol{q}_f} \\ & + \sum _{i=1}^{N} \Bigg [ \underbrace {- \boldsymbol{u}_0 \boldsymbol{\cdot }\Big ( \left (\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0\right ) \big ( (\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\xi }_{e_i}) \boldsymbol{\text{I}} - \hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i}\big ) \Big ) - \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\Big [ p_0 \boldsymbol{\text{I}} \big ( (\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\xi }_{e_i}) \boldsymbol{\text{I}} - \hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i}\big )^T \Big ] \Bigg ] y_i }_{\boldsymbol{A}^{(1)}_{ui} y_i} \\ & + \sum _{i=1}^{N} \Bigg [\! \underbrace {-\mu \hat {\boldsymbol{\nabla }} \!\boldsymbol{\cdot }\!\left ( \!\left (\!\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0 \!\right )\! \left (\!\hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i} \!\right )\! + \!\left (\!\hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i}\!\right )^T \!\left (\!\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0 \!\right )^T\right ) -\! 2 \mu \hat {\boldsymbol{\nabla }} \!\boldsymbol{\cdot }\!\Big ( \!\tilde {\boldsymbol{D}}(\boldsymbol{u}_0) \big ( (\hat {\boldsymbol{\nabla }} \!\boldsymbol{\cdot }\boldsymbol{\xi }_{e_i}) \boldsymbol{\text{I}} - \hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i}\!\big )^T \Big ) \!\Bigg ] y_i }_{\boldsymbol{A}^{(2)}_{ui} y_i}, \end{aligned} \end{equation}
\begin{equation} 0 = \underbrace {\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol{u}}}_{\boldsymbol{A}_{pu} \boldsymbol{q}_f} + \sum _{i=1}^{N} \underbrace { \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\Big ( \big ( \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\xi }_{e_i} \boldsymbol{\text{I}} - \hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_i}\big ) \boldsymbol{u}_0 \Big ) y_i}_{\boldsymbol{A}_{pi} y_i}. \end{equation}
The boundary conditions on the object’s surface are
$$\begin{align} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \begin{cases} {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{e}_y = z_i \quad \text{on} \quad \hat {\varGamma }_i, \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \nonumber\\ {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{e}_y = 0 \quad \text{ on} \quad \hat {\varGamma }_{j, i \neq j} .\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \,\,\,\kern0.5pt \end{cases} \end{align}$$
They are symbolically noted as
$-i\omega \sum \boldsymbol{B}_{\textit{fi}}^* y_i + \boldsymbol{A}_{\textit{ff}}^* {\boldsymbol{q}}_f=0$
, where
$\boldsymbol{B}_{\textit{fi}}^*$
and
$\boldsymbol{A}_{\textit{ff}}^*$
are restriction operators extracting the degrees of freedom localised along the boundaries of the cylinders.
The linearisation of (2.8) and (2.10) thus introduces the operators
$\boldsymbol{A}_{\textit{ff}} = \boldsymbol{A}_{uu} + \boldsymbol{A}_{up} + \boldsymbol{A}_{pu} + \boldsymbol{A}_{\textit{ff}}^*$
and
$\boldsymbol{B}_{\textit{ff}}$
that are purely driven by fluid variables, and the operators
$\boldsymbol{A}_{\textit{fi}} = \boldsymbol{A}^{(1)}_{ui} + \boldsymbol{A}^{(2)}_{ui} + \boldsymbol{A}_{pi}$
and
$\boldsymbol{B}_{\textit{fi}} = \boldsymbol{B}_{ui} + \boldsymbol{B}_{\textit{fi}}^*$
that arise from the interaction of fluid and ALE variables. In this way, (2.8) and (2.10) can be symbolically written with the previously defined operators as
\begin{equation} -i\omega \left ( \boldsymbol{B}_{\textit{ff}}\, \boldsymbol{q}_f + \sum _{i=1}^{N} \boldsymbol{B}_{\textit{fi}} y_i \right ) = \boldsymbol{A}_{\textit{ff}}\, \boldsymbol{q}_f + \sum _{i=1}^{N} \boldsymbol{A}_{\textit{fi}}y_i, \quad \text{ in } \hat {\varOmega }_{f}. \end{equation}
2.2.2. The cylinder’s equations
The lift force
$F_{y_i}$
acting on the
$i$
th cylinder was defined previously in primitive coordinates by (2.5). Using the ALE ansatz and the definition (2.11) of the stress tensor, one is led to an expression of the form
\begin{equation} F_{y_i} = \boldsymbol{F}_{\textit{if}}\, \boldsymbol{q}_f + \sum _{j=1}^{N} {F}_{\textit{ij}}^* y_j .\end{equation}
This expression is composed of two terms. The first is found by integrating on the boundary the stress which is purely linked to the fluid motion,
\begin{equation} \boldsymbol{F}_{\textit{if}}\, \boldsymbol{q}_f = \int _{\hat {\varGamma }_i} \big (-p \boldsymbol{n}+2\mu {\tilde {\boldsymbol{D}}}({\boldsymbol{u}}) \boldsymbol{n}\big ) \boldsymbol{\cdot }\boldsymbol{e}_{y} {\,\text{d}} \hat {\varGamma }_i. \end{equation}
The second component contains the effect of the deformation of the domain associated with the ALE method, and thanks to the discrete-ALE ansatz, it depends only upon the elementary extension fields
$\boldsymbol{\xi }_j$
associated with each of the cylinders,
\begin{eqnarray} {F}_{\textit{ij}}^* &= \displaystyle \int _{\hat {\varGamma }_i} \big (-p_0 \boldsymbol{\text{I}} + 2\mu {\tilde {\boldsymbol{D}}}(\boldsymbol{u}_0)\big ) \big ( \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\xi }_{e_j} \boldsymbol{\text{I}} - \hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_j} \big )^T \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{e}_{y} {\,\text{d}} \hat {\varGamma }_i \nonumber \\&\quad - \displaystyle \int _{\hat {\varGamma }_i} \mu \left ( (\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0) (\hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_j}) + (\hat {\boldsymbol{\nabla }} \boldsymbol{\xi }_{e_j})^T (\hat {\boldsymbol{\nabla }} \boldsymbol{u}_0)^T \right ) \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{e}_{y} {\,\text{d}} \hat {\varGamma }_i. \end{eqnarray}
Introducing (2.19) in (2.4), we obtain the following system for
$i=1,\ldots , N$
:
\begin{align} -i\omega y_i &= z_i, \end{align}
\begin{align} -i\omega z_i &= - \dfrac {4 \pi \gamma _i}{U_i^*} z_i - \Big ( \dfrac {2 \pi }{U^\ast _i} \Big )^2 y_i + \frac {4}{\pi m^*_i} \left ( \boldsymbol{F}_{i,f} \boldsymbol{q}_f + \sum _{j=1}^{N} {F}_{i,j}^* y_j \right ). \end{align}
2.2.3. Eigenvalue formulation for the coupled problem
Considering the coupled problem formulated in terms of the state-vector
$\boldsymbol{q}$
containing both the fluid part
$\boldsymbol{q}_f$
and the solid part
$[y_1,\ldots y_N,z_1,..,z_N]$
, the equations detailed in the two previous subsections can be written in the following matrix system:
\begin{equation} -i\omega \boldsymbol{B} \boldsymbol{q} = \boldsymbol{A} \boldsymbol{q}, \end{equation}
with matrices
\begin{equation} \boldsymbol{B} = \begin{bmatrix} \boldsymbol{B}_{\textit{ff}} & \boldsymbol{B}_{f1} & \ldots & \boldsymbol{B}_{fN} & 0 & \ldots & 0 \\ & 1 & & & & & \\ & &\ddots & & & (0) & \\ & & & 1 & & & \\ & & & & 1 & & \\ & (0) & & & &\ddots & \\ & & & & & & 1 \end{bmatrix}, \end{equation}
and
\begin{equation} \boldsymbol{A} = \begin{bmatrix} \boldsymbol{A}_{\textit{ff}} & \boldsymbol{A}_{f1} & \ldots & \boldsymbol{A}_{fN} & 0 & \ldots & 0 \\ & & & & 1 & & \\ & & (0) & & & \ddots & \\ & & & & & & 1 \\ \tfrac {4} {\pi m^*} \boldsymbol{F}_{1f} & \tfrac {4} {\pi m^*} F_{11}^*-\left( \tfrac {2 \pi }{U^*_1}\right)^2 & & \tfrac {4} {\pi m^*} F_{1N}^* & - \tfrac {4 \pi \gamma _1}{U_1^*}& & \\ \vdots & \vdots & \ddots & \vdots & & \ddots & \\ \tfrac {4} {\pi m^*} \boldsymbol{F}_{1N} & \tfrac {4} {\pi m^*} F_{N1}^* & & \tfrac {4} {\pi m^*} F_{NN}^*-\left( \tfrac {2 \pi }{U^*_N}\right)^2 & & & - \tfrac {4 \pi \gamma _N}{U_N^*} \end{bmatrix}. \end{equation}
2.3. Forced problem and impedance
Besides the resolution of the coupled problem as an eigenvalue problem as just described, we will also make use of an alternative method which consists of first considering the forced problem in which the motion of the cylinders are imposed to behave harmonically, i.e.
$Y_i(t) = y_i e^{-i\omega t}$
with imposed amplitudes
$y_i$
and real frequency
$\omega$
. Thanks to the linearity of the problem (2.25), we can express its solution in compact form as
\begin{equation} \boldsymbol{q}_f = \sum _{j=1}^{N} \boldsymbol{q}_{\textit{fj}} y_j, \quad \mbox{with} \quad \boldsymbol{q}_{\textit{fj}} = - \left [ \boldsymbol{A}_{\textit{ff}} + i\omega \boldsymbol{B}_{\textit{ff}} \right ]^{-1} \left [ \boldsymbol{A}_{\textit{fj}} + i\omega \boldsymbol{B}_{\textit{fj}} \right ]. \end{equation}
In practice, each
$\boldsymbol{q}_{\textit{fj}}$
is the solution of the forced problem (2.21) considering a unitary displacement of the
$j$
th cylinder, namely
\begin{equation} y_j = 1 ; \quad \quad z_j = -i\omega ; \quad y_{i\ne j} =z_{i\ne j} =0. \end{equation}
One can now introduce the decomposition (2.34) into the definition (2.26), of the lift forces acting on the
$i$
th body. Using the operators defined in (2.27) and (2.28), this leads to
\begin{equation} F_{y_i} = \sum F_{\textit{ij}} y_j \quad \text{where} \quad F_{\textit{ij}} = \boldsymbol{F}_{\textit{if}}\, \boldsymbol{q}_{\textit{fj}} + F_{\textit{ij}}^* .\end{equation}
Each of the terms
$F_{\textit{ij}}$
can be considered as a transfer function, corresponding to the ratio of the force exerted on the body
$i$
to the displacement of the body
$j$
. Note that each term
$F_{\textit{ij}}$
depends only upon the elementary displacement field
$\boldsymbol{\xi }_{e_j}$
and the elementary solution
$\boldsymbol{q}_{\textit{fj}}$
of the forced problem, both calculated by imposing
$y_j = 1$
,
$z_j = -i\omega$
,
$y_{i \neq j} = 0$
and
$z_{i \neq j} = 0$
. In other words, we impose the movement of the
$j$
th cylinder and fix all others in order to calculate
$F_{\textit{ij}}$
.
Rather than this definition as a transfer function, it turns out to be more physically significant to define an impedance
$Z_{\textit{ij}}$
relating the force on body
$i$
to the opposite of the velocity of body
$j$
, hence the definition
\begin{equation} Z_{\textit{ij}} = - \frac {F_{\textit{ij}}}{z_j} \equiv \frac {F_{\textit{ij}}}{i \omega }. \end{equation}
Note that this definition is equivalent to the one used in Sabino et al. (Reference Sabino, Fabre, Leontini and Jacono2020), except for a factor of 2 due to the fact that they founded their definition upon the lift coefficient
$C_y$
instead of the dimensionless force
$F_y = C_y/2$
.
In a compact form, the impedance
$Z_{\textit{ij}}$
can also be expressed in terms of the previously introduced operators as
\begin{equation} Z_{\textit{ij}} = (i \omega )^{-1} \left \{ \boldsymbol{F}_{\textit{if}} \left [ \boldsymbol{A}_{\textit{ff}} + i\omega \boldsymbol{B}_{\textit{ff}} \right ]^{-1} \left [ \boldsymbol{A}_{\textit{fj}} + i\omega \boldsymbol{B}_{\textit{fj}} \right ] + F_{\textit{ij}}^* \right \} . \end{equation}
2.4. Generalised impedance criterion for a tandem of cylinders
We will now focus on the case of a tandem of cylinders (
$N=2$
). Solving the forced problem for the front and rear cylinder will, respectively, give the impedances
$Z_{11}$
,
$Z_{21}$
and
$Z_{12}$
,
$Z_{22}$
. We can plug (2.28) along with the definition of the impedances into the harmonic oscillator (2.29) of the fluid–structure problem, which yields
\begin{equation} \begin{cases} \left (-\omega ^2 - \dfrac {4\pi \gamma _1}{U_1^*} i \omega + \left ( \dfrac {2\pi }{U_1^*} \right ) ^2 \right )y_1 = \dfrac {i \omega }{\pi m_1^*} 4\left ( Z_{11}y_1 + Z_{21}y_2 \right )\!, \\[3ex] \left (-\omega ^2 - \dfrac {4\pi \gamma _2}{U_2^*} i \omega + \left ( \dfrac {2\pi }{U_2^*} \right ) ^2 \right )y_2 = \dfrac {i \omega }{\pi m_2^*} 4\left ( Z_{12}y_1 + Z_{22}y_2 \right )\!. \end{cases} \end{equation}
Building the matrix
\begin{equation} Z_T = \begin{bmatrix} - \omega ^2 - \dfrac {4\pi \gamma _1}{U_1^*} i \omega + \left ( \dfrac {2\pi }{U_1^*} \right ) ^2 - \dfrac {i \omega 4 Z_{11}}{\pi m_1^*} & - \dfrac {i \omega 4 Z_{21}}{\pi m_1^*} \\ - \dfrac {i \omega 4 Z_{12}}{\pi m_2^*} & -\omega ^2 - \dfrac {4\pi \gamma _2}{U_2^*} i \omega + \left ( \dfrac {2\pi }{U_2^*} \right ) ^2 - \dfrac {i \omega 4 Z_{22}}{\pi m_2^*} \end{bmatrix}, \end{equation}
the (2.39) can be condensed as the following system:
\begin{equation} Z_T \boldsymbol{\cdot }\begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = 0, \end{equation}
Finally, we define a generalised impedance function as
\begin{equation} H(\omega ) = \text{det}(Z_T). \end{equation}
It is an analytical function of the complex frequency
$\omega =\omega _r+i\omega _i$
.
At this point, we can remark that complex roots of (2.42), corresponding to non-trivial solutions of the two-dimensional system (2.41), also correspond to solutions of the fluid–structure eigenvalue problem (2.31). We thus have replaced the resolution of a matrix eigenvalue problem of large dimension by the sole inspection of a
$2\times 2$
matrix. This is, however, a nonlinear eigenvalue problem since
$\omega$
appears quadratically in
$Z_T$
and most importantly because
$Z_{\textit{ij}}$
depends on
$\omega$
. In practice, the computation of all physically relevant eigenvalues requires the knowledge of the functions
$Z_{\textit{ij}}$
in the whole complex
$\omega$
-plane. However, if one is only interested in localising the marginally stable states, it is only required to know the values of these functions along the real
$\omega$
axis. This property is at the origin of a very efficient method which will be explained and validated in § 3, and subsequently used to perform parametric studies in § 4.
The impedance-based method developed in the present work is an exact realisation of the classical
$p{-}k$
method (Hassig Reference Hassig1971), as it provides a fully coupled description of the linear fluid–structure interaction. Instead of relying on approximate aerodynamic transfer functions, the fluid response is obtained directly by solving a series of forced problems, in which the bodies are imposed to oscillate harmonically at prescribed real frequencies. This procedure yields frequency-dependent impedance functions that rigorously encapsulate the hydrodynamic feedback of the flow on the moving structures. The resulting impedance matrix, combined with the structural parameters, defines a compact stability criterion from which the onset of instability can be predicted without explicitly solving the coupled eigenvalue problem. The approach thus retains the interpretability of classical GAF-based methods while extending the concept to the exact hydrodynamic coupling derived from the full linearised Navier–Stokes equations, enabling accurate and efficient parametric stability analyses across a wide range of configurations.
The advantage of the criterion described here is that a limited number of calculations is required in order to get the stability prediction of a vast number of different cases. Once a set of forced problems are calculated for a fixed set of {
${\textit{Re}}$
and
$L$
}, the stability of systems for any set of {
$U^*_1$
,
$U^*_2$
,
$m^*_1$
,
$m^*_2$
,
$\gamma ^*_1$
,
$\gamma ^*_2$
} is acquired by the simple inspection of the determinant of
$2\times 2$
matrices.
2.5. Numerical implementation
The equations are rewritten in a variational formulation, spatially discretised and solved with the FreeFem++ open-source software (Hecht Reference Hecht2012). The problem is necessarily formulated along a truncated domain of sufficient dimension (see figure 1), and thus the equations are complemented with suitable boundary conditions for the external boundaries (
$\varGamma _{in}$
,
$\varGamma _{l}$
,
$\varGamma _{out}$
). For the fluid variables, a Dirichlet boundary condition is imposed at the inflow boundary
$\varGamma _{in}$
:
$\boldsymbol{u}_0 |_{\varGamma _{in}} = U_\infty \boldsymbol{e}_x$
with
$ U_\infty \equiv 1$
, and a stress-free condition is imposed on the lateral and outflow boundaries. For the ALE variables, homogeneous Dirichlet boundary conditions are imposed on all outer boundaries. Following the classical procedure, the base-flow is computed using a Newton method and the eigenproblems are solved using a shift-invert method as implemented in the SLEPc library.
Two additional tricks are employed to lighten the resolution and/or improve the precision. First, mesh adaptation is intensively used to increase the mesh density in regions of strong gradients while decreasing it in other regions. Second, for the linearised problem, we employ the complex mapping method (Sierra, Fabre & Citro Reference Sierra, Fabre and Citro2020) to characterise the stability properties of the problem and to suppress artificial unstable modes arising due to the strongly convective nature of the wake. Such a method has been successfully employed in the past in other fluid configurations, for instance, the jet flow past a circular aperture (Sierra-Ausin et al. Reference Sierra-Ausin, Citro and Giannetti2022), a flow of two coaxial jets (Corrochano et al. Reference Corrochano, Sierra-Ausín, Martin, Fabre and Le Clainche2023) or the wake flow past a rotating particle (Sierra-Ausín et al. Reference Sierra-Ausin, Citro and Giannetti2022). When using the complex mapping method, the spatial structure of the global mode near the boundary becomes evanescent and does not have an influence on the stability properties of the problem.
Monitoring of all computations and postprocessing is done thanks to the StabFem interface (Fabre et al. Reference Fabre, Citro, Ferreira Sabino, Bonnefis, Sierra, Giannetti and Pigou2018). Following to the philosophy of this project, sample codes reproducing key results of the present paper are available on the website of the project (https://stabfem.gitlab.io/StabFem/).
Regarding the method for threshold detection in terms of the impedance concept explained in § 2.4, the numerical resolution procedure consists, in a first step, of generating a tabulation of the impedance functions
$Z_{\textit{ij}}$
as function of
$\omega$
and
${\textit{Re}}$
. Then, zeros of the generalised impedance (2.42) are computed by considering it as two functions (real and imaginary parts of
$H(\omega )$
as function of the two variables
$\omega$
,
$U^*$
, and a Levenberg–Marquardt method is used to solve it, using interpolation along the range of tabulated
$\omega$
to evaluate the impedances
$Z_{\textit{ij}}$
and their
$\omega$
-derivatives.
If one looks at the computational cost of traditional LSA (Arnoldi method) on a standard desktop computer on a single core, solving for
$20$
eigenvalues in our configuration takes approximatively
$20$
s. Let us assume that one needs a total of
$10$
computations in order to find the thresholds, which brings a total of
$200$
s to find the thresholds for a selected set of parameters {
${\textit{Re}}$
,
$L$
,
$U^*_1$
,
$U^*_2$
,
$m^*_1$
,
$m^*_2$
,
$\gamma ^*_1$
and
$\gamma ^*_2$
}. The impedance-based method on the other hand, requires the computation of a set of forced problems, one of them taking approximately
$5$
s. Considering that one needs tabulated values of
$\omega =[0.4:0.01:1.2]$
, that brings the computational cost to
$405$
s. These calculations, however, are only needed for a set of {
${\textit{Re}}$
and
$L$
}. From there, finding the zeros of the impedance function for any set of {
$U^*_1$
,
$U^*_2$
,
$m^*_1$
,
$m^*_2$
,
$\gamma ^*_1$
and
$\gamma ^*_2$
} is instantaneous.
3. Validation
Throughout this article, we will consider a tandem of spring-mounted cylinders (with the exception of Appendix B). The Reynolds numbers investigated range up to
${\textit{Re}}=100$
and the damping parameters of the cylinders are considered to be the same and noted
$\gamma =\gamma _1=\gamma _2$
. The damping ratios are set to
$0$
, except for § 4.2.3. The reduced velocity and reduced mass of both cylinders will be considered to be the same and will be noted as
$U^*=U_1^*=U_2^*$
and
$m^*=m_1^*=m_2^*$
, except for part of § 4.2.5.
3.1. Linear coupled problem
Figure 3 shows the real and imaginary parts of the leading eigenvalues against
$U^*$
for
$L=1.5$
at a Reynolds number of
${\textit{Re}}=100$
and for
$m^*=2.546$
and
$m^*=20$
. Two leading unstable modes are found and the evolution of their growth rate and frequency is in very good agreement with Zhang et al. (Reference Zhang, Lu and Zhang2024b
). For
$m^*=2.546$
, we observe a slight discrepancy with results from Tirri et al. (Reference Tirri, Nitti, Sierra-Ausin, Giannetti and de Tullio2023). These authors used an immersed boundary method which was shown to yield incorrect results when the added-mass effect is not properly taken into account (Suzuki & Inamuro Reference Suzuki and Inamuro2011).
Figure 3. Real and imaginary parts of the leading eigenvalues (LSA) with respect to
$U^*$
at
${\textit{Re}} = 100$
and
$L=1.5$
for
$m^* = 2.546$
(a,c) and
$m^* = 20$
(b,d). Plain lines are the results of the current study. Results from Tirri et al. (Reference Tirri, Nitti, Sierra-Ausin, Giannetti and de Tullio2023) are shown by 
and 
. Results from Zhang et al. (Reference Zhang, Lu and Zhang2024b
) are shown by 
and 
. The unstable region is depicted as the grey zone. The natural frequency of a spring-mounted cylinder in vacuum
$\omega _n=({2\pi }/{U_n^*})$
is shown as 
. The growth rate and frequency of the fluid mode behind two fixed cylinders are displayed as 
. The predictions from the impedance criterion are shown as 
.
3.2. Impedance-based stability predictions
The generalised impedance provides a criterion of instability: the system is unstable if there exists a non-trivial solution of (2.41), for which
$\omega _i\gt 0$
, as described in Fabre et al. (Reference Fabre, Longobardi, Citro and Luchini2020). The link between impedance and instability can be formulated using Nyquist diagrams. Figure 4 shows an example of the threshold detection based on impedance criterion for
$L=1.5$
,
${\textit{Re}}=100$
and
$m^*=2.546$
. Figure 4(a) shows the zero isolines of both the real and imaginary parts of the function
$H$
in the
$\omega - U^*$
plane. The intersection of these lines gives the predictions of
$U^\ast$
and
$\omega$
at which the system is neutrally stable (noted as 
). These predictions are reported in figures 3(a) and 3(c) as points P1 and P2, and are in perfect agreement with the results from linear stability. Figure 4(b) shows the Nyquist diagrams, i.e. the imaginary part of the function
$H$
against its real part for different values of
$U^*$
. For the first point of neutral stability (at
$U^*=4.38$
), one eigenvalue switches from stable to unstable when increasing the reduced velocity
$U^*$
, as the Nyquist curve transitions from encircling the origin to having the origin lying at the right-hand side of the curve’s trajectory, relative to its parametrisation from
$\omega = 0$
to
$\omega \to \infty$
(see figure 4b i). However, for the second point of neutral stability (at
$U^*=5.94$
), one eigenvalue goes from unstable to stable when increasing the reduced velocity
$U^*$
, and the Nyquist curve transitions from having the origin lying to the right-hand side of the curve’s trajectory to the curve encircling the origin (see figure 4b ii). The impedance-based detection for
$m^*=20$
is also in perfect agreement with LSA results as shown in figures 3(b) and 3(d) by the points P3 and P4.
Figure 4. Illustration of the threshold detection based on the impedance criterion. (a) Zero isolines of the real (
) and imaginary (
) parts of the impedance function
$H$
in the
$\omega - U^*$
plane. The symbols 
show the zeros of
$H$
. (b) Plot of the imaginary part of the impedance function
$H$
with respect to its real part (Nyquist curve) for several values of
$U^*$
. Nyquist curves around the points of neutral stability at
$U^*=4.28$
(bi) and
$U^*=5.94$
(b ii).
4. Results and discussion
4.1. Description of the modes
Let us first come back to the cases
$\{L=1.5; Re=100; m^*=2.546\}$
and
$\{L=1.5; Re=100; m^*=20\}$
which were previously used for validation and plotted in figure 3, and comment on them from a physical point of view. To explain the physical significance of these results, we will also plot in figure 5 the structure of the modes for selected sets of parameters. Each eigenmode has been normalised by setting the highest velocity
$z_i$
to
$1$
. The transversal fluid velocity
$u_y$
and the vertical velocities of each cylinder
$z_1$
and
$z_2$
are shown in figure 5 for
${\textit{Re}}=100$
and for
$m^*=2.546$
. The real and imaginary part of the velocities
$z_i$
are, respectively, plotted as an arrow with a triangle (
) and a circle head (
). For the sake of readability, these amplitudes have been doubled when plotted on the following figures. With this choice of normalisation, large values of the transversal velocity
$u_y$
in the bulk (indicated by the colour ranges of the panels) imply that the transversal motion of the cylinder associated with the mode is much weaker than the transversal fluid motion. On the other hand, a low maximum in the transversal fluid velocity implies a strong transversal motion of the cylinder.
Figure 5. Transversal velocity
$u_y$
fields of the flow past a tandem of cylinder spaced of
$L=1.5$
for
${\textit{Re}}=100$
and
$m^*=2.546$
: mode B (a–c) and mode A (d–f) at reduced velocities
$U^*=3$
(a,d),
$U^*=5$
(b,e) and
$U^*=8$
(d,f). See explanations in text for significance of the colour maps and representations of the cylinder’s displacement.
As reported in Tirri et al. (Reference Tirri, Nitti, Sierra-Ausin, Giannetti and de Tullio2023) and Zhang et al. (Reference Zhang, Lu and Zhang2024b
), the results of LSA for
$m^*=2.546$
introduces two leading eigenmodes classified as A (
) and B (
).
Mode A is unstable at low reduced velocities with its growth rate matching that of the leading mode behind the fixed tandem of cylinders
$\omega _{f_i}=0.039$
(
in figure 3). It becomes stable at a reduced velocity of
$U^*=5.94$
and its frequency varies little and matches the frequency of the fixed tandem mode
$\omega _{f_r}=0.76$
(
), with a small decrease towards the highest values of the reduced velocity. This mode has been classified as a `fluidic’ mode in the literature. Mode B has a more complex structure arising from the fluid–structure interaction. The mode becomes unstable at
$U^*=4.28$
with a maximum growth rate at
$U^*=7$
. For low reduced velocities, the frequency of the mode matches the natural frequency of a structure-only system (
$\omega _n={2\pi }/{U_n^*}$
, shown as 
in figure 3). For higher reduced velocities, its frequency tends to that of the fixed tandem mode. Both cylinders display transverse motion throughout the whole range of reduced velocities, confirming the structural nature of the mode. In accordance with the results of Tirri et al. (Reference Tirri, Nitti, Sierra-Ausin, Giannetti and de Tullio2023), the flow field of mode A resembles that of the wake of two fixed cylinders, as can be seen in figure 5(d–f). The large values of the transversal velocity (as seen for the colour ranges used in the panels) confirm the fluid-dominated nature of this mode. The region of high transversal velocity is shifted downstream with increasing reduced velocities.
On the other hand, the wake of mode B includes a high transversal velocity region localised around the bodies, especially at low reduced velocities (as seen in figure 5
a). At higher reduced velocity, however, the high transversal velocity region is shifted downstream, resembling the structure behind a fixed tandem of cylinders. Correspondingly, for low reduced velocity, the low values of the transversal velocity show that the transversal motion of the cylinder is dominant compared with the fluid’s motion (figure 5
a). Increasing the reduced velocity, the value of the transversal velocity increases slightly as the mode transitions to a more fluid-dominated nature (figures 5
b and 5
c). We can also note that generally, the motion of the rear cylinder is greater than of the front, except at low reduced velocity (figure 5
a). Navrose & Mittal (2016) found, for a single oscillating cylinder, that the high frequency of the fluid-elastic modes at low
$U^*$
was linked to a high vorticity region close to the body. On the other hand, a low frequency induced a shift of the high vorticity region downstream of the body. We observe the same link between frequency and structure of the wake for the tandem system.
For
$m^*=20$
at
${\textit{Re}}=100$
(see figures 3
b and 3
d), mode A remains unstable for all reduced velocities investigated. Its growth rate matches that of the fixed tandem mode only to depart from it around
$7\lt U^*\lt 10$
with a maximum at
$U^*=8.2$
. Its frequency is the same as that of the fixed tandem mode
$\omega _{f_r}$
and the displacement of both cylinders is null. Concerning mode B, the mode is unstable for
$5.425\lt U^*\lt 7.6$
with a maximum in growth rate at
$7.15$
. At
${\textit{Re}}=100$
, its frequency follows very closely that of the natural structural-only system
$\omega _n$
, as shown in figure 3. The structure of the transversal velocity of the different modes is very similar to that at a lower mass ratio. As the Reynolds number decreases, some exchanges of stability between the modes are observed, as reported in § 4.2.
4.2. Parametric study
4.2.1. Neutral curves for
$m*=2.5$
and
$20$
;
$L = 1.5$
and
$3$
Figure 6 details the neutral stability curves in the
${\textit{Re}}-U^*$
plane, obtained for two values of
$m^*$
and two values of
$L$
. Note that the figure displays both results obtained through the resolution of the eigenvalue problem, i.e. LSA, (
, 
and 
), and results obtained with the impedance-based method (
, 
and 
). The excellent agreement of the two methods gives a further validation of the impedance-based method, which will be mostly used in the subsequent section for further parametric studies in a larger range of parameters. Indeed, as already explained, once the impedance functions have been previously calculated and tabulated, one is able to generate results for all values of the structural parameters with no additional cost than simply solving a
$2\times 2$
linear system.
Figure 6. Regions of instability (shaded colours) in the
${\textit{Re}}-U^*$
plane for (a,b)
$L=1.5$
and (c,d)
$L=3$
; and for (a,c)
$m^*=2.5$
and (b,d)
$m^*=20$
. Full lines are the results from the impedance-based predictions and markers are results from LSA. Light grey indicates regions where one unstable mode exists and dark grey regions where two unstable modes exist.
Consider, first, the situation for
$\{L=1.5; m^*=2.5\}$
(figure 6
a). At this spacing and for a low reduced mass of
$m^*=2.5$
, mode A (
) is stable for all
$U^*$
below
${\textit{Re}}=75$
. Above that Reynolds number, the mode is unstable for values of the reduced velocity below
$U^*=6$
. Concerning mode B (
), it is unstable for Reynolds numbers above
${\textit{Re}}=18$
, spanning a broad range of reduced velocities from
$U^*\approx 4$
to
$U^*\approx 18$
. For Reynolds numbers above
${\textit{Re}}=80$
, mode B is unstable for all reduced velocities above
$U^*\approx 4$
. Above
${\textit{Re}}\approx 78$
, the two unstable modes coexist in a region which is indicated by darker grey shading in figure 6(a).
Consider, now, the case
$\{L=1.5,m^*=20\}$
(figure 6
b). For this set of parameters, the mapping of the instability regions is more complex to describe, since an exchange occurs between the two leading branches. Namely, for the largest values of
${\textit{Re}}$
considered in the figure, an A mode exists for the whole range of
$U^*$
and a B mode exists in a range of
$U^*$
centred around the value
$U^*=7$
. This is consistent with the results previously shown in figure 3(b) for
${\textit{Re}}= 100$
where the threshold of the mode A branch is displayed in blue while the one of the mode B branch is displayed in red. When decreasing the Reynolds number, a topological transition occurs, leading to a situation where the mode A branch for small
$U^*$
becomes connected with what was previously the mode B branch and vice versa, leading to a situation similar to what was displayed in figure 3(a). This transition occurs for
${\textit{Re}} \approx 84$
, at which value one finds a double root of the eigenvalue problem. This exceptional point occurs at
$U^* \approx 6.8$
and is identified by a marker in the figure 6(b). It is codimension-two but does not correspond to a bifurcation point, since the coincidence occurs in the unstable region. Due to this exchange of branches, it becomes difficult to distinguish the two modes in the whole range of parameters. When decreasing further the Reynolds number, another peculiar feature appears in the stability map. Namely, the right-hand part of the neutral curve displays a small loop, between the values
${\textit{Re}} = 80$
and
${\textit{Re}}= 75$
. Within this loop, two unstable modes exist. A detailed inspection shows that in this region there exists a codimension-two point where a topological transition occurs. Both exceptional points are represented as 
in figure 6(b). This complex behaviour is detailed in the Supplementary material (see Appendix A).
Similarly, neutral curves for
$\{L=3;m^*=2.5\}$
and
$\{L=3;m^*=20\}$
are, respectively, plotted in figures 6(c) and 6(d). Three leading eigenmodes are found: mode A (
), mode B (
) and mode C (
). Figure 7 shows the transversal velocity field of the leading eigenmodes at
$m^*=2.5$
and
${\textit{Re}}=100$
at different reduced velocities, following the normalisation explained previously. For
$m^*=2.5$
, mode A is unstable for
${\textit{Re}}\gt 80$
in the low reduced velocity range. At
${\textit{Re}}=100$
, the frequency of the mode is that of the fixed tandem configuration and the shape of the mode also resembles the wake mode behind the fixed tandem configuration (see figure 7
a–c). Moreover, the displacement of both of the cylinders is negligible, as can be seen from the values of the transversal velocity. Mode B, on the other hand, is unstable over a wide range of reduced velocities (above
$U^*\approx 4$
), for
${\textit{Re}}\gt 18$
. At
${\textit{Re}}=100$
and for low reduced velocities, the high transversal velocity region is localised around the front body. Both cylinders exhibit high displacement as seen in figure 7(d). At the same time, the frequency of the mode follows that of the natural frequency of the structure-only spring-mounted system. Increasing the reduced velocities, the displacement of the rear cylinder becomes greater than the front one and the high transversal velocity region is shifted downstream (see figures 7
b and 7
c). The frequency of the mode then tends to that of the fixed tandem configuration. Mode C is unstable for
${\textit{Re}}\gt 48$
, over a limited range of reduced velocities, between
$6\lt U^*\lt 11$
. At
${\textit{Re}}=100$
, its frequency is that of the fixed tandem configuration. For low reduced velocities, the shape of the mode also resembles the wake mode behind the fixed tandem configuration (see figure 7
g). However, towards higher reduced velocities, the rear cylinder exhibits large vertical displacement as can be seen from the low values of the transversal velocity in figures 7(h) and 7(i). For
$m^*=20$
, the branch of neutral stability of mode A is shifted towards higher reduced velocities and the ranges of reduced velocities over which modes B and C are unstable are reduced.
Figure 7. Transversal velocity
$u_y$
fields of the flow past a tandem of cylinder spaced of
$L=3$
for
${\textit{Re}}=100$
and
$m^*=2.5$
: mode A (a–c), mode B (d–f) and mode C (g–i) at reduced velocities
$U^*=4$
(a,d,g),
$U^*=8$
(b,e,h) and
$U^*=12$
(d,f,i). See explanations in text for the significance of the colour maps and representations of the cylinder’s displacement.
4.2.2. Effect of the mass
After having computed the forced problem for
$L=1.5$
and
$L=3$
over a range of Reynolds number
${\textit{Re}}=[5-100]$
(by steps of
${\textit{Re}}=1$
), we apply the impedance-based criterion for a range of reduced masses going from
$m^*=0.01$
to
$m^*=100$
. The damping parameter of both cylinders is set to zero. Figure 8 shows the neutral curves in the
${\textit{Re}}-U^*$
plane predicted by the impedance-based criterion. For
$L=1.5$
(see figure 8
a), increasing the reduced mass of the bodies decreases the range of reduced velocities over which mode B is unstable. The stability threshold of mode A is shifted towards higher reduced velocities when increasing the reduced mass, extending the range over which the mode is unstable. The exchange of stability observed at
${\textit{Re}} \approx 84$
for
$m^*=20$
described in § 4.2.1 is also observed for higher masses. For
$L=3$
(see figure 8
b), the stability threshold of mode A is also shifted towards higher reduced velocities and the range over which mode B is unstable is decreased. Concerning mode C, increasing the reduced mass also restricts the range of reduced velocities over which the mode is unstable. Generally speaking, increasing the mass has a stabilising effect on modes B and C and a destabilising effect on mode A.
Figure 8. Curves of neutral stability in the
${\textit{Re}}-U^*$
plane for (a)
$L=1.5$
and (b)
$L=3$
, computed from the impedance-based method. The colour codes for the different modes are similar to figure 6. The colour intensity is going from light to dark for increasing masses whose values are reported in corresponding colours.
4.2.3. Effect of the damping ratio
From the same calculations of the forced problem, the impedance-based criterion is used for a range of damping ratios going from
$\gamma =0$
to
$\gamma =1$
, keeping the reduced mass fixed to
$m^*=2.5$
. Figure 9 shows the neutral curves in the
${\textit{Re}}-U^*$
plane predicted by the impedance-based criterion. For
$L=1.5$
(see figure 9
a), mode A remains largely unaffected by light damping (
$\gamma =0.1$
). On the other hand, the range of reduced velocities over which mode B is unstable gets narrower with increasing damping. The onset of instability of mode B is also shifted towards higher Reynolds numbers, going from
${\textit{Re}}=18$
for
$\gamma =0$
to
${\textit{Re}}=60$
for
$\gamma =1$
. The system exhibits a critical transition at between
$\gamma =0.1$
and
$\gamma =0.13$
at around
$U^*\approx 6$
and
${\textit{Re}}\approx 81$
, where the branches coalesce. Figures 10(a) and 10(b), respectively, show the growth rates and frequencies of the modes before and after the higher codimension point, for a fixed
${\textit{Re}}$
number. The frequency of the coalescing modes being the same at the exceptional point, the latter corresponds to a codimension-three point (double Hopf with strong resonance) and is plotted as 
in figure 9(a).
Figure 9. Curves of neutral stability in the
${\textit{Re}}-U^*$
plane for a low reduced mass of
$m^*=2.5$
and for (a)
$L=1.5$
and (b)
$L=3$
, computed from the impedance-based method. The colour codes for the different modes are similar to figure 6. The colour intensity is going from light to dark for increasing damping ratios, whose values are reported in corresponding colours.
Figure 10. Growth rates (c) and frequencies (d) of the leading eigenmodes (LSA) for
$L=1.5$
and
$m^*=2.5$
at
${\textit{Re}}=81$
. The modes for
$\gamma =0.1$
are plotted with (
) and (
) and the coalesced mode for
$\gamma =0.13$
is plotted with (
). The unstable region is depicted as the grey zone. The natural frequency of a spring-mounted cylinder in vacuum
$\omega _n={2\pi }/{U_n^*}$
is shown as 
. The growth rate and frequency of the fluid mode behind two fixed cylinders are displayed as 
.
For
$L=3$
(figure 9
b), the threshold of mode A is slightly shifted towards higher reduced velocities by light damping and the ranges over which modes B and C are unstable are reduced by increasing the damping. Two transitions are observed, which are similar to the one observed for
$L=1.5$
. Between
$\gamma =0.1$
and
$\gamma =0.25$
, around
${\textit{Re}}\approx 95$
and
$U^*\approx 5$
, mode A and the low-
$U^*$
branch of mode C coalesce. Similarly, between
$\gamma =0.45$
and
$\gamma =0.5$
, around
${\textit{Re}}\approx 75$
and
$U^*\approx 7$
, we observe the coalescence of mode B with the low-
$U^*$
branch of the previously merged modes A–C. These exceptional points are also of codimension-three and are plotted as 
in figure 9(b).
4.2.4. Effect of the spacing
Figure 11 maps the stability of the tandem of cylinders in the
$L-U^*$
plane for
${\textit{Re}} = 80$
and
$m^*=2.5$
. Increasing the distance between the bodies shows dynamics similar to what has been noticed for the fixed-tandem configuration by Zdravkovich (Reference Zdravkovich1987) or by Papaioannou et al. (Reference Papaioannou, Yue, Triantafyllou and Karniadakis2008) for the 2DOF tandem configuration. Mainly, the evolution of mode A, which is foremost a fluid mode, seems to follow different wake interference regimes. For
$1.5\lt L\lt 1.8$
, mode A is unstable in the low
$U^*$
range, as described in § 4.2.1. This coincides with the `slender body’ regime of the fixed tandem, where the free shear layer of the upstream body does not reattach to the downstream one. The vortex shedding comes from the detachment of the upstream body shear layers. For
$1.8\lt L\lt 3$
, mode A is stable and mode C becomes unstable. For the fixed tandem, the shear layers of the front body reattach to the rear one. The vortex street is only formed in the wake of the rear body. In the
$3\lt L\lt 4.5$
region, mode A is unstable. This region could correspond to the `intermittent-regime’ where the vortex street of the upstream body intermittently reattaches to the rear one. In this region, the eigenmodes were found to be highly sensitive to variations in the base flow within the gap. Since the computational mesh had to be adapted for each value of the length
$L$
in order to solve the forced problem, the resulting neutral curve exhibited some oscillations. To address this, a quadratic fit was applied to smooth the curve. As of
$L\gt 4.5$
, the dynamics correspond to the binary vortex street regime where vortex streets are being shed in the wake of both bodies. Up to
$L=4.5$
, the modes behave in the same way described for
$L=3$
and
${\textit{Re}}=100$
in § 3.2. Mainly, mode A induces no displacement of the bodies and modes B and C induce the displacement of both bodies, depending on
$U^*$
. For spacings above
$L=4.5$
, however, a change of dynamics is observed. For mode A, both cylinders display motions that are of comparable amplitudes. Mode B induces higher displacement of the front cylinder than of the rear. Finally, mode C induces higher displacement of the rear cylinder than of the front. Modes B and C start to exhibit a decoupling of the body’s motion, whereas mode A arises from the coupling of both cylinders’ motion. This phenomenon is confirmed for higher distances, as is explained in the next subsection.
Figure 11. Curves of neutral stability in the
$L-U^*$
plane for a low reduced mass of
$m^*=2.5$
at
${\textit{Re}}=80$
, computed with the impedance criterion. The colour codes for the different modes are similar to figure 6. Light grey indicates regions where one unstable mode exists and dark grey regions where two unstable modes exist.
A possible explanation for the disappearance of mode A in the range
$1.8 \lt L\lt 2.8$
might be that for this range of spacing, the second cylinder lies in the recirculation region of the first one, which corresponds to the region of absolute instability responsible for the instability. A similar explanation was given by Hosseini et al. (Reference Hosseini, Griffith and Leontini2021) in their investigation of the fixed tandem and three fixed-cylinder configurations for
${\textit{Re}}=200$
. Accordingly, vortex shedding from the first cylinder is suppressed and the global flow reorganises. In the present case, which involves elastically mounted cylinders, the nature of mode A below
$L=4.5$
is dominated by the fluid, so one might expect this scenario to remain essentially valid.
4.2.5. Mode decoupling at high spacing
Figure 12 shows results from LSA and from the impedance-based method for
${\textit{Re}} = 60$
and
$m^* = 2.5$
at
$L=10$
. Figures 12(a) and 12(b), respectively, show the growth rate and frequency with respect to
$U^*$
for different configurations. Calculations where both bodies are spring-mounted are represented with plain lines. Calculations where either the front or rear cylinder was fixed are, respectively, shown with circles and rectangles. As can be seen, mode B (
) of the full problem (both cylinders free to oscillate) is identical to the mode arising from the configuration where the rear cylinder is fixed (
). On the other hand, mode C (
) is almost identical to the mode arising from the configuration where the front cylinder is fixed (
). Finally, mode A (
) only follows the mode of the rear-fixed cylinder case (
) for low
$U^*$
. For
$U^*\gt 5.5$
, neither mode 
nor 
accounts for mode A’s dynamics.
Figure 12. Results from LSA of cases with both freely oscillating cylinders (
, 
and 
), front cylinder fixed (
, 
and 
) and rear cylinder fixed (
, 
and 
). Real and imaginary parts of the leading eigenvalues with respect to
$U^*$
at
${\textit{Re}} = 60$
and
$L=10$
for
$m^* = 2.5$
(a,b). Correspondingly, zero isolines of the real and imaginary parts of
$H_1$
(respectively, 
and 
) and
$H_2$
(respectively, 
and 
) in the
$\omega - U^*$
plane (c). The points of neutral stability of the different modes for both cylinders freely oscillating are reported in both figures as 
, 
and 
.
The coupled and decoupled nature of the different modes can also be seen by examination of the terms present in the impedance matrix. Let us define the diagonal terms of
$Z_T$
(2.40) as
\begin{equation} \begin{cases} H_1= - \omega ^2 - \dfrac {4\pi \gamma _1}{U_1^*} i \omega + \left ( \dfrac {2\pi }{U_1^*} \right ) ^2 - \dfrac {i \omega 4 Z_{1,1}}{\pi m_1^*}, \\[3ex] H_2 = -\omega ^2 - \dfrac {4\pi \gamma _2}{U_2^*} i \omega + \left ( \dfrac {2\pi }{U_2^*} \right ) ^2 - \dfrac {i \omega 4 Z_{2,2}}{\pi m_2^*}. \end{cases} \end{equation}
Figure 12(c) shows the zero isolines of the real and imaginary parts of
$H_1$
(respectively, 
and 
) and
$H_2$
(respectively, 
and 
) in the
$\omega - U^*$
plane. The points of neutral stability of the full linear problem (LSA) are reported in figure 12(c) as dashed circles. One can note that the zeros of
$H_1$
approximate relatively well points P1 and P2. As shown from the LSA results, modes B and A are linked to the rear-fixed configuration at these reduced velocities. Similarly, P5 and P6 are relatively well approximated by the zeros of
$H_2$
. Correspondingly, mode C is linked to the front-fixed configuration at these reduced velocities. On the other hand, the neutral points P3 and P4 are not detected by the sole zeros of
$H_1$
or
$H_2$
, showing there the coupled nature of mode A. Accordingly, the LSA results show that in that range of reduced velocities, mode A cannot be described from the front-fixed or rear-fixed configurations solely.
5. Conclusion and perspectives
The purpose of this paper was to give a thorough description of the linear stability properties of the flow past a tandem of elastically mounted cylinders, involving both a parametric analysis in a large range of parameters and a physical description of the unstable eigenmodes.
Due to the large number of parameters involved in the problem, one was led to design an efficient method to perform parametric analyses. For this sake, we designed a generic method which consists of solving in a first step a series of forced problems in which the cylinders are imposed to harmonically oscillate at a given real frequency
$\omega$
. This allows us to define impedance functions
$Z_{\textit{ij}}(\omega )$
, condensing in a synthetic way the retroaction of the fluid onto the oscillating bodies. This allows, in a second step, to predict the instability criteria by the sole inspection of a
$2\times 2$
matrix involving the parameters of the structure (reduced mass, velocity and damping) along with these impedance functions.
The impedance-based method was compared with the resolution of the eigenvalue problem for the coupled fluid–body problems, yielding very good agreement. Using the impedance-based stability criterion, an extensive parametric study was then conducted. The effect of increasing the mass and the damping ratio are found to be generally stabilising, except for mode A which is destabilised by the increase of the mass. The stability analysis in the
$L-U^*$
plane for
$ Re = 80$
and
$m^* = 2.5$
reveals that the tandem cylinder dynamics evolve with spacing, showing similarities to previously observed wake interference regimes. Mode A dominates in certain regions and aligns with fluid instabilities, while transitions between regimes reflect changes in wake interactions, from slender body behaviour to binary vortex shedding. For
$L \gt 4.5$
, a shift in dynamics occurs, with modes B and C showing decoupled motions of the cylinders, while mode A remains a coupled mode. This is confirmed by further analysis at
$ Re = 60$
and
$ L = 10$
, where comparison with fixed-body configurations shows that modes B and C align with the rear-fixed and front-fixed cases, respectively. Mode A, however, cannot be captured by either configuration alone at higher reduced velocities, underscoring its coupled nature, as supported by the impedance matrix analysis.
A number of perspectives are opened for the continuation of this study. First, the formalism developed here applies to an arbitrary number of rigid bodies, with almost identical computational cost compared with the fluid problem owing to the use of the discrete-ALE ansatz. An illustration for a system of three cylinders is given in the Appendix B. Extension to a higher number of bodies is straightforward. A particularly interesting perspective would be to consider periodic configurations, for which one can expect to apply Floquet–Bloch approaches allowing one to predict the large-scale dynamics of such configurations based on the resolution of problems simply formulated into an elementary cell.
Aside from increasing the number of bodies, extending to additional degrees of freedom, including streamwise oscillation and rotational motion, could also be considered. According to the discussion in the introduction, such degrees of freedom are not expected to lead to novel dynamics in the case of tandem bodies for symmetry reasons, as in-line oscillations would be coupled to symmetric vortex shedding which is unlikely to occur. On the other hand, streamwise motions could come into play for configurations such as the side-by-side cylinders or the three-cylinder in pinball configuration investigated in Appendix B. In such a case, streamwise oscillations in an antisymmetric manner could happen, especially in configurations where the cylinders are close enough to lead to a global antisymmetric vortex shedding similar to that of a single composite body (Carini, Giannetti & Auteri Reference Carini, Giannetti and Auteri2014).
Another perspective would be to investigate the dynamics of a set of oscillating cylinders located below a free surface. Recent studies (Patel et al. Reference Patel, Sun, Yang and Zhu2025) have investigated the effect of a free surface on the vortex-shedding instability behind a fixed cylinder. Interesting dynamics can be anticipated by considering spring-mounted objects. This topic is currently the object of ongoing investigations in our team.
Lastly, since the present study was restricted to the linear regime, one should extend the study by considering nonlinear dynamics. Considering the nonlinear dynamics is especially interesting for two main reasons. First, in the course of the parametric study, a number of exceptional points of codimension-two or -three have been identified. Rich dynamics can be expected in the vicinity of such points, which could be elucidated thanks to nonlinear simulations or alternative approaches such as weakly nonlinear expansions, time spectral method or spectral submanifolds. Secondly, considering the application to energy harvesting which was the original motivation for this work, evaluation of the performances and energy efficiencies necessarily requires investigation of the nonlinear regimes. A key question for such applications will be to identify the optimal range of parameters, reduced mass and velocity, and most importantly, the parameter
$\gamma _i$
identified here with a damping parameter, but which for an energy harvester would rather be identified as an energy extraction coefficient.
Acknowledgements
T.M. would like to thank E. Boujo for sharing unpublished results and exchanging thoughts during the validation process of the linear ALE method presented in this paper.
Data availability
The codes used to produce the results in this paper are available in the StabFEM open-source project.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Exchange of stability between mode branches at
$m^*=20$
In § 4.2 we showed that for
$m^*=20$
, multiple exchanges of stability between the modes occur. To clarify the dynamics taking place, we show the growth rates, frequencies as well as root loci from of the different modes from LSA (figure 14) for different relevant Reynolds numbers (plotted as 
figure 13, which is a close-up of figure 6
b). The points of neutral stability have been plotted figure 14 as 
or 
, according to the colours of the neutral curves in figure 13.
Figure 13. Curves of neutral stability in the
${\textit{Re}}-U^*$
plane for
$m^*=2.5$
and
$L=1.5$
, from the impedance-based method (zoom of figure 6
b). The exceptional points are plotted as 
.
Figure 14. Growth rates (a,d,g,j), frequencies (b,e,h,k) and root loci (c,f,i,l) of the leading eigenmodes (LSA) for
$m^*=20$
and
$L=1.5$
:
${\textit{Re}}=100$
(a–c),
${\textit{Re}}=78$
(d–f),
${\textit{Re}}=75$
(g–i) and
${\textit{Re}}=60$
(j–l). The natural frequency of a spring-mounted cylinder in vacuum
$\omega _n={2\pi }/{U_n^*}$
is shown as 
. The growth rate and frequency of the fluid mode in the wake of the corresponding fixed configuration are displayed as 
.
When decreasing the Reynolds number, the first exchange occurs around
${\textit{Re}}\approx 84$
. The A branch for small
$U^*$
becomes connected with what was previously the B branch and vice versa, as can be seen comparing figure 14(a) for
${\textit{Re}}=100$
and figure 14(d) for
${\textit{Re}}=78$
. This exceptional point occurs for
$U^*\approx 6.8$
and is identified figure 13 as 
.
The second exchange appears in the loop feature of the neutral curve, around
${\textit{Re}}\approx 75.5$
. For
$U^*\approx 9.75$
(second symbol 
in figure 13), the two unstable branches connect and exchange stability behaviour as can be seen comparing figures 14(d) and 14(g). The asymptotic limit for low
$U^*$
at
${\textit{Re}}=75$
corresponds to the critical Reynolds number for which the mode of the fixed tandem configuration is neutral (see 
figure 14
g).
Note that all neutral points plotted 
as occur when the modes are structure-dominated, as can be seen by looking at the matching frequency. On the other hand, all neutral points plotted as 
arise from fluid-dominated modes.
Appendix B. Impedance-based criterion of a three-body system
The impedance-based stability prediction can be applied to any number
$n$
of bodies. We show here the validity of the method for the prediction of the stability of a cluster of three cylinders that are centred on the vertices of an equilateral triangle of side length
$3D/2$
. The cylinders are free to oscillate in the transverse direction of the flow. For fixed cylinders, Chen et al. (Reference Chen, Ji, Alam, Williams and Xu2020) showed that the flow dynamics were highly sensitive to the spacing
$L$
and Reynolds number
${\textit{Re}}$
. This configuration is also known as the fluidic pinball when the cylinders are rotatable. The LSA at
${\textit{Re}}=60$
for
$m^*=2.5$
introduces five leading eigenmodes, two being fluid-dominated modes while the three others are of structure-dominated nature. Figure 15 shows the real and imaginary parts of the leading eigenvalues against
$U^*$
(respectively, figures 15
a and 15
b) as well as the vorticity field of the different modes at a reduced velocity of
$U^*=7$
(figure 15
c–g). Modes A (
) and B (
) are both unstable over the whole range of reduced velocities investigated and their frequencies, respectively, match the frequencies of the two unstable modes found in the corresponding fixed case (displayed 
as and 
). These modes do not induce displacement of any of the three bodies. Their vorticity fields are, respectively, plotted in figures 15(d) and 15(f) at reduced velocity
$U^*=7$
. Modes C (
), D (
) and E (
), respectively, become unstable at reduced velocities
$U^*=5.88$
,
$U^*=5.13$
and
$U^*=7.67$
. All three modes follow the natural frequency of the structure-only spring-mounted system (
$\omega _n= {2\pi }/{U_n^*}$
, shown as 
). Modes C and E induce the displacement of all cylinders. Both eigenmodes are symmetric and the zone of high vorticity is localised around the bodies as well as in the outer regions of the wake in proximity to the bodies, as seen in figures 15(c) and 15(g). Mode D, on the other hand, induces the displacement of both rear cylinders, leaving the front one stationary. The structure of the mode is similar to that of the other structure-dominated modes but it is antisymmetric, as can be seen in figure 15(e). We computed the forced problem for that configuration and generalised the impedance-based criterion described in § 2.4 to a three-body problem. Practically, the thresholds are found by the inspection of the determinant of a
$3\times 3$
matrix. The stability predictions from the impedance-based criterion are plotted as 
in figures 15(a) and 15(b). The model correctly predicts the reduced velocities at which the growth rates of the modes become neutral. The prediction of the modes’ frequencies is also in very good agreement with the results from LSA.
Figure 15. Real part of the leading eigenvalues from LSA (a) and their frequencies (b) with respect to
$U^*$
at
${\textit{Re}} = 60$
for
$m_1^* = m_2^* = m_3^* = 2.5$
. The unstable region is depicted as the grey zone. The natural frequency of a spring mounted cylinder in vacuum
$\omega _n=({2\pi }/{U_n^*})$
is shown as 
. The two leading eigenvalues of the corresponding fixed tandem case are displayed as 
and 
. The predictions from the impedance criterion are shown as 
. Vorticity fields of the eigenmodes (c–g) are showed for
$U^*=7$
.
Open access
