On the development of a nonlinear time-domain numerical method for describing vortex-induced vibration and wake interference of two cylinders using experimental results

Abstract

A nonlinear mathematical model is developed in the time domain to simulate the behaviour of two identical flexibly mounted cylinders in tandem while undergoing vortex-induced vibration (VIV). Subsequently, the model is validated and modified against experimental results. Placing an array of bluff bodies in proximity frequently happens in different engineering fields. Chimney stacks, power transmission lines and oil production risers are few engineering structures that may be impacted by VIV. The coinciding of the vibration frequency with the structure natural frequency could have destructive consequences. The main objective of this study is to provide a symplectic and reliable model capable of capturing the wake interference phenomenon. This study shows the influence of the leading cylinder on the trailing body and attempts to capture the change in added mass and damping coefficients due to the upstream wake. The model is using two coupled equations to simulate the structural response and hydrodynamic force in each of cross-flow and stream-wise directions. Thus, four equations describe the fluid–structure interaction of each cylinder. A Duffing equation describes the structural motion, and the van der Pol wake oscillator defines the hydrodynamic force. The system of equations is solved analytically. Two modification terms are added to the excitation side of the Duffing equation to adjust the hydrodynamic force and incorporate the effect of upstream wake on the trailing cylinder. Both terms are functions of upstream shedding frequency (Strouhal number). Additionally, the added mass modification coefficient is a function of structural acceleration and the damping modification coefficient is a function of velocity. The modification coefficients values are determined by curve fitting to the difference between upstream and downstream wake forces, obtained from experiments. The damping modification coefficient is determined by optimizing the model against the same set of experiments. Values of the coefficients at seven different spacings are used to define a universal function of spacing for each modification coefficient so that they can be obtained for any given distance between two cylinders. The model is capable of capturing lock-in range and maximum amplitude.

1 Introduction

An array of bluff bodies, placed in proximity, is a frequent set-up in various fields of engineering such as production risers along a side of an FPSO or tendons of a Tensioned Leg Platform (TLP) in offshore engineering, as shown in Fig. 1. The dynamics of a body immersed in the wake of another structure is significantly different from the same body when it is in an undisturbed flow. Based on the orientation of structures concerning flow and each other, three configurations are possible: side by side, tandem and staggered (Fig. 2).

Fig. 1

Interference of TLP legs in proximity under action of waves and ocean currents

Fig. 2

Three possible arrangements of a pair of cylinders in close proximity. a Side by side, b tandem and c staggered

Zdravkovich [1] conducted extensive experiments on a pair of cylinders at different configurations and observed three regions in upstream wake based on its interference with trailing cylinder: “Proximity Interference” at distances less than 1.2D to 1.8D where the pair behave as a single body; “Wake Interference” in which trailing cylinder is fully or partially submerged in leading cylinder’s wake. “No-interference” is the third region where two cylinders are placed far away from each other enough to behave as two single bodies.

Numerous experiments have been conducted on two cylinders in proximity interference and wake interference regions. They mainly focus on the study of flow regime around structures and hydrodynamic coefficients (lift and drag). Sumner et al. [2] offered an extensive classification of flow pattern around a pair of cylinders in staggered arrangement. They observed that when two cylinders were attached to each other or placed at a very small distance, only one vortex street was formed and they acted as a single bluff body. Based on their observation, when the angle between two cylinders is $\psi <30^{\circ }$, the flow pattern could be divided into three groups: (i) at small pitch ratio and small angle of incident, the upstream shear layers reattached on the trailing cylinder; (ii) as $\psi$ grew larger, the reattachment could not be maintained so the shear layer was deflected into and rolled up in the gap between two cylinders which induced separation on the trailing cylinder; and (iii) while $\psi$ was still small, if the gap grew larger, the deflected shear layers in the gap could form a fully developed Kármán vortex street which was referred to as vortex impingement flow pattern. Furthermore, in arrangements with a large angle of incident, both cylinders developed separate vortex streets. The most common flow pattern was synchronised vortex shedding where vortex streets were synchronised and formed adjacent anti-phase streets. In the same region, Zhang et al. [3] conducted their experiment in which the leading cylinder was allowed to oscillate in cross-flow direction in front of a fixed counterpart. They observed that the leading cylinder underwent galloping at $0.3<L/D<1.2$ where due to lack of damping the oscillation amplitude continued to grow, unlike a typical VIV response. The oscillation amplitudes dropped dramatically at higher spacings $1.5<L/D$.

The alteration of flow pattern around two cylinders in wake interference region has a significant influence on the pressure gradient around the cylinders. Igarashi [4] observed that the base pressure of leading cylinder is proportional to the spacing. He argued that the shear layers had enough time to form vortices in larger spacings. As the first vortex was formed, pressure distribution was similar to a single cylinder. Such an observation was confirmed in studies where cylinders were allowed to oscillate; for instance, Armin et al. [5] observed that leading cylinder underwent vortex-induced vibration (VIV) response similar to an isolated cylinder at $L/D>4$ (where L is the centre-to-centre distance between the cylinders and D is the diameter, as shown in Fig. 2). On the contrary, Kim et al. [6] and Huera-Huarte and Gharib [7] observed dramatic variation in response of both cylinders where $L/D<4$. These different observations suggest that the flow regime changes at a spacing between $L/D=3$ to 4 which is referred to as critical spacing. The exact value of critical spacing depends on Re. The flow regime transformation also impacts hydrodynamic coefficients and Strouhal number [8].

The hydrodynamic coefficients meaningfully depend on the location of two bodies relative to each other. Zdravkovich [1] provided a comprehensive map of hydrodynamic coefficients and Strouhal number of the cylinders at different positions. Sumner et al. [9] reported a sudden jump in drag and lift coefficients as well as Strouhal number when leading cylinder forming an independent vortex street. Such a behaviour has been observed in many other studies, [4, 8, 10, 11]. These observations all confirm that hydrodynamic coefficients become a function of spacings (L/D) as well as Re.

As a result of alteration in hydrodynamic coefficients, the response of trailing cylinder is significantly different from a typical VIV response of an isolated cylinder. Assi et al. [12] conducted their experiment at “Wake Interference” region and placed a flexibly mounted cylinder in the wake of a stationary cylinder. They observed that at small spacings, trailing cylinder response is similar to that of an isolated cylinder. However, contrary to the VIV response of an isolated cylinder, the amplitude did not decrease at high velocities but displayed a galloping-like response. Thus, they recognized two different motions for trailing cylinder; one was vortex-induced vibration which was excited by vortices detached from the cylinder itself, and the other was wake-induced vibration (WIV) caused by buffeting vortices in the wake of leading cylinder. Moreover, they divided the oscillation response into three regions based on flow velocity: region of solely VIV response at low velocities, a region of solely WIV motion at very high velocities and a mid-region where the response was excited by the combination of VIV and WIV.

Armin et al. [5] considered a more general experiment set-up and focused on the amplitude and frequency response of two oscillating tandem cylinders rather than flow regime around them. They observed how flow velocity and spacing affected the response amplitude of both cylinders. It was concluded that the trailing cylinder response amplitude was not only a function of the undisturbed flow velocity but also the distance between two cylinders. Spacing was observed to influence the cylinder response in two different manners. It determined the flow velocity that the trailing cylinder was experiencing (shielding effect), as well as a secondary force, exerted to the structure by the buffeting upstream vortices.

The interaction of two cylinders in proximity is explored extensively, and despite its complexity, it is relatively understood. Nevertheless, attempts to offer a model for simulating this interaction are limited and almost non-existence. Having a time-domain model that can provide fast VIV and WIV simulations with good accuracy is important at initial steps of design. There are mathematical models simulating VIV, but they are limited to an isolated cylinder and do not capture wake interference.

These models typically utilize two differential equations to describe the structural response and hydrodynamic forces and couple them together to represent the fluid–structure interaction [13]. A common approach in the literature is to simulate the structural response by a simple equation of motion with the wake force as the excitation term. Additionally, the self-exciting and self-limiting nature of the fluid force could be simulated by a wake oscillator. The wake oscillator model is often represented by a van der Pol or Rayleigh equation and is related to the equation of motion by a coupling term. These are known as empirical models due to the inclusion of empirical coefficients. Empirical coefficients are usually determined by model calibration against experimental results.

Facchinetti et al. [14] conducted an extensive study on simulation of the fluid force on a rigid isolated cylinder by a van der Pol wake oscillator and simulated structure response with a mass–damping system. Moreover, they investigated different coupling terms. Three different coupling terms were considered which were proportional to cylinder displacement, speed and acceleration. It was concluded that acceleration coupling yielded the most agreeable results with experiment.

However, their study was limited to a constrained system with one degree of freedom (DoF) in cross-flow, whereas in engineering applications, structures generally have higher DoF. Thus, Zanganeh [15] tried to extend this model to a system oscillating in stream-wise and cross-flow. In this study, he suggested replacing the structure motion equation with a Duffing oscillator. He demonstrated that empirical coefficients could be determined as a function of the mass and/or damping to omit the need for calibration against experimental results.

Furthermore, different approaches were adopted to provide alternative time-domain models. Thorsen et al. [16] simulated the wake force by a modified form of Morison’s equation. The modification was done so that the drag term could simulate a controlling effect on vibration. The force obtained from this equation was then implemented to a finite element analysis software as an input to simulate the structural response. Bai and Qin [17] used Rayleigh equation to capture the wake force. However, rather than a simple linear coupling term, energy generated by wake was considered as the excitation term. Potential flow was used to simulate the excitation term in the Rayleigh equation. They divided the wake force into two regions, one close to the cylinder wall where vortices were generated and other further downstream where vortices were detached and simulated as discrete point vortices. They were able to simulate the vibration amplitude in both directions, frequency and trajectory of motion. Similarly, wake force was divided into two components by Skop and Balasubramanian [18], an excitation term and a stall parameter. They suggested a van der Pol equation to simulate wake excitation force. The proposed stall parameter was defined as a function of shedding frequency and cylinder velocity that provided negative values for large structural motion and could couple the wake force to the structural motion equation.

Nonetheless, these attempts are limited to an isolated cylinder. The current study aims to modify the mathematical model developed by Srinil and Zanganeh [19] so that it captures interference between two cylinders in tandem. One of the objectives in this study is to maintain the simplicity of wake oscillator models and avoid focusing on complex and varying interaction between two cylinders for each spacing. The new model will be able to simulate the trailing cylinder response due to VIV and WIV altogether. In this regard, the first few assumptions are made to obtain an initial model similar to Shiau [20]. Then, results from the initial model are compared with experiments done by Armin et al. [5]. Moreover, these experimental results are used as a benchmark to improve the initial model. Finally, an overall model is proposed which is able to capture the onset of lock-in, maximum amplitude and lock-in range width.

2 Modelling

To simulate the interaction between a pair of cylinders, it is assumed that cylinders are identical, meaning both have the same dimensions and structural properties. Each cylinder was modelled by a simple mass–spring and damper system similar to Fig. 3.

Fig. 3

Model of oscillating cylinders in tandem as simple mass, spring and damping systems

2.1 A wake oscillator to describe leading cylinder

If the stream direction is assumed to be from top to down, then F(t) is the lift force which induces motion response, Y(t). According to such a system, the structural response can be modelled by an equation of motion [Eq. (1)]. Dotted parameters in this equation and throughout this paper represent derivatives with respect to time. Lift (wake) force exerted on the structure is proportional to flow velocity and oscillating lift coefficient $(C_L)$ of the cylinder and can be obtained by Eq. (2).

$$\begin{aligned}&M{\ddot{Y}}+c{\dot{Y}}+kY=F_{Y}(t) x \end{aligned}$$

(1)

$$\begin{aligned}&F_{Y}=\frac{1}{2}\rho U^{2}DC_{L}. \end{aligned}$$

(2)

Here, c is the sum of viscous and fluid added damping $(c = c_\mathrm{s} + c_\mathrm{a})$ where $c_\mathrm{a}$ can be calculated using Eq. (3). ($\omega _\mathrm{s}$ is vortex shedding period, and $\gamma$ is stall parameter which is a function of mean drag coefficient [21].)

$$\begin{aligned} c_\mathrm{a}=\gamma \omega _\mathrm{s}\rho D^{2}. \end{aligned}$$

(3)

Also, mass (M) is the combination of structural mass (m) and fluid added mass $(m_\mathrm{a})$ which can be calculated from the following expression.

$$\begin{aligned} m_\mathrm{a}=\frac{\pi C_\mathrm{a} \rho D^{2}}{4}. \end{aligned}$$

(4)

Here, $C_\mathrm{a}$ is fluid added mass coefficient which is considered as unit for a smooth cylinder [21]. It is necessary to use non-dimensional parameters in mathematical modelling, so the model could be used irrespective of structure dimensions. Mass ratio is a non-dimensional parameter to represent the total mass which is defined as:

$$\begin{aligned} \mu =\frac{m+m_\mathrm{a}}{\rho D^{2}}. \end{aligned}$$

(5)

Simulation of fluid interaction with structure has been discussed in literature extensively [22]. It has been remarked in these studies that VIV response is a self-exciting and self-limiting phenomenon, and therefore, van der Pol equation [Eq. (6)] has been suggested for simulation of the oscillating lift coefficient $(C_L)$.

$$\begin{aligned} \ddot{C}_{L}+\epsilon \omega _\mathrm{s}\left( {C_{L}}^{2}-1\right) \dot{C}_{L}+{\omega _\mathrm{s}}^{2}C_{L}=T. \end{aligned}$$

(6)

$\epsilon$ is an empirical coefficient that should be determined case by case against experimental results.

Some studies, [14], replaced $C_{L}$ with q which is reduced vortex lift coefficient and is equal to twice the oscillating lift coefficient over a stationary cylinder lift coefficient $(C_{L_{0}})$. The reference value of $C_{L_{0}}$ can be considered 0.3 for a wide range of Re based on Blevins [21] and Pantazopoulos [23] studies.

T on the right-hand side of the wake oscillator is the coupling term. This term is defined to describe interaction between fluid and structure. Based on Facchinetti et al. [14] study on dynamic coupling terms, the term related to acceleration has the best agreement with experimental results. Therefore, a simple linear function of acceleration $(A{\ddot{Y}}_{1})$ is considered. Equation (7) can describe the structural vibration of leading cylinder in cross-flow.

$$\begin{aligned} {\left\{ \begin{array}{ll} M{\ddot{Y}}_{1}{+}(2\xi M\omega _{n}{+}\gamma \omega _\mathrm{s}\rho D^{2}){\dot{Y}}_{1}{+}kY_{1}{=}\frac{1}{2}\rho U^{2}DC_{L} \\ {\ddot{C}}_{{L}_{1}}{+}\epsilon \omega _\mathrm{s}\left( {C_{L_1}}^2{-}1\right) {\dot{C}}_{{L}_{1}}{+}{\omega _\mathrm{s}}^{2}{C}_{{L}_{1}}{=}A{\ddot{Y}}_{1} \\ \end{array}\right. } \end{aligned}$$

(7)

$\gamma$ can be assumed constant and equal to 0.8 in sub-critical region $(300<\mathrm{Re}<1.5\times 10^{5})$ for the sake of simplicity.

As mentioned before, to apply this model to any set-up regardless of structural dimensions, it is necessary for Eq. (7) to be in a dimensionless form. This is possible by introducing dimensionless time and space coordinates, $\tau =\omega _{n}t$ and $y =Y/D$, respectively. By replacing these dimensionless variables into Eq. (7), it becomes:

$$\begin{aligned} {\left\{ \begin{array}{ll} \ddot{y}_{1}+\left( 2\xi +\dfrac{\gamma }{\mu }\omega _{0}\right) {\dot{y}}_{1}+y_{1}=a{C}_{L_{1}}\\ {\ddot{C}}_{L_{1}}+\epsilon \omega _{0}\left( {{C}_{L_{1}}}^{2}-1\right) {\dot{C}}_{L_{1}}+{\omega _{0}}^{2}{C}_{L_{1}}=A\ddot{y}_{1}\\ \end{array}\right. }.\nonumber \\ \end{aligned}$$

(8)

Here, A is another empirical coefficient that should be determined by tuning the model against appropriate data. $\omega _0$ is the ratio between the vortex shedding period and the system natural period, $(\omega _0=\omega _\mathrm{s}/\omega _n )$ and $a=\dfrac{1}{8}\dfrac{\omega _{0}^{2}}{\pi ^{2}\mathrm{St}^{2}\mu }$. The complete mathematical steps to drive the dimensionless formula for both wake oscillators can be found in Armin [24]. Strouhal number (St) which is a dimensionless number and a function of vortex shedding frequency $(f_\mathrm{s})$ and free stream velocity is defined as $\mathrm{St}=\dfrac{f_{\mathrm{s}}D}{U}$. Strouhal number is also a function of Re and roughness and is assumed to be equal to 0.2 for a sub-critical range of Re [23]. This system should be solved simultaneously to simulate the response of a rigid cylinder with one degree of freedom. Velocity should also be stated in a non-dimensional form which is referred to as reduced velocity and is represented by $U_\mathrm{r} (=U/f_n D$, where $f_n$ is structure natural frequency).

2.2 A wake oscillator to describe trailing cylinder

Modelling two cylinders in tandem requires considering a similar system in the wake of the first cylinder plus the interaction between them. Two mechanisms of excitation were observed for the trailing cylinder response by Armin et al. [5]. The trailing cylinder response was observed to be induced by vortices from the cylinder itself and the buffeting vortices detached from leading cylinder. The response to vortex detachment from cylinder’s aft can be modelled by wake oscillators discussed earlier. Additionally, modifying this model to capture the effect of the buffeting upstream vortices is done by adding a force term $(P_{Y_{1}}(t))$ to excitation side of the structural motion equation [Eq. (9)].

$$\begin{aligned} M{\ddot{Y}}_{2}+c{\dot{Y}}_{2}+kY_{2}=F_{Y_{2}}(t)+P_{Y_{1}}(t). \end{aligned}$$

(9)

Shiau et al. [25] suggested to assume vortices convey the same energy to trailing cylinder as they do to the leading one during detachment. Therefore, they replaced $P_{Y_1 }$ with the wake force obtained from Eq. (8) plus a time delay to consider the time $(t_1)$ required by upstream vortices to reach trailing cylinder. Moreover, the time delay was defined as a function of the spacing between two cylinders (L) , spacing between vortices (d) and shedding frequency.

Solving the system of nonlinear differential equations discussed here is possible by making a few assumptions about response functions. Since response of a cylinder undergoing VIV is sinusoidal, it is valid to assume that the response functions have amplitudes of $y_1$ and $y_2$ and periods $\omega _1$ and $\omega _2$ [26] for leading and trailing cylinders, respectively. Furthermore, a force inducing a sinusoidal motion should be sinusoidal with the same frequency. Hence, $C_{L_1}$ and $C_{L_2}$, as the excitation forces, should have similar solutions with phase differences to motion amplitudes [Eq. (10)].

$$\begin{aligned} {\left\{ \begin{array}{ll} y_{1}={\bar{y}}_{1}e^{i\omega _{1}t}\\ C_{L_{1}}={\bar{C}}_{L_{1}}e^{i\omega _{1}t+\phi _{1}}\\ y_{2}={\bar{y}}_{2}e^{i\omega _{1}t+\theta _{2}}\\ C_{L_{2}}={\bar{C}}_{L_{2}}e^{i\omega _{1}t+\phi _{2}}.\\ \end{array}\right. } \end{aligned}$$

(10)

Armin et al. [5] observed that leading cylinder dictates the oscillation response of both cylinders up to high reduced velocities, and additionally, Okajima [27] and Tsutsui [28] observed that both cylinders have identical St where they were fixed; therefore, it is a valid assumption, for the sake of simplicity, that both cylinders are oscillating with similar frequencies.

$\phi _1$and $\phi _2$ are phase differences between leading cylinder motion and its wake force and trailing cylinders wake, respectively. Moreover, $\theta$ is the phase difference between leading and trailing cylinders motion.

3 2-DoF model

The relationship between stream-wise and cross-flow motion has been discussed previously in the literature [29]. When a structure is flexibly mounted and allowed to oscillate in both directions, a relative velocity appears between flow velocity and the oscillating structure, as shown in Fig. 4a. The direction of the fluid force acting on an oscillating cylinder rotates clockwise (Fig. 4a, b) or counterclockwise due to the relative motion of the cylinder with respect to flow. In other words, drag force $(F_{D})$ is not along the stream direction but creates an angle of $\beta$ which is time-dependent and is a function of the cylinder instantaneous velocity. Additionally, lift force $(F_L )$ is always perpendicular to the drag. Armin et al. [5] observed the trajectory crescents of both cylinders pointing downstream; therefore, stream-wise and cross-flow forces could be resolved (for counterclockwise rotation) in drag and lift directions as:

$$\begin{aligned} {\left\{ \begin{array}{ll} F_{X}=F_{D}\cos \beta - F_{L}\sin \beta \\ F_{Y}=F_{D}\sin \beta + F_{L}\cos \beta .\\ \end{array}\right. } \end{aligned}$$

(11)

Fig. 4

a Relative velocity and b force outcome for an oscillating cylinder with 2-DOF

By assuming that $\beta$ is small, it can be defined as $\sin \beta =-\frac{{\dot{Y}}}{U}$. Following the steps explained by Blevins and Coughran [30], Eq. (12) is obtained.

$$\begin{aligned} {\left\{ \begin{array}{ll} F_{X}=F_{D}+ F_{L}\dfrac{{\dot{Y}}}{U} \\ F_{Y}=F_{L}-F_{D}\dfrac{{\dot{Y}}}{U}.\\ \end{array}\right. } \end{aligned}$$

(12)

Srinil et al. [31] focused on this issue with a pendulum set-up for an isolated cylinder test. They suggested that due to geometry nonlinearity of the spring–mass system, structural motion equations should be in the form of a Duffing oscillator, Eq. (13), [32]. Two terms of $(x^3;\;y^3 )$ capture the axial nonlinear properties, and $(xy^2;\;yx^2)$ represent physical coupling between cross-flow and stream-wise motions. They also referred to Jian-Shu et al. [33] and Paul Raj and Rajasekar [34] as two other applications of such a coupled system.

$$\begin{aligned} {\left\{ \begin{array}{ll} \ddot{x}_{1}+\left( 2\xi +\dfrac{2\gamma }{\mu }\omega _{0}\right) {\dot{x}}_{1}+\left( x_{1}+\alpha _{x}{x_{1}}^{3}+\beta _{x}x_{1}{y_{1}}^{2}\right) \\ \quad =4aC_{D_{1}}+2\pi aC_{L_{1}}\dfrac{\dot{y}_{1}}{U_\mathrm{r}}\\ {\ddot{C}}_{D_{1}}+2\epsilon \omega _{0}\left( {C_{D_{1}}}^{2}-1\right) {\dot{C}}_{D_{1}}+4{\omega _{0}}^{2}C_{D_{1}}=A\ddot{x}_{1}\\ \ddot{y}_{1}+\left( 2\xi +\dfrac{\gamma }{\mu }\omega _{0}\right) {\dot{y}}_{1}+\left( y_{1}+\alpha _{y}{y_{1}}^{3}+\beta _{y}y_{1}{x_{1}}^{2}\right) \\ \quad =aC_{L_{1}}-8\pi aC_{D_{1}}\dfrac{\dot{y}_{1}}{U_\mathrm{r}}\\ {\ddot{C}}_{L_{1}}+\epsilon \omega _{0}\left( {C_{L_{1}}}^{2}-1\right) {\dot{C}}_{L_{1}}+{\omega _{0}}^{2}C_{L_{1}}=A\ddot{y}_{1}.\\ \nonumber \end{array}\right. }\!\!\!\!\!\!\\ \end{aligned}$$

(13)

Coefficients $\alpha _x,\;\beta _x,\;\alpha _y$ and $\beta _y$ are empirical coefficients which are determined by tuning against experimental data. In this study, these coefficients are assumed to be identical and equal to 0.7 [31].

4 Hydrodynamic force of the upstream wake

Isolated cylinder model requires significant modification to capture the galloping-like response of trailing cylinder. Spacing effect is not significant in proposed model by Shiau et al. [25]. Such an observation confirms that simply accounting and acknowledging the force of detaching vortices are not sufficient; the loss of vortices energy due to viscous resistance should be considered as well which results in a complex model.

4.1 Upstream wake influence

Any input from the upstream model adds mathematical complication to the model. The modification should be applied to eliminate any direct input from the leading cylinder and necessity to measure the vortices energy loss as they travel downstream. Furthermore, knowing the position of the trailing cylinder on its trajectory at the time of collision is important to determine the damping or excitation effect of each vortex. Thus, eliminating inputs from upstream cylinder mitigates these constraints and simplifies the model.

The upstream model can be used as a benchmark to modify the downstream model. The mathematical model can simulate the leading cylinder behaviour accurately. Thus, it requires an amendatory term which can simulate the difference between trailing and leading cylinders hydrodynamic forces ($F_1$ and $F_2$, respectively).

Wake force exerted on each cylinder can be calculated through Eq. (14) from experimental data obtained from Armin et al. [5].

$$\begin{aligned} {\left\{ \begin{array}{ll} F_{X}=\left( m_\mathrm{s}+m_\mathrm{a}\right) {\ddot{X}}+2\xi \omega _{n}\left( m_\mathrm{s}+m_\mathrm{a}\right) {\dot{X}}+kX\\ F_{Y}=\left( m_\mathrm{s}+m_\mathrm{a}\right) {\ddot{Y}}+2\xi \omega _{n}\left( m_\mathrm{s}+m_\mathrm{a}\right) {\dot{Y}}+kY.\\ \end{array}\right. } \end{aligned}$$

(14)

Velocity and acceleration of each cylinder are calculated by differentiating the corresponding displacement with respect to time. Then, solving Eq. (14) for $F_{D}$ and $F_L$ yields Eq. (15) through which oscillating drag and lift coefficients can be calculated.

$$\begin{aligned} {\left\{ \begin{array}{ll} C_{D}=\dfrac{F_{Y}\sin \beta -F_{X}\cos \beta }{\dfrac{1}{2} \rho _{w}DU^{2}L_{c}}\\ C_{L}=\dfrac{F_{Y}\cos \beta +F_{X}\sin \beta }{\dfrac{1}{2} \rho _{w}DU^{2}L_{c}}.\\ \end{array}\right. } \end{aligned}$$

(15)

Figure 5 demonstrates drag and lift coefficients at the corresponding centre-to-centre distance between two cylinders. These results are validated against drag force obtained through Vandiver expression [35]. The sum of point-by-point errors between the two sets of results is $14\%$.

Fig. 5

Oscillating drag and lift force exerted on leading (a, c, respectively) and trailing (b, d, respectively) cylinders by passing fluid versus drag force obtained from Vandiver equation

Figure 5a, b confirms that the increase in oscillation amplitude amplifies the oscillating drag coefficient for both cylinders [21, 23]. On the other hand, the drag amplification is not significant for the trailing cylinder. Contrary to previous observations for an isolated cylinder, the increase in trailing cylinder response amplitude does not result in drag and lift coefficients magnification and they are relatively constant. Both coefficients of trailing cylinder experience a sharp increase at the end of leading cylinder synchronization $(U_\mathrm{r}\approx 10)$ which is more significant in smaller spacings and non-existence for larger spacings $(L/D = 15$ to 20). This jump is due to change in excitation mechanism when the leading cylinder wake does not dictate the trailing cylinder motion which was discussed extensively by Armin et al. [5]. Trailing cylinder drag and lift display a direct dependency to spacing during upstream synchronisation. Moreover, the drag coefficient stays relatively constant at high reduced velocities, while the lift coefficient exhibits more sensitivity to the changes in spacing with an inverse correlation.

4.2 Modification coefficients

It is accepted in the literature that two mechanisms of excitation govern the trailing cylinder response, firstly, VIV motion due to the fluid current in the gap and secondly, buffeting vortices and turbulent flow regime in the wake of the leading body. Hence, the model requires a modification term that takes into account the effect of chaotic flow regime in the upstream wake, since the current model can simulate the VIV motion of the cylinders.

A simple method to identify the influence of the upstream wake is to deduct downstream wake force from upstream $(C_{D_2}-C_{D_1},~ C_{L_2}- C_{L_1})$. It should be emphasised that due to shielding effect, trailing body experiences a lower flow velocity which will be considered later through an additional modification coefficient.

Wake force could be divided into three different components, mean drag, oscillating lift and drag. Figure 6 displays oscillating drag and lift components for the spacing of $L/D = 4$. It is possible to examine this difference against several variables and find what parameter describes it the best. To avoid algebraic loops, the modification term should be considered simple so that it is readily transferable to the left-hand side, which means only a first-order polynomial function can be considered. Thus, three different variables were examined, leading cylinder cross-flow amplitude, trailing cylinder acceleration and velocity, see Fig. 6. It is evident that upstream displacement provides the best fit, however, as discussed before any input from the leading cylinder adds mathematical complications. Trailing cylinder acceleration provides satisfactory results as well.

Fig. 6

Curve fitting to hydrodynamic force (a drag, b lift) variance between up and downstream cylinder at spacing of $L/D=4$ with a polynomial function

4.2.1 Added mass modification term

Based on the previous discussion, the modification terms are defined as functions of acceleration. However, simple expressions of $A_{X}{\ddot{X}}$ and $B_{Y}{\ddot{Y}}$ would not yield non-dimensional terms after applying dimensionless time and distance. Therefore, the modification terms should be defined $A_{X}\dfrac{{\ddot{X}}_{2}}{D{\omega _\mathrm{s}}^{2}}$ and $B_{Y}\dfrac{{\ddot{Y}}_{2}}{D{\omega _\mathrm{s}}^{2}}$. ($A_{X}$ and $B_{Y}$ coefficients are functions of spacing to be determined using experimental results.) Hence, motion equations, with dimensionless time and space, become:¹

$$\begin{aligned} {\left\{ \begin{array}{ll} {\omega _{n}}^{2}D\ddot{x}_{2}+\left( 2\xi {\omega _{n}}^{2}D+\dfrac{2\gamma \omega _\mathrm{s}\omega _{n}\rho D^{3}}{M}\right) {\dot{x}}_{2}\\ +{\omega _{n}}^{2}Dx_{2}=\dfrac{\frac{1}{2}\rho U^{2}DC_{D}}{M}+\dfrac{\frac{1}{2}\rho U^{2}D}{M}A_{X}\dfrac{{\omega _{n}}^{2}D\ddot{x}_{2}}{D{\omega _\mathrm{s}}^{2}} \\ {\omega _{n}}^{2}D\ddot{y}_{2}+\left( 2\xi {\omega _{n}}^{2}D+\dfrac{\gamma \omega _\mathrm{s}\omega _{n}\rho D^{3}}{M}\right) {\dot{y}}_{2}\\ +{\omega _{n}}^{2}Dy_{2}=\dfrac{\frac{1}{2}\rho U^{2}DC_{L}}{M}+\dfrac{\frac{1}{2}\rho U^{2}D}{M}B_{Y}\dfrac{{\omega _{n}}^{2}D\ddot{y}_{2}}{D{\omega _\mathrm{s}}^{2}}.\\ \nonumber \end{array}\right. }\!\!\!\!\!\!\!\\ \end{aligned}$$

(16)

It should be emphasised that modification terms were added to structure equations to improve how the model is describing the wake force. Thus, the curve fitting was applied to the difference between downstream and upstream force coefficients. They are describing the change in inertia force due to the chaotic upstream wake. Furthermore, these terms are hydrodynamic coefficients, so they should be multiplied by free stream dynamic pressure. In this way, the increase in turbulence with a rise of Re in the wake will be considered. By applying non-dimensional time and distance to Eq. (16), the structural motion equation becomes:

$$\begin{aligned} {\left\{ \begin{array}{ll} \ddot{x}_{2}+\left( 2\xi +\dfrac{2\gamma \omega _{0}}{\mu }\right) {\dot{x}}_{2}+x_{2}=\dfrac{1}{8}\dfrac{C_{D_{2}}}{\mu \pi ^{2}\mathrm{St}^{2}}+\dfrac{1}{2}\dfrac{A_{X}}{\mu \mathrm{St}^{2}}\ddot{x}_{2} \\ \ddot{y}_{2}+\left( 2\xi +\dfrac{\gamma \omega _{0}}{\mu }\right) {\dot{y}}_{2}+y_{2}=\dfrac{1}{8}\dfrac{C_{L_{2}}}{\mu \pi ^{2}\mathrm{St}^{2}}+\dfrac{1}{2}\dfrac{B_{Y}}{\mu \mathrm{St}^{2}}\ddot{y}_{2}.\\ \nonumber \end{array}\right. }\!\!\!\!\!\!\!\!\\ \end{aligned}$$

(17)

Rearranging this equation yields:

$$\begin{aligned} {\left\{ \begin{array}{ll} \left( 1-\dfrac{1}{2}\dfrac{A_{X}}{\mu \mathrm{St}^{2}}\right) \ddot{x}_{2}+\left( 2\xi +\dfrac{2\gamma \omega _{0}}{\mu }\right) {\dot{x}}_{2}+x_{2}\\ \quad =\dfrac{1}{8}\dfrac{C_{D_{2}}}{\mu \pi ^{2}\mathrm{St}^{2}} \\ \left( 1-\dfrac{1}{2}\dfrac{B_{Y}}{\mu \mathrm{St}^{2}}\right) \ddot{y}_{2}+\left( 2\xi +\dfrac{\gamma \omega _{0}}{\mu }\right) {\dot{y}}_{2}+y_{2}\\ \quad =\dfrac{1}{8}\dfrac{C_{L_{2}}}{\mu \pi ^{2}\mathrm{St}^{2}}.\\ \nonumber \end{array}\right. }\\ \end{aligned}$$

(18)

The new acceleration term can be interpreted as a change in the added mass coefficient due to the increase in turbulence in the gap between two cylinders. In what follows, this term is referred to as added mass modification coefficient. The chaotic flow regime in the upstream wake changes the added mass coefficient of the cylinder in comparison with unit value suggested by Blevins [21]. This term is a function of St which is governed by leading cylinder. Table 1 includes stream-wise and cross-flow added mass modification coefficient obtained by curve fitting to experiment results [5], see Fig. 6.

Table 1 Added mass modification coefficients obtained from curve fitting for different spacings

Analytical solution of such a complex system is a challenge; thus, it was solved numerically. MATLAB program provides a suitable tool through SimuLink Add-on.

Two Simulink models were initially developed so that Eqs. (13) and (18) could be solved simultaneously. Fourth-order Runge–Kutta algorithm with variable time step (to enhance convergence and stability) was adopted. Reduced velocity was increased gradually with increments of 0.2 from zero. Initial conditions for all reduced velocities were considered similar at $t = 0;\; C_L = C_D = 2$ and $x = y = 0$. It should be noted that several initial conditions were tested and it was concluded that this model is not sensitive to the initial conditions. SimuLink simulations were run for 400s for each reduced velocity so that a steady-state solution was obtained.

Figure 7 presents the simulation results from the modified model for a range of reduced velocities against experimental results at the corresponding spacing. Added mass modification coefficient widens the lock-in range at all spacings. It is established in the literature that variation of mass ratio has a significant effect on the width of the lock-in range response. Thus, the effect of the modification term appears to be similar to that of the mass ratio.

Fig. 7

Experimental response amplitude versus mathematical model simulation at various reduced velocity and spacings. a, b $L/D=3.5$, c, d $L/D=4$, e, f $L/D=5$, g, h $L/D=8$, i, j $L/D=10$, k, l $L/D=15$, m, n $L/D=20$

Nevertheless, variation in spacing has an insignificant effect on the oscillation amplitude which is in contrast with observations from experimental investigations. Additionally, the oscillation amplitude predicted by the model is much higher than the experimental results. The difference between the model and experiment results is not constant at various reduced velocities; at lower velocities, the model and experiment have a better agreement, whereas at intermediately high velocities, the model yields larger oscillation amplitudes. However, the amplitude predicted by the model drops below experimental results at very high velocities. Such a self-limiting characteristic is similar to VIV phenomenon itself. It suggests that a secondary modification term that can damp the amplitude at medium velocities and increase it at higher velocities can resolve this issue. Thus, a the second modification term can be defined as a function of the cylinder velocity alike the added mass modification term.

4.2.2 Added damping modification term

The damping term may be considered through a similar approach previously used. Therefore, it is introduced to the motion equation as a force coefficient. Additionally, since this term is reflecting the effect of the upstream vortices, it should be a function of their shedding frequency, and two new modification coefficients are defined as $E_{X}\dfrac{{\dot{X}}_{2}}{D{\omega _\mathrm{s}}}$ and $F_{Y}\dfrac{{\dot{Y}}_{2}}{D{\omega _\mathrm{s}}}$. ($E_{X}$ and $F_{Y}$ coefficients are functions of spacing to be determined using experimental results.) Moreover, the difference between upstream and downstream mean drag must be considered as well. Therefore, the term $C_1$ is added to the equation due to shielding effect. So the dimensionless equations of motion become:

$$\begin{aligned} {\left\{ \begin{array}{ll} \left( 1-\dfrac{1}{2}\dfrac{A_{X}}{\mu St^{2}}\right) \ddot{x}_{2}+\left( 2\xi +\dfrac{\omega _{0}}{\mu }\left( 2\gamma -\dfrac{E_{X}}{2St^{2}}\right) \right) {\dot{x}}_{2}+x_{2}\\ =4a\left( C_{D_{2}}+C_{1}\right) +2\pi aC_{L_{2}}\dfrac{\dot{y}_{2}}{U_\mathrm{r}} \\ \left( 1-\dfrac{1}{2}\dfrac{B_{Y}}{\mu St^{2}}\right) \ddot{y}_{2}+\left( 2\xi +\dfrac{\omega _{0}}{\mu }\left( \gamma -\dfrac{F_{Y}}{2St^{2}}\right) \right) {\dot{y}}_{2}+y_{2}\\ =aC_{L_{2}}-8\pi aC_{D_{2}}\dfrac{\dot{y}_{2}}{U_\mathrm{r}}. \end{array}\right. }\nonumber \\ \end{aligned}$$

(19)

This equation includes geometrical nonlinearity term as well $\left( 2\pi aC_{L_{2}}\dfrac{\dot{y}_{2}}{U_\mathrm{r}}, 8\pi aC_{D_{2}}\dfrac{\dot{y}_{2}}{U_\mathrm{r}}\right)$. These new constants $\left( E_{X},F_{Y},C_{1}\right)$ should be determined through an optimisation process in such a way that the difference between experimental results and model simulation becomes minimum.

4.2.3 Optimisation

The objective of the optimisation is to determine $E_{X},F_{Y}$ and $C_{1}$ so that the difference between mathematical model simulation and experimental results becomes minimum. The optimisation function [Eq. (20)] was determined in such a way that the accumulative error between the experiment and the model in cross-flow direction is minimised. The error function was limited to the transverse direction, firstly for the sake of simplicity and secondly, it was observed by Srinil et al. [31] that stream-wise simulation is heavily influenced by cross-flow results, and hence, every change in cross-flow model could significantly alter the simulation in either direction, while changes in stream-wise direction have negligible effect on the predicted amplitude. Therefore, optimisation in cross-flow motion is considered sufficient.

$$\begin{aligned} e=\dfrac{\sum |A_{Y_{2}}-y_{2}|}{n} \end{aligned}$$

(20)

fminsearch command in MATLAB was used for optimisation. This command can calculate the local minimum of a discontinuous function with multi-variables using the derivative-free method. The algorithm used for fminsearch is Nelder–Mead simplex algorithm [36]. No constraint was set for the error function or any of the variables. Also, fminsearch command requires no constraint for optimisation function. Moreover, the initial guess for $\left( E_{X},F_{Y},C_{1}\right)$ was (0, 0, 0) for initial spacing of $L/D = 3.5$, and then, optimisation results from smaller spacings were used as the initial guess for consecutive spacings. Options chosen for optimisation can be seen in Table 2. It should be noted that fminsearch terminates optimisation process when conditional tolerances on variables and function value are satisfied simultaneously, see Table 2.

Table 2 Options and their designated values for optimisation

The result of optimisation can be seen in Table 3 with the corresponding spacing. Figure 8 demonstrates simulation results from Eq. (19) using Table 3 values.

Table 3 Optimisation output for three modification parameters of $E_{X}$, $F_{Y}$ and $C_{1}$

Fig. 8

Experimental response amplitude versus fully modified mathematical model simulation at various reduced velocity and spacings. a, b $L/D=3.5$, c, d $L/D=4$, e, f $L/D=5$, g, h $L/D=8$, i, j $L/D=10$, k, l $L/D=15$, m, n $L/D=20$

Fig. 9

Experimental VIV response amplitude versus fully modified mathematical model simulation at various reduced velocity and spacings. a, b $L/D=3.5$, c, d $L/D=4$, e, f $L/D=5$, g, h $L/D=8$, i, j $L/D=10$, k, l $L/D=15$, m, n $L/D=20$

It can be observed that the mathematical model has a good agreement with experimental data in cross-flow direction; however, their agreement with the stream-wise response is relatively poor. Moreover, if the capabilities of the model in simulation of downstream fluid-induced vibration (FIV) response are to be evaluated by its ability to accurately simulate four parameters of lock-in range width, lock-in onsets velocity, velocities at which maximum amplitude occurs and its magnitude, following observations can be made, see Fig. 8:

• Velocity onset of lock-in range is predicted correctly for both cylinders.

• Lock-in range width is simulated successfully for both cylinders. However, it is less accurate for the trailing cylinder in very large spacings.

• Reduced velocity at which maximum amplitude occurs is predicted accurately for both cylinders at all spacings.

• Oscillation amplitude at a given reduced velocity is simulated more accurately compared to the existing mathematical models for an isolated cylinder. Although it is evident that the model under-predicts the results, nevertheless, the error is less than that of existing models in the literature for a single cylinder.

• The first peak in the cross-flow response of the trailing cylinder is captured by the model as well (Fig. 8b, c, f). However, it disappears as the spacing grows large.

• The reduction of amplitude due to increase in spacing can be captured by the model successfully.

Figure 9 includes trailing cylinder response which consists of VIV motion and displacement induced by upstream wake interference. Moreover, Armin et al. [5] observed that structural motion due to VIV and WIV excitation mechanisms can be separated based on their frequency. They concluded that the motion related to irregular collisions of the upstream vortices has higher frequency. Therefore, if the irregular motion with the higher frequency was to be eliminated (for more information on separation procedure of high-frequency from low-frequency motions please, see Armin et al. [5]), the model simulation could be appreciated more. Figure 9 depicts a comparison between the model results and experimental result with the high-frequency amplitude removed. It is clear that the under-prediction problem in simulation results is no longer an issue. Moreover, simulation in the stream-wise direction appears more successful. Although it fails to capture the maximum amplitude value in this direction, the overall behaviour of the cylinder is captured in all spacings.

As mentioned in the previous section, the cross-flow model has a significant effect on the stream-wise response. Hence, the large pick in stream-wise simulation at approximately $U_\mathrm{r} = 6$ at small spacings occurs due to existence of a peak in cross-flow response at the corresponding reduced velocity even though it has no corresponding peak in the experiment results. It is clear that this peak disappears at large spacings where the cross-flow response peak fades.

Universal functions that can determine the modification coefficients $(A_X, B_Y, E_X, F_Y$ and $C_1)$ at different spacings, without depending on experimental results, are necessary. An attempt to develop such functions with the space between the two cylinders as their variables is presented in “Appendix.” These functions are obtained through a series of curve fitting to the values given in Tables 1 and 3 .

5 Conclusion

The well-received concept of a coupled system of a wake oscillator and a structural motion equation was employed to develop a time-domain model for simulation of interference between two cylinders in tandem. It was assumed that cylinders were rigid and flexibly mounted with identical structural stiffness and mass ratios in cross-flow and stream-wise directions. Hydrodynamic coefficients for both cylinders were considered to be equal for the sake of simplicity and compensating for the lack of experimental measurements for the trailing cylinder hydrodynamic coefficients.

The model was developed based on a coupled system of a van der Pol wake oscillator and a Duffing equation. Duffing equation was considered to describe the structural response of the cylinders in both directions, and it is capable of capturing the structural nonlinearity of the system. Van der Pol wake oscillator is well received in the literature as it can capture the self-exciting and self-limiting nature of VIV. The wake oscillator was coupled to the structural motion through a linear function of structural acceleration. The excitation term in Duffing equation was the wake force which was obtained from a van der Pol equation.

Furthermore, the excitation term of trailing cylinder was modified to consider the buffeting impact of the vortices in upstream wake. Any input from the upstream cylinder was avoided, and two modification terms were added to adjust the added mass coefficient and added fluid damping as functions of Strouhal number so that the effect of upstream wake instability could be incorporated in the model.

The hydrodynamic force exerted on each cylinder was calculated from experimental data through a simple motion equation for a mass–spring and damper system. Since the model simulated upstream cylinder behaviour with a good agreement, the difference between upstream and downstream hydrodynamic forces was calculated to determine the additional force due to wake interference. Functions of upstream displacement, downstream acceleration and velocity were fitted to these values. It was concluded that trailing cylinder acceleration governs this extra force. Two added mass modification coefficients were defined as $A_{X}\dfrac{{\ddot{X}}_{2}}{D{\omega _\mathrm{s}}^{2}}$ and $B_{Y}\dfrac{{\ddot{Y}}_{2}}{D{\omega _\mathrm{s}}^{2}}$, and a linear equation was fitted to upstream and downstream wake force variance. These modification coefficients were considered as forces in Duffing oscillator equations and were multiplied by dynamic pressure of the free stream. $A_X$ and $B_Y$ were determined by curve fitting to the variance.

It was observed from the simulation results obtained from the new model that predicted maximum amplitude is significantly higher than that from the experiment. It was concluded that this discrepancy was due to the increase in damping caused by the turbulent flow in the wake. The damping was adjusted by a secondary force modification coefficient as $E_{X}\dfrac{{\dot{X}}_{2}}{D{\omega _\mathrm{s}}}$ and $F_{Y}\dfrac{{\dot{Y}}_{2}}{D{\omega _\mathrm{s}}}$ which was multiplied by dynamic pressure of the free stream as well and added to the right-hand side of the Duffing equation. $E_X$ and $F_Y$ constants were derived through optimisation through which the error between experimental and simulation results was minimised. It should be mentioned that optimisation was carried out by fminsearch command in MATLAB with no constraint on the function value or any of the variables.

Based on the performance of the modification coefficients, it was concluded that wake interference reduces the added mass coefficient and increases viscous damping (added damping coefficient).

The final model was capable of:

• Predicting upstream and downstream lock-in onset velocity

• Simulating the lock-in range width for both cylinders

• Predicting the reduced velocity at which maximum amplitude occurs for both cylinders at all spacings

• Predicting the oscillation amplitude at a given reduced velocity more accurately compared to the existing mathematical models for an isolated cylinder.

Energy transfer between fluid and structure is another important aspect of such investigations that will be addressed in the future by investigating added mass and other hydrodynamic coefficients.

Appendix

At this stage, the simulation has been compared against experiment at seven spacings and it is possible to introduce universal coefficients for added mass and damping which can evaluate them at different spacings in the form of Eq. (21).

$$\begin{aligned} {\left\{ \begin{array}{ll} A_{X}=f(L/D) \\ B_{Y}=f(L/D) \\ E_{X}=f(L/D) \\ F_{Y}=f(L/D) \\ C_{1}=f(L/D) \end{array}\right. } \end{aligned}$$

(21)

Such a task can be accomplished by curve fitting to seven values of each parameter obtained in the previous section (Tables 2, 3). Figure 10 demonstrates the results of Gaussian function fitted to those four values and a linear function to $C_{1}$. The results of the curve fitting process can be found in Table 4.

Table 4 Modification coefficients as a function of spacing

Fig. 10

Modification coefficients curve fitting results to Gaussian function. a $A_{X}$, b $B_{Y}$, c $E_{X}$, d $F_{Y}$ and e linear function to $C_{1}$

Figure 10 shows the variation of each coefficient as the gap between two cylinders changes. Based on the experimental data, it is not clear after what distance the influence of upstream wake becomes insignificant; however, the fact that all functions describing different coefficients reach zero mathematically at large spacings (approximately 32D) is endorsing.