An analysis method of the vortex-induced vibrations of a tethered sphere

Abstract

Vortex-induced vibrations (VIV) in systems with more than one degree of freedom often present complex synchronization among the motion components, also hidden by the randomness that characterizes the motion itself. A phase average method has been here developed and applied to the displacements of a tethered sphere, at low mass and damping, to analyze its xy trajectories over a wide range of reduced velocities, 5 ≤ U* ≤ 25 (Reynolds numbers, 5.1 × 10³ ≤ Re ≤ 2.67 × 10⁴). This method has allowed the identification of both the periodic and chaotic contribution of each motion component, accurately reconstructing the underlying trajectory periodic pattern. The two classical vibration modes, I and II, have been also observed. The method developed here was able to better rebuild the experimental data compared to other methods found in the relevant literature, providing useful insights into the study of the dynamic response of a freely-oscillating tethered sphere immersed in a steady flow.

1 Introduction

The dynamic response analysis of bodies completely or partly immersed in steady flows is crucial in marine engineering, when analyzing the stability of various off-shore structures including production and drilling risers, pipelines, moorings, tethers of tension leg platforms, spar platforms and the members of jacket structures. Interest in this field of research has increased, given its relevance in green-energy exploitation techniques (e.g., [1]).

Most of the existing works focus on bodies with multiple degrees of freedom, and mainly cylinders. Only a few of them investigate the flow interactions with a sphere, due to the difficulties encountered in measuring, modelling, and analyzing its complex three-dimensional wake [2,3,4].

Early studies evaluated the action of surface waves on tethered buoyant spheres, empirically obtaining drag and inertia coefficients through Morison’s equation [5, 6]. Later, the focus moved to analyzing the tethered sphere oscillations, taking into account the mass ratio, m* (ratio between the mass of the sphere, m, and the mass of displaced fluid, m_d), which defines a body as ‘light’ (m* < 1) or ‘heavy’ (m* > 1).

In 1997, Williamson and Govardhan [7] and Govardhan and Williamson [8] mainly investigated the VIV of light tethered spheres, finding large crossflow oscillation amplitudes with typical figure-of-eight or crescent topologies trajectories. When experimenting with heavier spheres, the ratio between crossflow and streamwise oscillations was observed to be larger than that of light spheres. In both types of sphere, they detected the classical resonance condition and a synchronisation regime, which was reduced for heavy spheres, with a saturation amplitude close to one diameter.

In 2001, Jauvtis et al. [9] named the resonance condition and the saturation amplitude as Mode I and II, respectively, and found two new modes of dynamic response through wind tunnel experiments on heavy tethered spheres (80 ≤ m* ≤ 940) at reduced velocity, U* (ratio between the mean streamwise velocity, U, and the product of the sphere diameter, D, with natural frequency of the sphere in water, f_n), between 0 and 300: one periodic Mode III, at 20 ≤ U* ≤ 40 for m* = 80, and a Mode IV, at U* > 100, characterized by intermittent bursts of large-amplitude vibrations.

In 2005, Govardhan and Williamson [10] extended their previous experimental research on tethered spheres, which were shown to exhibit non-periodic vibrations when passing from Mode I to Mode II. In the case of light spheres, such transition regime was clearly distinct and characterized by a difference in the phase of the vortex force, consistent with a large difference in the timing of the vortex formation between the two modes.

Firstly in 2003, and later in 2005, Provansal et al. [11, 12] discussed how a different configuration of a tethered sphere towards the flow direction modified the shape of the trajectories. In particular, they analyzed the VIV of a heavy sphere (m* = 2.433), at Re = UD/ν (where ν is the cinematic viscosity of the fluid) between 600 and 800, in a flow parallel to the thread, observing different modes of body oscillation along circular, straight, and elliptical trajectories with increasing reduced velocity.

Studies on a heavy tethered sphere (m* = 7.87) were also carried out by van Hout et al. [2] who, at 2.8 ≤ U* ≤ 31.1, discovered three different bifurcation regions: in the first region, the sphere did not exhibit any crossflow oscillations but only a slight streamwise deflection; the second region overlapped with Mode I and II described by Jauvtis et al. [9] and Govardhan and Williamson [10]; in the third region, the sphere exhibited large crossflow non-periodic vibrations, exceeding those in the second region and comparable to the Mode III of Jauvtis et al. [9]. Similar results were obtained by the same authors [13] for even heavier tethered spheres (m* = 22.20) and for a narrower range of reduced velocities (3.18 ≤ U* ≤ 14.1). Eshbal et al. [14] experimentally analyzed the three-dimensional wave structures of a heavy tethered sphere for low Reynolds numbers (382 < Re < 1392), identifying different vortex configurations: at U* = 3.6, the sphere did not oscillate and hairpin vortices with vertical plane of symmetry were periodically shed (similarly to the stationary sphere); at U* = 5.8, the sphere motion was highly periodic and vortices with horizontal plane of symmetry were alternatively shed from the two sides of the sphere; for larger U*, the sphere oscillations were intermittent and the vortex shedding less organized.

Parallel experiments on light tethered spheres (m* = 0.392) in uniform flow at reduced velocities, 2.2 ≤ U* ≤ 10.1, corresponding to 430 ≤ Re ≤ 1925, were carried out by Eshbal et al. [15]. They detected two different bifurcation regions: the first one in which the sphere did not oscillate, coinciding with the onset of an asymmetric shedding of hairpin vortices observed for a stationary sphere, and the second (lock-in) region, which started at U* ≈ 2.2. In 2013, the same authors [3] analyzed both light (m* = 0.392) and heavy (m* = 7.87) tethered spheres for 2.2 ≤ U* ≤ 13.5 which, after an initially stationary state, started to oscillate at U* > 2.2 (Re = 430) for m* = 0.392, as showed in the previous study [15], and at U* > 3.5 (Re = 673) for m* = 7.87.

In the last decade, more attention has been paid to the study of the dynamic response of elastically mounted spheres in steady fluid flows. Among the first researchers to introduce elastically mounted spheres were Govardhan and Williamson [10]. They experimentally analyzed bodies of different mass ratio forced to move only in the flow crossflow direction, noting three responses (Mode I, II, and III), similar to those found for tethered spheres. For lower m*, the transition between the first two modes was quite distinct, exhibiting a jump, whereas for higher m* a more continuous transition regime was detected.

Computationally, Behara et al. [16] studied the vibrations of an elastically mounted sphere of m* = 2 with three degrees of freedom at fixed Re = 300 and reduced velocity in the range of U* = 4-9. They observed two distinct response branches corresponding to two different vortex shedding (hairpin and spiral) modes and to two different sphere oscillation trajectories (a nearly straight line on the crossflow plane and a perfectly circular orbit). Later, the same authors [4] expanded the range of reduced velocities (U* = 0-13) at Re = 300 and the range of Reynolds numbers to 300 ≤ Re ≤ 1000, fixing U* = 9. They showed how the vortex shedding mode and the sphere motion were both strongly influenced by the Reynolds number. In particular, for Re ≤ 400, in the spiral mode of the vortex shedding, the sphere moved along circular orbits on the crossflow plane; for Re = 500-600, however, passing from the spiral to the hairpin mode, the trajectory changed from circular to elliptical orbits.

Further interesting numerical results were obtained by Rajamuni et al. [17] on an elastically mounted sphere (m* = 2.865) forced to move only transversally, through two sets of simulations, i.e. at Re = 300 with 3.5 ≤ U* ≤ 100 and at Re = 800 with 3 ≤ U* ≤ 50. For Re = 300, large-amplitude (peaking at 0.4D) and highly periodic sphere oscillations were observed, showing a resonance response of amplitude 0.6D in the U* = 4.5–13 range with increasing Reynolds numbers up to 800. At both Re = 300 and 800, two sets of hairpin-type vortex loops were detected around the sphere, in agreement with Govardhan and Williamson [10] and Behara et al. [16].

Differently from the above-mentioned studies, Negri et al. [18] experimentally investigated an elastically mounted sphere (m* = 1.24) under a range of reduced velocities between 1.6 and 13.6. The sphere showed a behaviour similar in amplitude and frequency to the tethered sphere, exhibiting Mode I (5 ≤ U* ≤ 8), with oscillations between 0.4D and 0.65D, and Mode II (10 ≤ U* ≤ 14), characterized by higher crossflow amplitudes between 0.8D and 1D. In addition, the authors detected other trajectories different from the classical 8-shaped one, such as the shape of a segment, a “rain-drop”, and a D-shape.

Lately, the research focus moved more and more to analyzing the influence of the free surface on the movement of elastically mounted spheres. Experimental studies of Mirauda et al. [19,20,21,22,23,24] investigated a heavy elastically mounted sphere (m* = 1.34) at high Reynolds numbers (Re = 1.7·10⁴–7.0·10⁴) and for a range of U* between 1.90 and 7.58. When the sphere was fully submerged, it showed periodic crossflow oscillations and a transition trajectory between crescent topology and figure-of-eight shape. Moving towards the free surface, the sphere displayed non-periodic crossflow displacements and a trajectory of chaotic shape.

Sareen et al. [25] also experimentally investigated the effects of free surface on fully and partly submerged heavy elastically mounted spheres (m* = 7.8), free to oscillate only transversally. The fully submerged sphere displayed highly periodic vibration amplitudes, progressively increasing from Mode I to Mode II, and a Mode III characterized by small intermittent and clearly non-periodic vibrations. The partly submerged sphere, however, showed two different responses: Regime I, where the presence of the free surface widened the synchronization region and sometimes increased the vibration amplitudes, and Regime II in which, with the sphere moving nearer the free surface, a very sharp reduction of amplitude, and the two distinct modes, I and II, were observed.

Most of the above-mentioned studies underline how VIV trajectories are characterized by a large amount of periodicity but they can also contain a significant chaotic part (e.g. [26]). Moreover, when dealing with two-degree-of-freedom (xy) VIVs, the synchronization between the two components of motion may not be easy to detect. The identification of the periodic and chaotic contribution of each motion component, and hence of the trajectory pattern, provides useful insights into the fluid–structure interaction. “Classical” xy VIVs are 8-shaped, and they are encountered in elastically mounted circular cylinders as well as spheres and in tethered spheres (e.g. [9, 27]). Their VIVs are characterized by the presence of one frequency, f, in the crossflow component and one frequency, 2f, in the streamwise component of motion. Other kinds of periodic trajectories have been recently catalogued in circular cylinders with different natural frequencies in the streamwise (f_nx) and crossflow (f_ny) directions, particularly when f_nx/f_ny< 1 [28], as well as in square cylinders subjected to shear flow [29] and in elastically mounted spheres [18]. Resembling trajectories already observed in flexible cylinders [30], named “drop” or “egg” type in Kang and Jia [28], occur when the streamwise motion component includes the first harmonic (which is the frequency of the crossflow motion), besides the second harmonic. The topological description of these trajectories was presented in Negri et al. [18]. Trajectories with some degree of chaos were also shown in Dahl et al. [26] for a cylinder with different natural frequency in x and y, although the underlying periodic pattern could not be recognized. Kang and Jia [28] assessed that, for a circular cylinder, these non-classic trajectories, and hence the appearance of the first harmonic in the streamwise motion, are due to the different natural frequencies of the streamwise and crossflow motion. On the other hand, Sun et al. [29] showed that, for a square cylinder, the non-classic trajectories are due to the non-uniformity of the velocity profile. Finally, Negri et al. [18] observed these non-classical trajectories on a sphere fixed to a cantilever beam, showing that the shape of the trajectory mostly depends on the phase lags between the crossflow component and the two harmonics of the streamwise component.

Different methods exist to characterize the periodic pattern of a xy VIV trajectory and, in particular, to calculate the phase lag between the crossflow and streamwise component of motion. To the authors’ knowledge, they are:

• (1) Fourier analysis of the motion components for calculating average amplitudes and phases of each component.

• (2) Regression of the motion components on a sum of few sinusoids with frequencies multiple of the main one [28]. This method provides the average amplitudes and phases of each component.

• (3) Calculation of the instantaneous phases of x and y through the analytical signals [30] or through the linear interpolating phase (e.g. [31]).

• (4) Phase average of x and y on the basis of the periodic component (either x or y) and calculation of the first two Fourier coefficients of the phase-averaged quantities [18]. Even in this case, the method provides the average amplitudes and phases of each component.

To the best of the authors’ knowledge, the fourth method has been applied to VIV trajectories in only a few works, involving cylindrical and spherical structures [18, 32, 33], despite its suitability for this matter, since VIV trajectories often exhibit a sinusoidal component of motion, which can be used as phase reference. A possible reason for this infrequent use of the phase average is that classical and very periodic trajectories (8-shaped or crescent) are usually encountered in VIV studies, and methods 1-3 are also effective in recognizing the periodic pattern of these trajectories. Nonetheless, phase average is a common technique in analyzing periodic signals, and it has been used in many VIV works to represent the velocity and vorticity fields of the wake. For example, some studies such as Krakovich et al. [3], Govardhan and Williamson [10] and Eshbal et al. [15] applied this technique to experimentally investigate the wake of tethered spheres undergoing VIVs.

Phase-averaging methods are well summarized in the work of Ostermann et al. [34], where the experimental analysis of a fluidic oscillator was performed. They mainly subdivide into signal conditioning and mathematical methods. In detail, in the signal conditioning methods a time-resolved signal is used as phase reference, so that it identifies the oscillation periods. This reference can be an external signal or the signal itself to phase-average. There are two main ways to identify periods inside a signal: zero-crossing and autocorrelation. Mathematical methods include the Hilbert transform and the proper orthogonal decomposition (POD). In VIV studies, a signal conditioning method is usually adopted, since the periodic motion of the oscillating body provides a phase reference for the flow field and for the motion of the body itself.

In the work of Cagney and Balabani [32], the oscillation and the wake of a pivoted cylinder with one and two degrees of freedom (DOF) for low values of reduced velocity (U* < 3.5) were experimentally analyzed. The raw displacements of the cylinder were phase-averaged on the basis of the crossflow motion, and an 8-shaped pattern was recognized. Then, the phase-averaged streamwise and crossflow components of motion were fit with sinusoidal waves with frequency one (streamwise) double of the other (crossflow).

In 2015, Thomas et al. [33] experimentally analyzed the motion and the wake of a flexible cylinder. By applying the phase-average on the raw displacements of the cylinder, the authors managed to extract the periodic pattern of almost the sections of the cylinder, which was a slightly asymmetrical 8-shaped trajectory. Only for a particular section of the cylinder, the phase average applied by the authors could not reveal a precise periodical pattern.

In the work of Negri et al. [18], the phase average applied to the VIV trajectories of an elastically mounted sphere allowed the recognition of periodic patterns different from the classical 8-shaped or crescent type and validating a trajectory model based on the presence of (at most) two harmonics in the streamwise component of motion. With this model, even the less common trajectories encountered in the VIV literature were modelled: asymmetrical 8-shaped, D-shaped, and drop-shaped.

In this work, the phase-average method is extensively described and used for analyzing the VIVs of a tethered sphere. The method here presented is based on the self-synchronized phase average presented in Negri et al. [35] for the analysis of particle-streak-velocimetry (PSV) velocity fields and used in Negri et al. [18] for the analysis of the VIV trajectories of a cantilevered sphere. Here, the technical solutions developed for phase-averaging highly chaotic oscillations, and particularly for properly representing the average synchronization between the periodic component of motion (used as phase reference) and the secondary component of motion, are shown in detail. In addition, the present method proved its validity and accuracy through comparison with the other ones existing in literature (methods 1–3 listed above).

The paper is organised as follows: the experimental details are presented in Sect. 2; the method of data analysis is described in Sect. 3; the application of the analytical method on an experimental case and the dynamic response of a tethered sphere under VIVs are discussed in Sect. 4; the main findings of the present work are summarized in Sect. 5.

Main nomenclature
Axn	Amplitude of the nth harmonic of the streamwise motion
Ay	Crossflow oscillation amplitude
c	Sphere damping in still water
CA	Potential added mass coefficient (CA = 0.5 for the sphere)
D	Sphere diameter
f	Main oscillation frequency of the sphere
fx	Oscillation frequency of the sphere in x direction
fy	Oscillation frequency of the sphere in y direction
fn	Natural frequency of the sphere in water
fs	Sampling frequency
fvo	Vortex shedding frequency with stationary sphere
Fy0	Amplitude of the total crossflow force
Fy,v0	Amplitude of the vortex crossflow force
h	Water level
hs	Distance between the top of the sphere and the free surface (still liquid and sphere)
Jk	Number of samples in the kth phase bin
K	Number of phase bins
L	Tether length (distance between pivot and center of mass of the sphere)
m	Sphere mass
mA	Sphere added mass
md	Displaced fluid mass
M	Number of samples from each period in a phase bin
n’	Index of goodness of the modelled trajectory
N	Number of zero-crossing points
r	Reference signal
t	Time
t0	Instant of zero-crossing in the reference signal r
T	Reference period of oscillation (approximation of T0)
Tmin	Minimum period of oscillation
T0	Average period of oscillation
U	Mean streamwise velocity
x	Streamwise position of the sphere center
y	Crossflow position of the sphere center
$\theta_{K}$	Vector of bin centers for the phase average
ν	Kinematic viscosity of the fluid
$\xi$	Generic component of motion
$\xi_{ph}$	Phase average of $\xi$
$\sigma_{\xi }$	Standard deviation of $\xi$
$\sigma_{\xi ph}$	Phase standard deviation of $\xi$
$\sigma_{T}$	Standard deviation of the oscillation period
$\sigma_{{y^{{\prime }} }}$	Standard deviation of the random part of y
$\phi$	Instantaneous phase of a signal
$\phi_{x}$	Instantaneous phase of the streamwise position
$\phi_{y}$	Instantaneous phase of the crossflow position
$\phi_{xy1}$	Phase lag between streamwise and crossflow position (first harmonic)
$\phi_{xy2}$	Phase lag between streamwise and crossflow position (second harmonic)
$\varPhi_{y}$	Phase lag between crossflow force and displacement
$\varPhi_{y,v}$	Phase lag between crossflow vortex force and displacement

2 Experimental apparatus and measurement techniques

The experiments were performed in a re-circulating flume at the Hydraulics Laboratory of the Politecnico di Milano. The test section had a total length of 6 m, a width of 0.5 m and a height of 0.6 m and it was made of Plexiglas to ensure optical access from all sides. The mean streamwise velocity, U, was varied in a range of 0.09–0.40 m/s, corresponding to the Reynolds numbers between 5 × 10³ and 2.5 × 10⁴, while the water level was maintained at h/D = 7 for all experiments. The distance of the sphere from the free surface was h_s/D = 3 (considering the sphere at rest).

A Plexiglas sphere with diameter D = 0.06 m and m* equal to 1.38, covered in black paint to avoid reflection effects, was mounted far from the inflow section. A dyneema wire (braided fishing line) of diameter $D_{t}$ = 0.15 mm, two orders of magnitude smaller than the sphere diameter, acted as tether and was positioned from a fixed aluminium structure, mounted above the channel, down to the top of the body (Fig. 1). This way, the body was allowed to move on a spherical surface. The non-dimensionalized tether length, L* (ratio between the tether length, L, and the sphere diameter, D), was varied between 6.1 and 11.7 in order to increase the reduced velocity without increasing the water velocity too much (and causing free surface undulations), while still assuring the small-angle approximation.

Fig. 1

a Experimental scheme. b Snapshot of the channel with the sphere, the wire and its support

A charge-coupled device (CCD) camera, located below the flume, acquired the sphere oscillations in the xy plane at a frequency of 10 ÷ 20 Hz and with a resolution of 659 × 493 px. The technique used was to record the motion of a fluorescent circular marker placed at the bottom of the sphere (aligned with its center of gravity) in the darkness, the same way a bright point is detected on a black background. Subsequently, the centroid of the blob associated to the marker was evaluated after converting each frame into a binary image at a set brightness threshold [18].

Since the sphere, during its motion, significantly rose from its position at rest, the distance between the camera and the horizontal plane that included the sphere marker increased, so the image resolution [px/m] decreased. Therefore, a perspective correction was applied to the data, according to the procedure reported in the “Appendix”. This procedure basically considered (1) the camera image resolution [px/m] as a function of the distance between the camera and the sphere marker and (2) the constraint of the sphere marker moving on a spherical surface.

Before the flowing water tests, free decay oscillations in still water were performed, following the procedure described in Negri et al. [18], to obtain the natural frequency, f_n, and the damping ratio of the sphere in still water, $\zeta = \frac{c}{{4\pi f_{n} \left( {m + m_{A} } \right)}}$ (where c is the sphere damping in still water and m_A is the sphere added mass), assuming logarithmic decay of the oscillation amplitude.

3 Method of data analysis

The phase average method here presented is a signal conditioning method and deals with time-resolved signals. Nonetheless, many equations are expressed in the continuous-time formulation for the sake of simplicity, even though all the quantities involved in the analysis of the data are actually discrete-time functions.

Starting from the VIV trajectory analysis carried out preliminarily by Negri et al. [18] on an elastically-mounted sphere, the method is here furtherly developed and explained in details. Firstly, some general considerations about the phase average applied to VIV are deduced; secondly, the specific procedure of the developed phase-average method is presented.

Each component of motion $\xi$ (x or y) of a body under VIVs can be decomposed as the following:

$$ \xi \left( t \right) \equiv \bar{\xi } + \bar{\xi }\left( t \right) + \xi^{{\prime }} \left( t \right) $$

(1)

where $\bar{\xi }$ is the long-time average, $\bar{\xi }\left( t \right)$ is the periodic part, and $\xi '\left( t \right)$ is the random part. $\tilde{\xi }\left( t \right)$ and $\xi^{{\prime }} \left( t \right)$ have zero mean value.

The amplitude of $\xi \left( t \right)$ in terms of standard deviation, $\sigma_{\xi }$, is:

$$ \sigma_{\xi } = \sqrt {\frac{1}{t}\smallint \left( {\xi \left( t \right) - \bar{\xi }} \right)^{2} {\text{dt}}} = \sqrt {\sigma_{{\tilde{\xi }}}^{2} + \sigma_{{\xi^{{\prime }} }}^{2} + 2\sigma_{{\tilde{\xi }\xi^{{\prime }} }} } = \sqrt {\sigma_{{\tilde{\xi }}}^{2} + \sigma_{{\xi^{{\prime }} }}^{2} } $$

(2)

where $\sigma_{{\tilde{\xi }}}$ and $\sigma_{{\xi^{{\prime }} }}$ are the standard deviations of $\tilde{\xi }\left( t \right)$ and $\xi^{{\prime }} \left( t \right)$ respectively, while $\sigma_{{\tilde{\xi }\xi^{{\prime }} }}$ is the covariance between $\tilde{\xi }\left( t \right)$ and $\xi^{{\prime }} \left( {\text{t}} \right)$ and it is null, as $\tilde{\xi }\left( t \right)$ and $\xi^{{\prime }} \left( {\text{t}} \right)$ are independent.

If $\sigma_{\xi }$ represents the total variability of $\xi \left( t \right)$, the quantities $\sigma_{{\tilde{\xi }}}$ and $\sigma_{{\xi^{{\prime }} }}$ can be used to quantify the periodicity and the randomness of $\xi \left( t \right)$.

If the instantaneous phase $\phi \left( t \right)$ of the motion can be calculated, the phase average $\xi_{ph} \left( \theta \right)$ and its standard deviation $\sigma_{\xi ph} \left( \theta \right)$ can thus be derived, for each value of the phase $\theta$ included in the range $\left[ {0,2\pi } \right)$, by averaging (and computing the standard deviation of) the samples of $\xi \left( t \right)$ characterized by a same value $\theta$ of the phase. This process is indeed the phase average. $\phi \left( t \right)$ can be generated according to the leading component of motion, which is almost harmonic.

Since the motion trajectory must be closed, any phase-averaged motion component $\xi_{ph} \left( \theta \right)$ will be periodic with $\theta$. If $\xi$ is the leading harmonic component of motion, $\xi_{ph} \left( \theta \right)$ will be a sinusoid. If $\xi$ is the secondary component of motion,$\xi_{ph} \left( \theta \right)$ could theoretically be characterized by infinite harmonics n/2π (with $n = 1 \to \infty$) but, actually, in literature VIV cases, it is characterized by two harmonics at the most [18]. Therefore, the quantities $\sigma_{{\tilde{\xi }}}$ and $\sigma_{{\xi^{{\prime }} }}$ can be estimated as follows:

$$ \sigma_{{\tilde{\xi }}} \cong \sqrt {\frac{1}{2\pi }\mathop \smallint \limits_{0}^{2\pi } \xi_{ph}^{2} \left( \theta \right)d\theta } $$

(3)

$$ \sigma_{{\xi^{{\prime }} }} \cong \sqrt {\frac{1}{2\pi }\mathop \smallint \limits_{0}^{2\pi } \sigma_{\xi ph}^{2} \left( \theta \right)d\theta } $$

(4)

In order to quantify the harmonics that compose $\xi_{ph} \left( \theta \right)$, the Fourier coefficients are calculated:

$$ {\mathcal{A}}_{\xi n} = \frac{1}{\pi }\mathop \smallint \limits_{0}^{2\pi } \xi_{ph} \left( \theta \right)e^{{ - j\frac{n}{2\pi }\theta }} d\theta $$

(5)

where ${\mathcal{A}}_{\xi n}$ is the complex amplitude of the nth harmonic and $j = \sqrt { - 1}$. If $\xi$ is the main component of motion, $n$ = 1. If $\xi$ is the secondary component of motion, $n$ = 1 and/or 2 [18].

Let one consider that the motion components are simultaneously sampled at a constant frequency $f_{s}$. The instantaneous phase $\phi \left( t \right)$ of the motion is calculated on the basis of the main component of motion, which is usually harmonic. Similarly to the approach reported in Negri et al. [35], the reference signal $r\left( t \right)$ is generated by applying a low-pass filter on the harmonic component of motion, and zero-crossing instants $t_{0,i}$ with positive slope are detected in $r\left( t \right)$ through linear interpolation (with $i = 1,2, \ldots ,N$ where $N$ is the total number of zero-crossing points detected). The signal $r\left( t \right)$ is thus divided into $N - 1$ periods. The average oscillation period is defined as the average difference between consecutive zero-crossing instants:

$$ T_{0} = \frac{1}{N - 1}\mathop \sum \limits_{i = 1}^{N - 1} \left( {t_{0,i + 1} - t_{0,i} } \right) $$

(6)

The period variability can be described by the standard deviation $\sigma_{T}$ of the periods detected in $r\left( t \right)$:

$$ \sigma_{T} = \sqrt {\frac{1}{N - 2}\mathop \sum \limits_{i = 1}^{N - 1} \left( {t_{0,i + 1} - t_{0,i} - T_{0} } \right)^{2} } $$

(7)

Let one introduce the reference period $T$, approximation of the period $T_{0}$ such as it is a multiple of the sampling time $\frac{1}{{f_{s} }}$ or a multiple thereof (Eq. 8):

$$ T \equiv [\frac{{T_{0} }}{{M/f_{s} }} + 0.5]\frac{M}{{f_{s} }} $$

(8)

where $\left\lfloor \ldots \right\rfloor$ is the floor function and $M$ is an integer number greater than or equal to 1. The approximation error will be $\left( {T - T_{0} } \right) \le \frac{M}{{2f_{s} }}$.

In this work, two different techniques of assigning the instantaneous phase $\phi \left( t \right)$ to $r\left( t \right)$ are considered, depending on whether the time instants $t$ are non-dimensionalized by a unique period, which is the reference period $T$ [Eq. (9)], or by the current period $t_{0,i + 1} - t_{0,i}$ [Eq. (10)], variable from cycle to cycle. Even if $T_{0}$ is the best estimate of the oscillation period, its approximated value $T$ allows the optimization of the phase average when the phase is defined according to Eq. (9):

$$ \phi_{c} \left( t \right) = 2\pi \frac{{t - t_{0,i} }}{T},\quad t_{0,i} \le t < t_{0,i + 1} , \quad i = 1,2, \ldots ,N - 1 $$

(9)

or

$$ \phi_{s} \left( t \right) = 2\pi \frac{{t - t_{0,i} }}{{t_{0,i + 1} - t_{0,i} }}, \quad t_{0,i} \le t < t_{0,i + 1} , \quad i = 1,2, \ldots ,N - 1 $$

(10)

where $t_{0,i}$ are the $N$ instants of zero-crossing of the reference signal $r\left( t \right)$.

$\phi_{c} \left( t \right)$ and $\phi_{s} \left( t \right)$ are sine phases, since the 0-phase value corresponds to the positive zero-crossing. $\phi_{c} \left( t \right)$ and $\phi_{s} \left( t \right)$ are periodic functions made of line segments (indeed, they are linear-interpolating phases).

The first phase assignment [Eq. (9)] has the advantage that all time instants are non-dimensionalized for the same quantity, therefore the resulting phase spacing of the samples is the same for all the identified periods and is equal to $\frac{2\pi }{{T \cdot f_{s} }}$. On the other hand, it has the disadvantage that the phase within a cycle does not exactly range between 0 and $2\pi$, but it ranges from 0 to a value smaller or larger than $2\pi$, whether the current period $t_{0,i + 1} - t_{0,i}$ is shorter or longer than the reference period $T$. This issue is as great as the variability of the period $\sigma_{T}$ is large. If the phase is assigned according to Eq. (10), on the contrary the phase ranges from 0 to $2\pi$ in all the cycles but the phase spacing of the samples is different from cycle to cycle, since the time instants are non-dimensionalized on the basis of different periods.

In order to perform the phase average of the samples, the phase domain $\left[ {0,2\pi } \right)$ is subdivided into $K$ equally-spaced bins, whose center coordinates are:

$$ \theta_{K} = \frac{{\left( {k - 0.5} \right)}}{K}2\pi ,\quad k = 1,2, \ldots ,K $$

(11)

with $K$ being an integer number.

If the phase is calculated according $\phi_{c} \left( t \right)$, the phase bins can be arranged so that every bin includes the same number $M$ of samples from each period, except in the last bin(s), where, due to period variability, samples from the shorter periods are missing. This is possible only if the width of the phase bin $\frac{2\pi }{K}$ is a multiple of $\frac{2\pi }{{T \cdot f_{s} }}$, and, considering that $K$ must be an integer, if $T$ is a multiple of $\frac{M}{{f_{s} }}$ [see Eq. (8)]. Hence, the optimum number $K$ of phase bins results:

$$ K = \frac{{T \cdot f_{s} }}{M} $$

(12)

The smaller $M$, the larger $K$, which means that the lower the smoothing, the finer the discretization of the phase average. If $M$ = 1, the phase average has the same temporal resolution of the original signal.

If phase is calculated according $\phi_{s} \left( t \right)$, the arrangement of the phase bins cannot be optimized to evenly include the samples from the different periods.

The $k$th values of the phase average $\xi_{ph,k}$ and phase standard deviation $\sigma_{\xi ph,k}$ of the generic motion component are calculated as:

$$ \xi_{ph,k} = \frac{1}{{J_{k} }}\mathop \sum \limits_{j = 1}^{{J_{k} }} \left( {\xi \left( {t_{k,j} } \right) - \bar{\xi }} \right), \quad t_{k,j} = \left\{ {t |\left( {\theta_{k} - \frac{\pi }{K} \le \phi \left( t \right) < \theta_{k} + \frac{\pi }{K}} \right)} \right\} $$

(13)

$$ \sigma_{\xi ph,k} = \sqrt {\frac{1}{{J_{k} - 1}}\mathop \sum \limits_{j = 1}^{{J_{k} }} \left( {\xi \left( {t_{k,j} } \right) - \xi_{ph,k} } \right)^{2} } ,\quad t_{k,j} = \left\{ {t |\left( {\theta_{k} - \frac{\pi }{K} \le \phi \left( t \right) < \theta_{k} + \frac{\pi }{K}} \right)} \right\} $$

(14)

where $t_{k,j}$ are the time instants such that their instantaneous phase is included in the kth phase bin, and $J_{k}$ is the number of these time instants.

The type of phase assignment [Eq. (9) or (10)] significantly affects the phase average. Two properties are here considered for evaluating the accuracy of the phase average:

• Smoothness: the smoothness of the original oscillation periods is conserved in their phase average.

• Periodicity: the last element $\xi_{ph,K}$ well approaches the first element $\xi_{ph,1}$ of the phase average.

The phase average can be usefully interpreted as a moving average of the phase-sorted samples, which are obtained by reorganizing the samples of $\xi \left( t \right)$ according to their (increasing) phase ($\phi_{c}$ or $\phi_{s}$). $\xi \left( \phi \right)$ can be interpreted as the phase-sorted samples. Generally, $\xi \left( \phi \right)$ is a multivalued “function”, since more than one sample of $\xi \left( t \right)$ can be characterized by the same value of phase. Moreover, the elements of $\xi \left( \phi \right)$ are not equally-spaced.

If the phase is assigned according to a constant period [Eq. (9)], the sequence with which samples from different periods follow one another in $\xi \left( \phi \right)$ is constant and is repeated every $\frac{2\pi }{{T \cdot f_{s} }}$ radians: $N - 1$ samples (one sample from each period) are included, with the same arrangement, inside every phase interval of $\frac{2\pi }{{T \cdot f_{s} }}$ radians, as long as $\phi < \frac{{T_{min} }}{T}2\pi$, where $T_{min}$ is the shortest period detected in $r\left( t \right)$. Therefore, when performing the moving average of $\xi \left( \phi \right)$, a bin width of $M\frac{2\pi }{{T \cdot f_{s} }}$ [see Eq. (12)] ensures that the smoothness of the cycles is also preserved in their average $\xi_{ph,k}$, for $\frac{k + 0.5}{K} < \frac{{T_{min} }}{T}$. The number of samples $J_{k}$ in the kth phase bin will be constant and equal to $M\left( {N - 1} \right)$, as long as $\frac{k + 0.5}{K} < \frac{{T_{min} }}{T}$. In the last bin(s), that is for $\frac{k + 0.5}{K} \ge \frac{{T_{min} }}{T}$, the number of samples $J_{k}$ could be lower and neither the smoothness nor the periodicity of the phase average is guaranteed.

Instead, defining the phase according to the current period [Eq. (10)] has the advantage that all the periods have the same duration in terms of phase (2π), and thus the periodicity of the resulting phase average is guaranteed. On the other hand, samples from different periods do not occur with a repeating sequence in $\xi \left( \phi \right)$, hence the arrangement of the phase bins cannot be optimized to evenly group the elements of $\xi \left( \phi \right)$. Generally, the samples from each period will not be uniformly distributed in the phase bins. Therefore, the smoothness of the originally sampled periods is not conserved in their average.

In order to exploit the benefits of the two methods and to obtain both a periodic and smooth phase average, an additional method is implemented: the samples of every period are interpolated on $K$ equally-spaced points, corresponding to the phase values of $\theta_{K}$ [Eq. (12)], with phase being calculated according to the current period [Eq. (10)]. In this way, the new interpolated samples are uniformly distributed throughout the phase bins: $N - 1$ samples (one from each period) fall in every bin center. The interpolation instant with phase $\theta_{K,k}$ in the ith period is defined as $t_{k,i} = \frac{{\theta_{K,k} }}{2\pi }\left( {t_{0,i + 1} - t_{0,i} } \right) + t_{0,i}$. Let ${\check{\xi}}_{k,i}$ be the interpolated value of $\xi \left( t \right)$ at the instant $t_{k,i}$. Then, the kth elements of the phase average and phase standard deviation are calculated as [Eqs. (15) and (16), respectively]:

$$ \xi_{ph,k} = \frac{1}{N - 1}\mathop \sum \limits_{i = 1}^{N - 1} \left({{\check{\xi}}_{k,i} - \bar{\xi}} \right) $$

(15)

$$ \sigma_{\xi ph,k} = \sqrt {\frac{1}{N - 2}\mathop \sum \limits_{i = 1}^{N - 1} \left({{\check{\xi}}_{k,i} - \xi_{ph,k}} \right)^{2}} $$

(16)

4 Results

The results are subdivided into two sections: (1) the detailed application of the presented method to a single VIV case of a tethered sphere, and (2) the analysis of the VIVs of a tethered sphere with varying reduced velocity.

4.1 Application of the phase average to a single VIV case

The case analyzed in this section was characterized by $D$ = 0.060 m, L* = 11.7, m* = 1.38, $\zeta$ = 0.03, U* = 7.7, and Re = 7000. The sampling frequency $f_{s}$ was 10 Hz and 1499 samples were acquired. In Fig. 2, the streamwise (a) and crossflow (b) coordinates of the sphere are depicted. Figure 2b reports the reference signal $r\left( t \right)$ for the phase average, obtained by low-pass filtering the crossflow signal $y\left( t \right)$ (which is the harmonic component of motion), and the zero-crossing points detected in $r\left( t \right)$. The Fourier transform of the motion components is depicted in Fig. 3; the cut-off frequency for the filter used to obtain $r\left( t \right)$ from $y\left( t \right)$ is set at twice the main oscillation frequency, f, and is indicated by the dashed line. In Negri et al. [35], the value of the cut-off frequency slightly affected the results, once that the frequency peak was included in the band. The spectrum of the crossflow component clearly shows a frequency peak $f \approx 0.26$ Hz, while the spectrum of the streamwise component x(t) is noisier, although a frequency peak approximately at $2f$ is detectable.

Fig. 2

Time history of the sphere positions non-dimensionalized against the sphere diameter, D: a streamwise position; b crossflow position

Fig. 3

Discrete Fourier Transform (DFT) of the non-dimensionalized a streamwise and b crossflow component of Fig. 2. The cut-off frequency for generating the reference signal from the y-component of motion is 2f

The number of periods detected in the reference signal $r\left( t \right)$ is $N - 1 = 39$. The average period and period variability are $T_{0} = 3.806$ s and $\sigma_{T} = 0.194$ s, respectively. The variability $\sigma_{T}$ is about the 5% of $T_{0}$. As showed by Fig. 3b, the periodicity of the crossflow component of motion is quite high. In Fig. 4, the duration $t_{0,i + 1} - t_{0,i}$ of the $N - 1$ periods detected in $r\left( t \right)$ is showed.

Fig. 4

Period duration $t_{0,i + 1} - t_{0,i}$ throughout the signal (plotted at $t = \left( {t_{0,i} + t_{0,i + 1} } \right)/2$). The average period $T_{0}$ and the standard deviation $\sigma_{T}$ are also shown

According to Eq. (8) and Eq. (12), adopting M = 1 (maximum resolution of the phase average that allows the constancy of $J_{k}$ through the phase bins), $T$ = 3.80 s and K = 38 are obtained. In Fig. 5, the instantaneous phase $\phi \left( t \right)$ of $y\left( t \right)$, calculated according Eq. (9) ($\phi_{c} \left( t \right)$) and Eq. (10) ($\phi_{s} \left( t \right)$), is depicted. $\phi_{s} \left( t \right)$ ranges from 0 to 2π in each period, while $\phi_{c} \left( t \right)$ starts at 0 and generally does not end at 2π in each period.

Fig. 5

Instantaneous phase $\phi \left( t \right)$ of $y\left( t \right)$, calculated according to Eqs. (9) and (10)

Once the instantaneous phase $\phi \left( t \right)$ is known, the samples of $y\left( t \right)$ are reorganized according to increasing phase, in order to obtain $y\left( \phi \right)$. Figure 6 shows $y\left( \phi \right)$ using Eqs. (9) and (10). The limits of the phase bins are also reported.

Fig. 6

Samples of crossflow motion $y\left( \phi \right)$ sorted according to their phase $y\left( \phi \right)$: a instantaneous phase calculated according to a constant period [Eq. (9)]; b instantaneous phase calculated according to the current period [Eq. (10)]. The limits of the phase bins are indicated by the dashed lines

Figure 7 shows the number $J_{k}$ of samples of each phase bin. If the phase is calculated according to a constant period [Eq. (9)], $J_{k}$ is constant except for the last bins, where it decreases due to the different duration of the periods (Fig. 7). On the contrary, if phase is calculated according to the current period [Eq. (10)], $J_{k}$ is not constant and there is not a decrease of samples in the last bins.

Fig. 7

Number of samples $J_{k}$ inside the phase bins, depending on whether the instantaneous phase is calculated according to a constant period [Eq. (9)] or to the current period [Eq. (10)]. For the case of phase calculated according to Eq. (9), additional bins beyond 2π are considered

In Fig. 8, the interpolated periods are reported. The interpolant used is cubic, but it has been verified that linear interpolation gives almost identical results.

Fig. 8

Phase-sorted [according to phase of Eq. (10)] and interpolated samples of the crossflow motion $y\left( \phi \right)$

In Fig. 9 and 10, the phase average and standard deviation of the crossflow and streamwise component are reported, using the three techniques previously explained: (1) instantaneous phase calculated according to a reference period [Eq. (9)], moving average and standard deviation [Eqs. (13) and (14)], (2) instantaneous phase calculated according to the current period [Eq. (10)], moving average and standard deviation [Eqs. (13) and (14)], (3) instantaneous phase calculated according to the current period [Eq. (10)], interpolation of the periods, moving average and standard deviation [Eqs. (15) and (16)].

Fig. 9

Phase-averaged crossflow motion: a bin mean; b bin standard deviation

Fig. 10

Phase-averaged streamwise motion: a bin mean; b bin standard deviation

It can be seen that the first technique preserves the smoothness in the phase average (except that in the last bins) but not the periodicity. On the contrary, the second technique guarantees the periodicity but does not preserve the smoothness. The third technique, which combines the benefits of both the previous two techniques, guarantees both periodicity and smoothness. This is also evident in the 2D trajectories (Fig. 11): the trajectory calculated with the first technique is smooth but not closed; the trajectory calculated with the second technique is indented but closed; the trajectory calculated with the third technique is smooth and closed.

Fig. 11

Phase-averaged xy trajectory. The arrows indicate the trajectory direction. The flow is from left to right

The oscillation of the bin standard deviation through the period (Figs. 9, 10) indicates that the amplitude of the random component (but not the random component itself) is correlated to the periodic component. For the crossflow motion (which is the leading motion component in this case), the amplitude of the random component is higher where the module of the periodic component is higher, and thus it accomplishes two oscillation cycles within the period. This is intuitive for an almost harmonic signal. For the streamwise motion (which is the secondary motion component in this case), the amplitude of the random component is almost in phase with the periodic component.

Finally, the amplitude and phase of each motion component were obtained by calculating the first Fourier coefficients of $x_{ph,k}$ and $y_{ph,k}$ (two harmonics for $x_{ph,k}$ and one harmonic for $y_{ph,k}$, the two harmonics corresponding to the frequencies $1/\left( {2\pi } \right)$ and $2/\left( {2\pi } \right)$). In Tables 1 and 2, the amplitudes of the first Fourier coefficients are reported for the streamwise and crossflow motion, respectively (the number of harmonics reported exceeds by two the effective number of harmonics of $x_{ph,k}$ and $y_{ph,k}$). Moreover, the standard error (SE) and the coefficient of determination, R², of the Inverse Discrete Fourier Transform (IDFT) of the phase average are reported for the increasing number of harmonics considered. It can be seen that the streamwise motion is completely described by two harmonics, while the crossflow motion is completely described by one harmonic.

Table 1 Amplitude of the harmonics of the streamwise motion, A_xn*, calculated through DFT; Standard Error (SE), and determination coefficient R² of the IDFT of the phase average for an increasing number of harmonics considered

Table 2 Amplitude of the harmonics of the crossflow motion, A_yn*, calculated through DFT; Standard Error (SE) and determination coefficient R² of the IDFT of the phase average for an increasing number of harmonics considered

The original data $x\left( t \right)$ and $y\left( t \right)$ were also analyzed with the other methods used in the VIV literature and indicated in the introduction as 1, 2, and 3. Table 3 shows the following quantities detected by each method:

Table 3 Comparison among the different methods of trajectory analysis

• main oscillation frequency $f$ (main component of motion);

• crossflow first harmonic amplitude $A_{y}^{*}$;

• streamwise amplitudes $A_{x1}^{*}$ (first harmonic) and $A_{x2}^{*}$(second harmonic);

• phase lags $\phi_{xy1}$, between the first harmonics of x and y, and $\phi_{xy2}$, between the first harmonic of y and the second harmonic of x, respectively.

In Fig. 12, the trajectory is modeled using the methods reported in Table 3. The comparison between the modeled trajectory and the original data is proposed through two different plots: on the left hand side of the figure, the original data are plotted as time-resolved trajectories; on the right hand side, the original data are plotted as frequency distribution on xy intervals, with width 0.0134D and 0.0432D in the $x$ and $y$ directions, respectively. For an overall estimation of the methods’ accuracy, the following parameter was used and included in Table 3:

$$ n^{{\prime }} \equiv \frac{{\mathop \sum \nolimits_{q = 1}^{Q} n_{q} }}{s} $$

(17)

where $n_{q}$ is the number of samples belonging to the qth xy interval crossed by the modeled trajectory (out of a total number Q of intervals, which is different for each method considered) and $s$ is the length of the modeled trajectory (expressed in sphere diameters).

Fig. 12

Trajectory modeling with different methods. a The original data are represented as xy trajectory; b the original data are represented as frequency distribution (the color is linked to the number of samples in each xy interval)

$n^{\prime}$ represents the number of samples per unit length of modeled trajectory: the larger $n^{\prime}$ is, the more the modeled trajectory traces the original positions of the sphere. It must be noted that the parameter $n^{\prime}$ is dependent on the dimension of the xy intervals; moreover, this estimation of accuracy cannot account for the orientation of the trajectory.

The first row of Table 3 reports the amplitudes calculated through the standard deviation. Since the spectrum of y is concentrated on the main peak (Fig. 3a), the y-amplitude calculated this way should be indicative of the mean amplitude. On the contrary, the spectrum of x (Fig. 3b) is very dispersed in the range 0–0.6 Hz, therefore the x-amplitude indicates only the mean variability and not the amplitude of any of the harmonics. In the second row of Table 3, the y-amplitude is calculated as mean peak-to-peak distance, while the x-amplitude cannot be calculated because x is not harmonic.

The first method, DFT, is the most straightforward one to calculate the main oscillation frequency, f, but it is also the least accurate in reconstructing the sphere trajectory since the amplitude and period duration of both motion components vary with time (see Fig. 2 for the amplitudes and Fig. 4 for the period duration of the crossflow motion component).

The second method, used in Kang and Jia [28], is the regression of $x\left( t \right)$ and $y\left( t \right)$ on a limited Fourier series (5 harmonics) with fixed frequencies (multiple of f) and unknown amplitudes and phases. Amplitude coefficients calculated with this method are 0.0175, 0.0238, 0.0005, 0.0004, and 0.0006 for $x$ and 0.4090, 0.0026, 0.0007, 0.0021, and 0.0015 for $y$, respectively. It is clear that the significant amplitudes are the first two for $x$ and only the first one for $y$. The results of this method are similar to those of the DFT method: the y-amplitudes are underestimated and the modeled trajectory does not resemble the real motion (Fig. 12). The rotation direction of the trajectory is clockwise, in agreement with the literature data [18].

The third method to characterize the $xy$ trajectory is to identify the synchronization between x and y (see [30]). This can be done by (1) calculating the instantaneous phases $\phi_{x} \left( t \right)$ and $\phi_{y} \left( t \right)$ of $x\left( t \right)$ and $y\left( t \right)$, (2) hypothesizing a frequency ratio $f_{x} /f_{y} \equiv m/n$ (with $m$ and $n$ integer numbers), (3) calculating the instantaneous phase lag between the two motion components $\phi_{xy,m:n} \left( t \right) = n\phi_{x} \left( t \right) - m\phi_{y} \left( t \right)$, and (4) calculating the mean phase lag from the frequency distribution of the instantaneous phase lag. Two techniques for calculating the instantaneous phase of a signal are here considered: through the analytical signal (Hilbert transform) of $x\left( t \right)$ and $y\left( t \right)$ [30] and through the linear interpolating phase [31]. The latter was here calculated by assigning the phases $\pi /2$ and $3\pi /2$ to the relative maxima and minima respectively (i.e. sine phase) of each motion component, and then linearly interpolating between these points. Here, the “average” phase lag is taken as the mode of the frequency distribution of the instantaneous phase lag. The application of these two techniques allowed calculating $\phi_{xy2}$ (although with a certain degree of uncertainty) but not $\phi_{xy1}$ (the frequency histogram does not show a recognizable peak). $\phi_{xy2}$ calculated with the two techniques are very different: the analytical signal gives a counterclockwise trajectory (Fig. 12a), which is not expected for the investigated value $U^{*}$ = 7.7. With both techniques, the amplitudes derived from the standard deviations are used to represent the trajectory.

The method developed here allows the determination of all the parameters of the average trajectory, such as the DFT method. The mean period is similar to the one detected through the DFT. The phase lag $\phi_{xy2}$ is very close to that estimated with the linear interpolating phase. The mean y-amplitude $A_{y}^{*}$ is similar to the average peak-to-peak amplitude. Compared to the other methods found in the literature, the present method of phase average better calculates all properties and reconstructs the mean sphere trajectory together with its rotation direction, as shown in Fig. 12. This is confirmed by high value of the parameter $n^{\prime}$ in Table 3.

4.2 VIV of a tethered sphere with varying reduced velocity

In order to study the phenomenon of the VIVs of a tethered sphere, a sufficiently wide range of reduced velocities was investigated. The tests were performed changing the tether length, to modify the natural frequency of the system, and varying the mean streamwise velocity (0.09 < U < 0.4 m/s).

The non-dimensionalized sphere mass was equal to 1.38, while the damping ratio was included in the range of 0.023 < $\zeta$ < 0.026. This led to a correspondent range 0.043–0.049 of the mass-damping parameters, $\left( {m* + C_{A} } \right)\zeta$, where C_A is the potential added mass coefficient equal to 0.5 for a sphere. According to the Griffin plot, the mass-damping parameter affects the crossflow amplitude motion if it is larger than 0.02 and thus the motion amplitude was expected to be slightly reduced [10]. The dynamic response of the sphere in the crossflow direction, $A_{y}^{*} ,$ is influenced by the set of non-dimensional groups defined in Table 4 and is expressable with the relationship of Khalak and Williamson [36]:

$$ A_{y}^{*} = \frac{3}{{32\pi^{2} }}\frac{{C_{y} sin\varPhi_{y} }}{{\left( {m^{*} + C_{A} } \right)\zeta }}\left( {\frac{{U^{*} }}{{f^{*} }}} \right)^{2} f^{*} $$

(18)

where C_y is the total crossflow force coefficient, $\varPhi_{y}$ the phase difference between the sphere crossflow displacement and the total crossflow force, and f* = f/f_n the frequency ratio equal to:

Table 4 Non-dimensional groups affecting the crossflow dynamic response

$$ f^{*} = \sqrt {\frac{{\left( {m^{*} + C_{A} } \right)}}{{\left( {m^{*} + C_{EA} } \right)}}} $$

(19)

Here, C_EA is an “effective” added mass coefficient that includes an apparent effect due to the total crossflow fluid force in phase with the body acceleration ($C_{y} \cos \varPhi_{y}$):

$$ C_{EA} = \frac{3}{{16\pi^{2} A_{y}^{*} }}\left( {\frac{{U^{*} }}{{f^{*} }}} \right)^{2} C_{y} cos\varPhi_{y} $$

(20)

Figure 13 displays the trend of the crossflow amplitude oscillations with varying U*.

Fig. 13

Crossflow oscillation amplitude versus U*. The standard deviation of the random component $\sigma_{{y^{{\prime }} }} /D$ [Eq. (4)] is reported in the bands

The amplitude response of the sphere in the current study is similar to the trend of a heavy sphere with m* = 2.8 [9]. At U* ≈ 5, Mode I is detected, which corresponds to the vortex formation frequency lying close to the natural frequency of the tethered body (Fig. 14). For U* > 7, the sphere starts having oscillations higher than the amplitudes corresponding to the first response, until it reaches the maximum crossflow amplitude, A*_y≈ 0.9, corresponding to Mode II. This value is equal to the saturation value 0.9 found by Govardhan and Williamson [10] for values of mass–damping lower than 0.02.

Fig. 14

Crossflow oscillation frequency versus U*. The variability of the oscillation frequency $f^{*} \sigma_{T} /T$ is reported in the bands

The transition between the modes is continuous and gradual with U* and no separation is observed, unlike the case of a tethered sphere with a very low mass ratio, where the two vibration modes are distinct by a desynchronized region. The randomness of vibrations increases with increasing U*, as shown by the variability (represented with the bands) of oscillation amplitude (Fig. 13) and frequency (Fig. 14), similarly to van Hout et al. [2] and Sareen et al. [25].

Figure 14 shows that, for values of U* around 5, the oscillation frequency is very close to the natural frequency of the body. With increasing U*, the oscillation frequency increases, as does the variability of the frequency, expressed by the bands.

As observed in Fig. 13, the transition between Mode I and II is continuous with U*, and thus more difficult to detect than from the plot of the phase differences between the force and the sphere displacement signals.

Indeed, Fig. 15 displays the variations of $\varPhi_{y}$ and $\varPhi_{y,v}$ (phase difference between the sphere displacement and the vortex force) versus the reduced velocity. In the same figure, the trends of the total crossflow force coefficient, C_y, and the vortex force coefficient, C_y,v, are also reported. The phase of the vortex force, $\varPhi_{y,v}$, is equal to 90° when Mode I occurs and exceeds 90° as the Mode II response reaches the peak amplitude. The phase of the total crossflow force, however, remains constant at values around zero in Mode I and crosses 90° in Mode II. This difference might be due to a change in the shape of vortex structures in the wake of the vibrating sphere [10].

Fig. 15

Variation of the total force coefficient (a) and phase (c) and vortex force coefficient (b) and phase (d) with U*

Figure 16 presents the sphere trajectories on the transverse plane for several values of U*, on which the trajectories obtained through the phase average are overlapped. For U* > 5, the second harmonic appears in the streamwise component and the parameter which better represents the trajectories’ shape is the phase lag, $\phi_{xy2}$, between the streamwise and crossflow motions showed in Fig. 17. As evident from the figure, the phase lag increases with U* and is included in the range $\pi \le \phi_{xy2} \le 7\pi /4$. $\phi_{xy2}$ seems to have two plateaus, one in Mode I ($\phi_{xy2} \approx \pi /6$) and one in Mode II ($\phi_{xy2} \approx 5\pi /3$). According to the definition of Dahl et al. [37], Fig. 17 also indicates clockwise trajectories ($\phi_{xy2} < 3\pi /2$) in Mode I and in the transition region, and counterclockwise ones ($\phi_{xy2} > 3\pi /2$) in Mode II. Figure 16 displays how the clockwise trajectories are more chaotic than the counterclockwise ones. An explanation given in Dahl et al. [37] for cylinder VIVs is that, in counterclockwise trajectories, the body goes downstream after the vortices are shed, staying closer to them as they go downstream too, and hence exploiting the suction force more. As well as the rotation direction of the trajectory, the phase $\phi_{xy2}$ also describes the shape of the trajectory [18]. In Mode I and Mode II, the trajectory is 8-shaped. Between the two modes, approximately at U* ≈ 8.5, $\phi_{xy2}$ crosses $3\pi /2$ and the trajectory is crescent-type. In Mode I, the trajectory is characterized by central symmetry. For U* > 5, the tips of the trajectory are downstream the center of the trajectory.

Fig. 16

Experimental trajectories in the xy plane with the phase-averaged trajectories and the IDFT of the phase-averaged trajectories superimposed for: a U* = 5.68, b U* = 7.64, c U* = 8.52, d U* = 9.62, e U* = 12.72, f U* = 14.09, g U* = 15.96, h U* = 25.33

Fig. 17

Phase lag between the streamwise and crossflow motion component

In Fig. 18, the periodicity of the streamwise and crossflow motion component is shown and was calculated as the ratio between the amplitude of the periodic contribution (according to Eq. (3)) and the total amplitude. Since the periodicity of the crossflow component is always high, the periodicity of the streamwise component can be also read as the synchronization degree between the two motion components. In the investigated range, the periodicity of the crossflow components decreases as U* increases, while the periodicity of the streamwise components has a maximum corresponding to U* ≈ 8.5. In Mode II, both the crossflow and streamwise periodicity decreases as U* increases.

Fig. 18

Periodicity of the streamwise and crossflow motion component

In the present results, the sphere shows significant oscillation amplitude up to a normalized velocity (U*/f*)St ≈ 2.7 (where St = f_voD/U and f_vo is the vortex shedding frequency with a stationary sphere), therefore the end of Mode II is not detected. Govardhan and Williamson [10] (case m* = 2.8, $\left( {m* + C_{A} } \right)\zeta$ = 0.029) and Sareen et al. [25] (case m* = 7.8, $\left( {m* + C_{A} } \right)\zeta = 0.017$) found that Mode II extended up to (U*/f*)St ≈ 2.2. After this value, Sareen et al. [25] found a plateau with still significant oscillation amplitude, while Govardhan and Williamson [10] observed a strong decrease of the oscillation amplitude before approaching Mode III. Since the mass-damping parameter in this work is higher, $\left( {m* + C_{A} } \right)\zeta = 0.043 \div 0.049$, the persisting Mode II may be linked to the low value of the sphere mass m* = 1.38.

5 Conclusions

In this work, a phase average method was developed in order to analyze the motion of a tethered sphere with m* = 1.38, which is characterized by both periodicity and randomness. The method was applied to experimental VIVs considering a wide range of reduced velocities 4.96 ≤ U* ≤ 25.33 (Reynolds numbers 5x10³ ≤ Re ≤ 2.5x10⁴). Samples acquired at a constant frequency were analyzed. The component of motion showing periodic behavior (crossflow component) was used as reference phase and the trigger points were precisely detected through interpolation, dividing periods throughout the time history of that motion component. Due to the variability of the oscillation period within the crossflow component, the phase was assigned to each sample with respect to the current oscillation period. The resulting non-uniform distribution of samples against the phase was overcome through an interpolation of the samples on the basis of a phase vector of equally-spaced points.

The results showed how the mean and standard deviation of the phase average identified and quantified the periodic and random part, respectively. The Fourier transform of the periodic part revealed that the streamwise motion component was characterized by two harmonics (crossflow frequency and its double). The amplitude (standard deviation) of the random part was observed oscillating twice the crossflow frequency, both in the streamwise and crossflow direction, while it increased with increasing U* for the crossflow motion only. Analogously, the variability of the oscillation period increased with U*.

The method developed here proved to be efficient in the presence of a periodic component of motion and when the motion is stationary, being abler to extract the periodic trajectory pattern than the other methods used in the literature, especially with VIV trajectories characterized by a level of randomness. Further work is needed to adapt the present method to the analysis of motions characterized by synchronization changes.

Additional results were obtained regarding the dynamic response of the sphere. In particular, the crossflow oscillation amplitude closely followed the trend reported in the literature for similar mass ratios showing the two response modes, Mode I and Mode II. A saturation peak of about 0.9D, equal to the one found by Govardhan and Williamson [10], was also observed. The frequency response was close to unity in the resonance condition and increased with U*. For higher values of U*, significant variability in amplitude and frequency of vibration highlighted a less periodic motion of the sphere. The transition between Mode I and II was more evident in the plot of the vortex phase corresponding to a shift in the timing of vortex formation. The shape and rotation direction of the trajectory are linked to the phase lag between the crossflow and streamwise motion (2^nd harmonic), $\phi_{xy2}$, which increased with U*, at least in the range investigated here (U* = 5 ÷ 25). The sphere trajectories presented the classical 8-shape with a clockwise direction of rotation in Mode I ($\phi_{xy2} < 3\pi /2$). An increase of U* transformed the trajectory into a crescent-like shape (U* ≈ 8.5) as $\phi_{xy2}$ crossed $3\pi /2$; with higher U* (Mode II), $\phi_{xy2}$ increased too and the trajectory became 8-shaped again, but with counterclockwise rotation. The highest degree of the trajectory periodicity was found at U* ≈ 8.5. As U* increased within Mode II, the trajectory was more and more characterized by a higher level of randomness.

Appendix: a method for calculating the three-dimensional position of the sphere using a single framing

The motion of the sphere was recorded by a camera placed under the sphere, with its optical axis being vertical. The sphere had a marker at its bottom, aligned along the tether direction. During the sphere motion, the marker moved on a spherical surface, hence changing its distance from the camera. This caused a change in the resolution R [px/m] of the image of the current plane, where the marker lied. In order to obtain the true tridimensional position of the marker from each frame, it was thus necessary to perform some correction on the raw position of the marker in the image.

1.1 Perspective correction

The perspective scheme is considered in Fig. 19. The tether pivot is in C. The sphere is in a generic position and its marker is in point P with coordinates $\left( {x_{P} ,y_{P} ,z_{P} } \right).$ The distance between the tether pivot and the sphere marker is ℓ. Three planes are indicated in the scheme:

Fig. 19

Perspective scheme

• the sensor plane in the camera, where the image is created. The width of the sensor is $\varvec{w} = \varvec{n}_{\varvec{p}} \varvec{l}_{\varvec{p}}$, where $\varvec{n}_{\varvec{p}}$ and $\varvec{l}_{\varvec{p}}$ are the number and dimension of pixels along the $\varvec{x}$ direction, respectively;

• the reference plane, where the marker lies when the sphere is at rest (point O). The image of this plane has a resolution R₀ [px/m];

• the current plane, where the marker lies when the sphere moves (point P). The image of this plane has a resolution R(z) < R₀.

The point F is the focus of the camera. For now, the refraction of the optic rays (due to the different mediums they cross) may be ignored. Two coordinate systems are considered: the absolute reference system, whose origin and horizontal axis x lie in the reference plane, and the image reference system, whose origin and horizontal axis x’ lie in the sensor plane. Note that here, for convenience, the origin of the absolute reference system is different from the one in Fig. 1.

The coordinate $x_{{P^{\prime}}}^{{\prime }}$ of the marker in the image can be expressed as:

$$ x_{{P^{\prime}}}^{{\prime }} = x_{{O^{\prime}}}^{{\prime }} + O^{{\prime }} A^{{\prime }} + A^{{\prime }} P^{{\prime }} $$

(21)

where $x_{{O^{\prime}}}^{{\prime }}$ is the coordinate (in the image) of the marker at rest, which can be known by taking a picture when the sphere is still. The lengths of the segments ${\text{O}}^{{\prime }} {\text{A}}^{{\prime }}$ and ${\text{A}}^{{\prime }} {\text{P}}^{{\prime }}$ are easily calculable (in the subsequent equations, the resolution $R$ is a function of $z$, even though it is not indicated):

$$ O^{\prime}A^{\prime} = R \cdot O^{*} A = R\left( {O^{*} M - O\tilde{M}} \right) = R \cdot \left( {\frac{{O^{\prime}M^{\prime}}}{R} - \frac{{O^{\prime}M^{\prime}}}{{R_{0} }}} \right) = O^{\prime}M^{\prime}\left( {1 - \frac{R}{{R_{0} }}} \right) = \left( {x_{{M^{\prime}}}^{{\prime }} - x_{{O^{\prime}}}^{{\prime }} } \right)\left( {1 - \frac{R}{{R_{0} }}} \right) $$

(22)

$$ A^{\prime}P^{\prime} = R \cdot AP = R \cdot \left( {x_{P} - x_{O} } \right) $$

(23)

where $x_{O} = x_{{O^{\prime}}}^{{\prime }} /R_{0}$. The resolution $R_{0}$ of the image of the reference plane can be calculated from the picture acquired with the sphere at rest.

Finally, it is possible to make the real coordinate of the marker explicit:

$$ x_{P} - x_{O} = \frac{1}{{R\left( {z_{P} } \right)}}\left( {x_{{P^{\prime}}}^{{\prime }} - x_{{O^{\prime}}}^{{\prime }} - \left( {x_{{M^{\prime}}}^{{\prime }} - x_{{O^{\prime}}}^{{\prime }} } \right)\left( {1 - \frac{{R\left( {z_{P} } \right)}}{{R_{0} }}} \right)} \right) $$

(24)

In the same way, for the other horizontal coordinate:

$$ y_{P} - y_{O} = \frac{1}{{R\left( {z_{P} } \right)}}\left( {y_{{P^{\prime}}}^{{\prime }} - y_{{O^{\prime}}}^{{\prime }} - \left( {y_{{M^{\prime}}}^{{\prime }} - y_{{O^{\prime}}}^{{\prime }} } \right)\left( {1 - \frac{{R\left( {z_{P} } \right)}}{{R_{0} }}} \right)} \right) $$

(25)

However, by expressing the terms that depend on $z$ using the resolution $R\left( z \right)$, it is possible to bypass the analysis of the refraction of the optic rays, which may be intricate.

1.2 Calculating the image resolution $R\left( z \right)$ at various distances from the camera

Figure 20 is a scheme of the optical rays, from the camera towards the target. They undergo refraction when crossing surfaces of separation between the two mediums, that is between air and channel wall and between channel wall and water. Note that here, for convenience, the origin of the absolute reference system is different from the one in Fig. 1.

Fig. 20

Refraction of the optic rays

As long as there is no change of medium, the width of the framing, $W$, linearly increases with the distance from the camera. In the area with water, $W$ can be expressed as:

$$ W\left( z \right) \equiv pz + q $$

(26)

The values of p and q are:

$$ p = 2\tan \alpha_{3} $$

(27)

$$ q = 2\left( {d\tan \alpha_{1} + t_{w} \tan \alpha_{2} + d_{0} \tan \alpha_{3} } \right) $$

(28)

where $d$ is the distance between the focal point and the dry side of the channel wall, $t_{w}$ is the thickness of the channel wall, and $d_{0}$ is the distance between the wet side of the channel wall and the reference plane. $\alpha_{1} , \;\alpha_{2} ,\; {\text{and}}\;\alpha_{3}$ are the angles included between the most external optical rays and the vertical direction in air, channel wall, and water respectively.

The angle $\alpha_{1}$ depends only on the camera sensor and lens:

$$ \alpha_{1} = tan^{ - 1} \left( {\frac{w}{{2f_{l} }}} \right) $$

(29)

where $f_{l}$ is the focal length of the lens.

The values of the angles $\alpha_{2}$ and $\alpha_{3}$ can be calculated through Snell’s law:

$$ \alpha_{2} = \sin^{ - 1} \left( {\frac{{n_{1} }}{{n_{2} }}\sin \alpha_{1} } \right) = \sin^{ - 1} \left( {\frac{{n_{1} }}{{n_{2} }}\frac{w}{{\sqrt {4f_{l}^{2} + w^{2} } }}} \right) $$

(30)

$$ \alpha_{3} = \sin^{ - 1} \left( {\frac{{n_{2} }}{{n_{3} }}\sin \alpha_{2} } \right) = \sin^{ - 1} \left( {\frac{{n_{1} }}{{n_{3} }}\sin \alpha_{1} } \right) = \sin^{ - 1} \left( {\frac{{n_{1} }}{{n_{3} }}\frac{w}{{\sqrt {4f_{l}^{2} + w^{2} } }}} \right) $$

(31)

where $n_{1}$, $n_{2}$, and $n_{3}$ are the refractive indexes of the three mediums.

Therefore, $p$ only depends on the characteristics of the camera (and lens) and of the mediums where the camera and the target are. Therefore, it does not depend on the distance of the camera from the channel or on the geometrical characteristics of the channel. On the contrary, the coefficient $q$ depends on the distance between the camera and the channel and on the channel thickness, so its value must be updated every time the position of the camera is changed.

The resolution of the image can be calculated as:

$$ R\left( z \right) = \frac{{w/l_{p} }}{W\left( z \right)} = \frac{{w/l_{p} }}{pz + q} \equiv \frac{{p_{1} }}{{ z + q_{1} }} $$

(32)

where $p_{1} \equiv w/l_{p} /p$ and $q_{1} \equiv q/p$. Again, $p_{1}$ only depends on the characteristics of the camera and the properties of the mediums where the camera and the target are. $q_{1}$ depends also on the specific arrangement of the experiment $(d$, $t_{w}$, $d_{0}$, etc.)

Therefore, by taking pictures at different values, $z_{i}$, and measuring the resolution, $R_{i}$, it is possible to find the coefficient $p_{1}$ through the regression of the points $\left( {z_{i} ,R_{i} } \right)$. $p_{1}$ will be valid for all the experiment arrangements, provided that the camera, the lens and the mediums of the camera, and the target (the sphere) are the same for all the experiments.

Subsequently, only a reference image is necessary to find the coefficient $q_{1}$, which depends on the specific arrangement of the camera and of the channel. Choosing $z = 0$ (the plane where the marker of the sphere at rest lies), the following is straightforward:

$$ R\left( 0 \right) = R_{0} = \frac{{p_{1} }}{{ q_{1} }} $$

(33)

The final expression for the resolution will be:

$$ R\left( z \right) = \frac{{p_{1} }}{{ z + p_{1} /R_{0} }} $$

(34)

1.3 System of equations

The last equation needed to find the real three-dimensional position of the marker $\left( {x_{P} ,y_{P} ,z_{P} } \right)$ is that expressing the three-dimensional surface $F\left( {x_{P} ,y_{P} ,z_{P} } \right) = 0$ on which the marker moves. In the case of the tethered sphere (i.e. pendulum), it is a spherical surface:

$$ F\left( {x_{P} ,y_{P} ,z_{P} } \right) = \left( {x_{P} - x_{O} } \right)^{2} + \left( {y_{P} - y_{O} } \right)^{2} + \left( {\ell - z_{P} } \right)^{2} - \ell^{2} = 0 $$

(35)

The system of equations is thus:

$$ \left\{ {\begin{array}{*{20}c} {x_{P} - x_{O} = \frac{1}{R}\left( {x_{{P^{\prime}}}^{{\prime }} - x_{{O^{\prime}}}^{{\prime }} - \left( {x_{{M^{\prime}}}^{{\prime }} - x_{{O^{\prime}}}^{{\prime }} } \right)\left( {1 - \frac{R}{{R_{0} }}} \right)} \right)} \\ {y_{P} - y_{O} = \frac{1}{R}\left( {y_{{P^{\prime}}}^{{\prime }} - y_{{O^{\prime}}}^{{\prime }} - \left( {y_{{M^{\prime}}}^{{\prime }} - y_{{O^{\prime}}}^{{\prime }} } \right)\left( {1 - \frac{R}{{R_{0} }}} \right)} \right)} \\ {R = R\left( {z_{P} } \right)} \\ {F\left( {x_{P} ,y_{P} ,z_{P} } \right) = 0} \\ \end{array} } \right. $$

(36)

If the functions $R\left( {z_{P} } \right)$ and $F\left( {x_{P} ,y_{P} ,z_{P} } \right)$ are known, the unknowns of the system are $x_{P}$, $y_{P} , z_{P}$, and R and the system is solvable. Therefore, the coordinates of the sphere center can be easily calculated through trigonometry from the coordinates of the marker.