Resonance in vortex-induced in-line vibration at low Reynolds numbers

We present simulations of a circular cylinder undergoing vortex-induced vibration in-line with a free stream in conjunction with a theory for the fluid dynamics. Initially, it is shown that increasing the Reynolds number from 100 to 250 results in a 12-fold increase of the peak response amplitude at a fixed mass ratio of $m^*=5$. Subsequently, we vary $m^*$ from 2 up to 20 at a fixed Reynolds number of 180. The response amplitude as a function of the reduced velocity $U^*$ displays a single excitation region with peak amplitudes of approximately 1\% of the cylinder diameter, irrespectively of the $m^*$ value. The vibration is always excited by the alternating mode of vortex shedding. We develop a new model for the in-line fluid force, which comprises an inviscid inertial force, a quasi-steady drag, and a vortex drag due to the vorticity periodically shed in the wake. Our analysis shows that the vortex drag appropriately captures a gradual shift in the timing of vortex shedding in its phase variation as a function of $U^*$ while its magnitude displays a resonant amplification within the excitation region. We use the theory to illustrate why peak amplitudes, which occur when the vibration frequency is equal to the structural frequency in still fluid, do not depend on $m^*$, in agreement with our simulations as well as previous experimental observations at higher Reynolds numbers than considered here. This new theory provides physical insight which could not be attained heretofore by employing semi-empirical approaches in the literature.


Introduction
The flow periodicity due to the vortex shedding from a bluff body exposed to a fluid stream can cause structural vibration if the body is flexible or it is elastically mounted, which is referred to as 'vortex-induced vibration'. In practical applications, compliant structures usually have degrees of freedom to move both along and across the incident flow. Much of the fundamental research on the problem has dealt with rigid circular cylinders as bluff bodies, constrained elastically so as to have a single degree of freedom to oscillate either in-line with a free stream (the streamwise direction) or transversely (the cross-stream direction). The circular cylinder spawns the characteristic that is not prone to galloping vibration and vortex-induced vibration occurs in its purest form. In an early review on vortex shedding and its applications, King (1977) noted that maximum inline amplitudes are approximately 0.2 diameters, peak-to-peak, or about one-tenth of the corresponding maximum cross-stream amplitudes. As a consequence, subsequent research has mostly concentrated on purely transverse vortex-induced vibration as attested in later reviews on the topic (see, e.g., Bearman 1984;Sarpkaya 2004;Williamson & Govardhan 2004;Gabbai & Benaroya 2005;Bearman 2011;Païdoussis et al. 2011). Insight into the fundamentals of vortex-induced vibration can be gained by the complementary study of either in-line or transverse vortex-induced vibrations since both share the same excitation mechanism. The present study deals with the in-line case.

Characteristics of in-line free response
The condition for the onset of vortex-induced in-line vibration can be broadly expressed as U * a ≈ 1/(2S), whereas the corresponding value for cross-stream vibration is U * a ≈ 1/S, where U * a = U ∞ /f n,a D is the reduced velocity and S = f v0 D/U ∞ is the Strouhal number; here U ∞ is the velocity of the free stream, D is the dimension of the body in the crossstream direction, f v0 is the frequency of vortex shedding for a stationary cylinder, and f n,a is the natural frequency of the structure in still fluid, i.e. including the 'added mass'. The reduced velocity is a non-dimensional parameter, which is typically employed to present experimental data. Table 1 displays various definitions of the reduced velocity, which are used in the present study. Other non-dimensional parameters governing the structural response are: the mass ratio, m * , the damping ratio, ζ, as well as the Reynolds number, Re, which determines the flow regime. The definitions of the mass ratio and the damping ratio are also not uniform in the literature but depend on whether the added fluid mass is taken into account.
Typically, the response amplitude of cylinder vibration is magnified in distinct ranges of the reduced velocity, which mimic the classical resonance of a single degree-of-freedom oscillator to external harmonic forcing. These distinct regions of high-amplitude response have been given various names such as 'instability regions' (King 1977), 'excitation regions' (Naudascher 1987), or 'response branches' (Williamson & Govardhan 2004). It has been established as early as the 1970's that there exist two distinct excitation regions of free in-line vibration: the first one appears at U * a 2.5 and has been associated with symmetrical shedding of vortices simultaneously from both sides of the cylinder, whereas the second one appears at U * a 2.5 and has been associated with alternating shedding of vortices from each side of the cylinder (Wootton et al. 1972;King 1977;Aguirre 1977). The value of U * a ≈ 2.5 corresponds to 1/(2S) assuming a Strouhal number of 0.20, at which reduced velocity the frequency of vortex shedding from a stationary cylinder becomes equal to half the natural frequency of the structure in still fluid, i.e. f v0 ≈ 1 2 f n,a . A factor of 2 arises in the denominator from the fact that two vortices shed from alternate sides of the cylinder each one induces a periodic oscillation of the fluid force in the streamwise direction. It should be remembered that the Strouhal number is a function of the Reynolds number, defined as Re = ρU ∞ D/µ, where ρ is the density and µ is the dynamic viscosity of the fluid. Therefore, the above value of U * a ≈ 2.5 should generally be replaced by U * a ≈ 1/(2S) in dealing with response data at different Reynolds numbers. In experimental studies with elastically mounted rigid cylinders, the structural frequency depends on the oscillating mass and the stiffness of the supporting springs, which provide elastic restoring forces. In an early study, Aguirre (1977) conducted more than a hundred tests in a water channel to investigate vortex-induced in-line vibration. He concluded that the structural mass and the stiffness affected the cylinder response independently and in different ways: for a given value of f /f v0 , the density ratio (i.e. the mass ratio here) did not affect the response amplitude when normalized with the cylinder diameter, nor did the stiffness affect the normalized frequency of vibration f /f n,a , where f is the actual vibration frequency. Beyond that early study, the effects of the structural mass and stiffness, which are embodied in the mass ratio and the reduced velocity values, have not been systematically addressed in more recent studies. Yet, Okajima et al. (2004) found that the response amplitude of in-line oscillation decreases in both excitation regions with increasing the reduced mass-damping, or Scruton number -a nondimensional parameter that is proportional to the product m * ζ and is often employed to compile peak amplitude data as a function of a single parameter. The later findings possibly illustrate the influence of structural damping alone since the mass ratio was constant in those tests. The existence of two distinct response branches of free in-line vibration and corresponding modes of vortex shedding were confirmed in more recent experimental studies at Reynolds numbers in the range approximately from 10 3 to 3.5 × 10 4 (Okajima et al. 2004;Cagney & Balabani 2013b). The drop in response amplitude in-between the two branches, i.e. at U * a ≈ 2.5, has been attributed to the phasing of alternating vortex shedding, which provides a positive-damping, or negative-excitation force with respect to the oscillation of the cylinder (Konstantinidis et al. 2005;Konstantinidis 2014). More recently, a mixed mode of combined symmetric and alternating vortex shedding was also reported to exist in-between the two branches (Gurian et al. 2019).

Energy transfer and harmonic approximation
For self-excited vibrations to be possible, energy must be transferred from the fluid to the structural motion in an average cycle so as to sustain the oscillations, i.e. E 0 where E = F x dx c ; F x is the instantaneous fluid force driving the body motion and x c is the displacement of the body. A rationale is to determine the energy transfer from forced vibrations of the body in order to predict whether free vibrations can occur for the corresponding case where the cylinder is elastically constrained. This is typically based on the approximation that both the motion of the cylinder x c (t) and the driving fluid force per unit length F x (t) can be expressed as single-harmonic functions of time t, e.g.
where X 0 is the mean streamwise displacement of the cylinder, A is the amplitude and f is the frequency of body oscillation, F x0 and F x1 respectively are the magnitudes of the mean and unsteady in-line fluid forces, and φ x is the phase lag between the displacement and the driving force at the oscillation frequency. Using the above harmonic approximations, it can be readily shown that Hence, free vibration is possible in the general case only if sin φ x 0, or equivalently the phase lag be within the range 0 • φ x < 180 • . Two independent studies employing forced harmonic in-line vibrations at fixed amplitudes of oscillation in the range from 0.05 to 0.28 diameters peak-to-peak, have shown that energy is transferred from the fluid to the cylinder motion in two excitation regions separated by approximately U * r ≈ 2.5 for Re values higher than 10 3 (Tanida et al. 1973;Nishihara et al. 2005). It should be noted here that the reduced velocity can only be defined using the actual frequency of forced oscillation (see Table 1). This is also critical in correlating forced-vibration studies with the response from free-vibration studies, in which case the vibration frequency is not necessarily equal to the structural frequency (Williamson & Govardhan 2004;Konstantinidis 2014). Overall, predictions using forced harmonic vibration agree well with the excitation regions found in free vibration at relatively high Reynolds numbers, including the wake modes responsible for free vibration (Tanida et al. 1973;Nishihara et al. 2005). On the contrary, Tanida et al. found that energy transfer was always negative for all reduced velocities at Re = 80. They stated that results obtained for Re = 80 were representative over the range 40 Re 150, which indicated that free vibration may not be possible for relatively low Reynolds numbers.
More recently, detailed results from two-dimensional numerical simulations for a cylinder placed in an oscillating free stream and the equivalent case of a cylinder oscillating in-line with a steady free stream both showed that E < 0 for all reduced velocities at fixed Reynolds numbers of Re = 150 (Konstantinidis & Bouris 2017) and Re = 100 (Kim & Choi 2019). The lowest amplitudes of forced oscillation in those studies were 0.1 and 0.05 diameters, respectively. The findings from both studies have also indicated that free in-line vibration may not be feasible for Reynolds numbers in the laminar regime, which is consistent with the earlier experimental study of Tanida et al. (1973). Nevertheless, it is plausible that energy transfer may become positive at lower amplitudes of oscillation than those employed in previous studies, which has not received attention to date since free in-line vibration has scarcely been studied at low Reynolds numbers; the authors are aware of only few numerical studies where in-line vibration of circular cylinder rotating at prescribed rates was investigated at Re = 100 (Bourguet & Lo Jacono 2015;Lo Jacono et al. 2018). These studies showed that an elastically-mounted rotating cylinder can be excited into large-amplitude galloping-type vibrations as the reduced velocity increases. However, the response amplitudes were negligible in the case of a non-rotating cylinder compared to the rotating cases and the former results were not discussed.
Apart from addressing the question whether in-line vortex-induced vibration is possible at low Reynolds numbers, a more fundamental issue is to clarify what are the flow physics causing variations in the amplitude and frequency of response when self-excited vibration does occur, not only at low Reynolds numbers. This issue impacts our understanding of vortex-induced vibration as well as its modelling and prediction using semi-empirical codes in industrial applications. To address that issue, it is essential to formulate a theoretical framework with the aid of which results can be interpreted. In this study, we maintain that the in-line free vibration offers a convenient test case because it allows different dynamical effects, which are associated with fluid inertia, fluid damping, and the unsteady wake, to be segregated.

Previous theoretical-empirical approaches
A long-standing approach is to represent the in-line force per unit length F x (t) based on the equation proposed by Morison et al. (1950). For a cylinder of circular cross section oscillating in-line with a steady free stream, the equation can be written as where ρ is the density of the fluid. The coefficients C dh and C mh are often referred to as drag and added mass (or inertia) coefficients, respectively, and their values are empirically determined from measurements or simulations. Inherent to this approach is the harmonic approximation since the coefficients C dh and C mh are often determined from tests where the cylinder is forced to vibrate harmonically. Even when the cylinder motion is selfexcited, it is still necessary to characterize the vibration in terms of the least number of appropriate non-dimensional parameters for compiling fluid forcing data; this usually boils down to the use of two parameters, i.e. the normalized amplitude and normalized frequency of oscillation, which can fully characterize only single-harmonic oscillations. Another similar approach is to decompose the fluctuating part of the in-line force into harmonic components in-phase with the displacement (or alternately acceleration) and in-phase with the velocity of the oscillating cylinder, in addition to a steady term for the mean drag. By using harmonic approximations, the steady-state response can be predicted as we present in appendix A. This approach is analogous to using Morison et al.'s equation and linearising the drag term as  Konstantinidis & Bouris (2017) demonstrated the inability of the thus reconstructed force acting on a cylinder in non-zero-mean oscillatory flows to capture fluctuations due to the vortex shedding in the drag-dominated regime, where the vortex shedding and the cylinder motion (the wave motion in that study) are not synchronized. When the primary mode of sub-harmonic synchronization occurs, Morison's equation provides a fairly accurate fit to the in-line force but subtle differences still exist, which may have detrimental effects when the model equation is used for predicting the free response. A disadvantage of previous approaches based on Morison et al.'s equation as well as on the harmonic representation, is that the values of the force coefficients show some dependency on the best-fitting method, e.g. Fourier averaging vs. least-squares method (see Konstantinidis & Bouris 2017). Another disadvantage, more important, is that force coefficients are empirically determined and as a consequence it is difficult to decipher the flow physics from the variation of the force coefficients, which is our primary goal in this work.

Force decomposition and added mass
Despite the empirical use of Morison et al.'s equation, its inventors stated that it originates from the summation of a quasi-steady drag force and the added-mass force resulting from 'wave theory' (Morison et al. 1950). A comprehensive discussion of the theory can be found in Lighthill (1986). For a body accelerating rectilinearly within a fluid medium, there is an ideal 'potential' force acting on the body, which can be expressed as F x,potential = −C a m dẍc , where C a is the added mass coefficient, m d is the mass of fluid displaced by the body, andẍ c is the acceleration of the body. Thus, the body behaves as if it has a total mass of m + m a , where m a = C a m d is the added mass of fluid. According to the theory of inviscid flow, in which case the velocity field can be defined by the flow potential, the added mass coefficient C a of any body is assumed to depend exclusively on the shape of the body. For a circular cylinder, C a = 1. Free-decay oscillation tests in quiescent fluid have shown that C a is quite close to the ideal value of unity. However, the applicability of the ideal C a value in general flows, including cylinders oscillating normal to a free stream, has been criticized (Sarpkaya 1979(Sarpkaya , 2001(Sarpkaya , 2004. On the other hand, Khalak & Williamson (1996) argued that removing the ideal added-mass force from the total force will leave a viscous force that may still comprise a component in-phase with acceleration, i.e. the decomposition does not have to separate all of the accelerationdepended forces as done in empirical approaches. Ever since the separation of 'potential' (inviscid) and 'vortex' (viscous) components has been widely employed to shed light into the vortex dynamics around oscillating bodies transversely to a free stream (see, e.g., Govardhan & Williamson 2000;Carberry et al. 2005;Morse & Williamson 2009;Zhao et al. 2014Zhao et al. , 2018Soti et al. 2018).
For body oscillations in-line with a free stream, one may also split the streamwise force as (1.5) in order to explore the link between fluid forcing and vortex dynamics. This was previously done for the case of a fixed cylinder placed normal to a free stream with smallamplitude sinusoidal oscillations superimposed on a mean velocity (Konstantinidis & Liang 2011). This case is kinematically equivalent to the forced vibration of the cylinder in-line with a steady free stream. In that study, large-eddy simulations corresponding to Re = 2150 showed that alternating vortex shedding provides positive energy transfer for U * r > 2.5, in very good agreement with previous experimental studies discussed earlier. It was also observed that the streamwise vortex force diminished in magnitude while the instantaneous phase of the vortex force with respect to the imposed oscillation drifted continuously near the middle of the wake resonance (synchronization) region. This was considered to be inconsistent with the flow physics in the following sense: within the synchronization region the vortex shedding and the oscillation are strongly phaselocked and the corresponding wake fluctuations are resonantly intensified. Therefore, the magnitude of the vortex force would have been expected to increase and its instantaneous phase to remain fairly constant in this region. The irregular phase dynamics observed in that study indicates that the vortex force remaining from subtracting the ideal inertial force from the total force may not fully represent the effect of the unsteady vortex motions on the fluid forcing.

New theory and outline of the present approach
In this paper, we develop a new theoretical model for representing the streamwise force on a cylinder oscillating in-line with a free stream. The model stems from some recent observations. In particular, it was recently shown that Morison et al.'s equation based on the sum of a quasi-steady viscous drag force and an inviscid inertial force represents the in-line force with comparable accuracy as does the equation with best-fitted coefficients over a wide range of parameters from the inertia to drag-dominated regimes (Konstantinidis & Bouris 2017). However, neither method could capture fluctuations at the vortex shedding frequency in the drag-dominated regime, as noted earlier. Thus, the idea here is to introduce an independent term for the 'vortex drag', i.e. to express the total force as where the first term represents the quasi-steady drag, the second term represents the inviscid added-mass force, and the third term represents the force due to the vorticity shed in the wake. In this new approach, there are two viscous contributions: a quasi-steady drag, which is an 'instantaneous' reaction force, and a vortex drag, which represents the 'memory' effect in a time-dependent flow. These contributions may be thought of as related to the vorticity bounded around the cylinder and that shed in the wake, respectively. Both are affected by the rate of diffusion of the vorticity, which is finite. However, at sufficiently high Reynolds numbers for which separation occurs, the region of bounded vorticity is thin and the diffusion in this region occurs fast enough, almost 'instantaneously'. Then, the force required to supply the rate of increase of the kinetic energy of the rotational motion in the region of bounded vorticity may be taken to be proportional to the square of the relative velocity U ∞ −ẋ c , which gives rise to a quasi-steady drag (Lighthill 1986). On the other hand, the diffusion of vorticity in the separated wake region is a very complex process (see, e.g., Konstantinidis & Bouris 2016); as a consequence the resulting vortex force acting on the body depends on the history of the vortex motions.
Assuming that a single periodic mode of vortex shedding occurs, the vortex force can be modelled as a single-harmonic function of time, i.e.
where the coefficient C dv represents the magnitude of the unsteady vortex force, and φ dv the phase between the vortex force and the displacement of the cylinder. The frequency f dv excited by the unsteady wake depends on the vortex shedding mode, e.g. f dv = 2f vs for the alternating mode, whereas f dv = f vs for the symmetrical mode. The variation of C dv and φ dv as functions of the normalized amplitude and normalized frequency of oscillation, as well the Reynolds number, is expected to be instructive of the flow dynamics in the wake of vibrating cylinders.
For this study, we conducted numerical simulations of the vortex-induced vibration of a circular cylinder elastically constrained so as to oscillate only in-line with a free stream. Simulations were restricted to the two-dimensional regime to keep computer time within reason so as to examine the influence of the reduced velocity and the mass ratio over wide ranges and with a good resolution. The two-dimensional approach prohibits the realistic simulation of cylinder-wake flows at Reynolds numbers above the threshold for which three-dimensional instabilities occur, which may be assumed approximately the same as for a stationary cylinder, i.e. Re c ≈ 190 (Williamson 1996). However, we have conducted a number of simulations at slightly higher Reynolds numbers up to 250 to check the trends in the results. The numerical simulations provided complete sets of data for the flow fields and induced fluid forces, and their interaction with the resulting free motion of the cylinder, which allowed us to address the issues raised in the foregoing paragraphs and hopefully make a contribution to the understanding of the complex phenomena encountered in vortex-induced vibration.

Governing equations
The flow is assumed incompressible and two-dimensional while physical properties of the fluid are constant. The fluid motion is governed by the momentum (Navier-Stokes) and continuity equations, which can be written in non-dimensional form using the pressure-velocity formulation as Coordinates have been normalized with D, fluid velocities have been normalized with U ∞ , time with D/U ∞ , and of pressure with 0.5ρU 2 ∞ . The acceleration of the cylinderẍ * c appears on the right-hand side of equation (2.1) because the Navier-Stokes equations are applied in a non-inertial frame of reference that moves with the vibrating cylinder. Instead of explicitly enforcing the continuity equation, the the pressure field was computed by solving the following Poisson equation at each time step, where D = ∂u x /∂x + ∂u y /∂y is the dilation. Although the dilation is zero by default in incompressible flows (Eq. 2.3), the term ∂D/∂t is kept in Eq. (2.4) to avoid the propagation of numerical inaccuracies (Harlow & Welch 1965). On the cylinder surface, the no-slip boundary condition gives and the following condition for the normal pressure gradient at the wall where n refers to the component normal to the cylinder surface pointing to the fluid side. At the far field, a potential flow field is assumed so that where u x,pot and u y,pot is the velocity field from the known potential of irrotational flow. The corresponding condition for the far-field pressure is The initial field corresponds to the potential flow around a circular cylinder. The dimensionless form of the equations illustrates that the fluid motion depends solely on the Reynolds number Re given the boundary and initial conditions. The displacement of a cylinder elastically constrained so that it can oscillate only inline with a uniform free stream is governed by Newton's law of motion, which can be written in non-dimensional form as where x * c ,ẋ * c , andẍ * c respectively are the non-dimensional displacement, velocity and acceleration of the cylinder, normalized using D and U ∞ as length and velocity scales; F * x (t) is the sectional fluid force on the cylinder normalized by 0.5ρU 2 ∞ D; U * = U ∞ /(f n D) is the reduced velocity based on the natural frequency of the system in vacuum; m * is the ratio of the cylinder mass to the fluid mass displaced by the cylinder; ζ is the ratio of the structural damping to the critical damping at which the mechanical system can exhibit oscillatory response to external forcing (the dimensional form of the equation of motion and the definition of non-dimensional parameters are given in appendix A). Here, the forcing is provided by the ambient flow through balancing the normal and shear stresses on the cylinder surface. The cylinder is initially at rest. The dimensionless form of equation (2.9) illustrates that the cylinder motion depends on the reduced velocity U * , the mass ratio m * , and the damping ratio ζ. The full set of independent dimensionless parameters of the problem comprises Re, U * , m * , and ζ.

Numerical code
An in-house code based on the finite difference method was used to solve the equations of fluid motion (Baranyi 2008). The flow domain is enclosed between two concentric circles: the inner circle is the boundary fitted to the cylinder surface while the outer circle represents the far field boundary. The polar physical domain is mapped into a rectangular computational domain using linear mapping functions. The computational mesh of the 'physical domain' is fine in the vicinity of the cylinder and coarse in the far field while the corresponding mesh of the transformed domain is equidistant. Space derivatives are approximated using a fourth order finite-difference scheme except for the convective terms for which a third-order modified upwind difference scheme is employed. The pressure Poisson equation is solved using the successive over-relaxation (SOR) method and the continuity equation is implicitly satisfied at each time step. The Navier-Stokes equations are integrated explicitly using the first-order Euler method and the fourth-order Runge-Kutta scheme is employed to integrate the equation of cylinder motion in time. At each time step the fluid forces acting on the cylinder are calculated by integrating the pressure and shear stress at the cylinder surface, which are obtained from the flow solver. The streamwise force is supplied to the right-hand side of Eq. (2.9) which is integrated to advance the cylinder motion. At the next time step, the cylinder acceleration is updated and the equations of fluid motion are integrated to complete the fluid-solid coupling. For all simulations reported here, the cylinder is initially at rest and the initial field around the cylinder satisfies the potential flow.

Domain size, grid resolution and time-step dependence studies
In this section, we present results from preliminary simulations to check the the dependence of main output parameters on a) the size of the computational domain in terms of the radius ratio R 2 /R 1 , b) the grid resolution ξ max × η max , and c) the dimensionless time step ∆t. Here ξ max and η max are the number of grid points in peripheral and radial directions, respectively. During these computations, the Reynolds number, the reduced velocity, the mass ratio, and the structural damping ratio were fixed at Re = 180, U * = 2.55, m * = 10, and, ζ = 0 respectively. The main output parameters of interest are the normalized amplitude A * and the normalized frequency f * of cylinder response as well as the in-line and transverse force coefficients, C x and C y , respectively.
First, three different values of the radius ratio R 2 /R 1 of the inner and outer circles defining the computational domain were tested and the corresponding results are shown in Table 2 in each case in order to keep the grid equidistant in the transformed domain. For these computations, the time step was fixed at ∆t = 10 −4 U * ∼ = 0.0002. Table 2 shows that the most sensitive quantity on the domain size is C x displaying a relative difference of 1.7% between the small and large domains, whereas the corresponding differences for A * , f * and C y are below 0.7%. The relative differences for all quantities of interest become less than 0.6% between the medium and large domains. Thus, the medium-sized domain with a radius ratio of R 2 /R 1 = 160 was chosen for the rest of the computations. Next, we tested three grids with different resolutions ξ max × η max where the number of peripheral and radial grid points was increased so that the grid remains equidistant in the transformed plane. For these computations, the radius ratio and dimensionless time step values were fixed at R 2 /R 1 = 160 and ∆t = 10 −4 U * ∼ = 0.0002, respectively. The results of these tests are shown in Table 3. Again, C x is the most sensitive quantity displaying a relative difference of 0.7% between coarse and fine grids. All quantities of interest display relative differences less than 0.3% between medium and fine grids. Thus, the medium-resolution grid was chosen for further computations.
Finally, we tested the dependence on the dimensionless time step, ∆t for three different values corresponding to ∆t ∼ = 2 × 10 −4 U * , 10 −4 U * and 5 × 10 −5 U * . For these computations, the radius ratio and grid resolution values were fixed at R 2 /R 1 = 160 and 360 × 292, respectively. The results are shown in Table 4. In this tests, the quantity that is most sensitive to the time step is A * , which shows a relative difference of 0.4% between the largest and smallest time steps whereas the corresponding relative differences for f * , C x and C y are all below 0.15%. The relative differences of all quantities of interest are below 0.1% between the intermediate and the smallest time steps. Thus, a dimensionless time step of ∆t = 10 −4 U * was chosen for the main computations.

Code validation
The flow code used in the present study was previously employed in several studies of flows about stationary and oscillating cylinders and results have been extensively compared against data from the literature (see Baranyi 2008;Dorogi & Baranyi 2018. For instance, Baranyi (2008) has found good agreement with the study of Al-Mdallal et al. (2007) in terms of the time history of the lift coefficient and Lissajous patterns of lift vs. cylinder displacement at comparable situations for the case of a cylinder forced to oscillate in the streamwise direction. In addition, the extended code handling flow-structure interaction has been validated for the case of a cylinder undergoing vortexinduced vibration with two degrees of freedom with equal natural frequencies in the streamwise and transverse directions against results from Prasanth & Mittal (2008) for m * = 10, ζ = 0, and Reynolds numbers in the range from 60 to 240 (for details see Dorogi & Baranyi 2018). Furthermore, Dorogi & Baranyi (2019) showed that results obtained with the present code compare well with published results in Navrose & Mittal (2017) for purely transverse free vibration, as well as in (Prasanth et al. 2011) andBao et al. (2012) for free vibration with two degrees of freedom with equal or unequal, respectively, natural frequencies in the streamwise and transverse directions, for similar conditions in each case.
In addition to the validation tests presented in previous studies, we compare in figure 1 results obtained with the present code against the study of Bourguet & Lo Jacono (2015) for purely in-line free vibration. Here, we employed a finer step in the reduced velocity to capture the peak amplitude as well as the minimum of the unsteady in-line force. There is excellent agreement of results for the response amplitude A * and the fluctuating in-line force coefficient C x but there are some minor deviations regarding the response frequency f * and the transverse force coefficient C y of 0.2% and 1.2%, respectively, which might be attributable to different numerical methods employed in those studies (finite difference vs. spectral element). Overall, previous and present validation tests show that the numerical code employed in the present study provides accurate solutions.

Results and discussion
3.1. Effect of Reynolds number on in-line response at a fixed mass ratio To start with, we present in figure 2 the significant effect of the Reynolds number on the in-line response of a cylinder with a mass ratio of m * = 5. The structural damping was set to zero so as to allow for the highest possible amplitude response to take place. For all Reynolds numbers there exists a single excitation region in which the response amplitude A * displays a marked peak. The peak amplitude over the entire U * range, denoted as A * max , increases from 0.002 at Re = 100 to 0.024 at Re = 250, i.e. an increase in Re by a factor of 2.5 results in an increase in A * max by a factor of 12. This is a remarkable increase of the order of magnitude, which contrasts the constancy of peak amplitudes of purely transverse free vibration in the corresponding range of Reynolds numbers (a compilation of A * max data as a function of Re from several studies can be found in Govardhan & Williamson 2006). The U * value at which peak amplitudes occur decreases with Re, which can be attributable to the corresponding increase of the Strouhal number since peak amplitudes occur at U * a ≈ 1/(2S). The normalized frequency f * is approximately twice the corresponding Strouhal number at each Reynolds number. However, f * displays a trough within the excitation region that becomes more pronounced as the Reynolds number increases; it can be hardly discerned for Re = 100. The variations of A * and f * appear to be strongly correlated. The characteristics of free response remain similar for all Reynolds numbers considered here. However, a sudden drop in A * appears just after the peak amplitude for Re = 250, a feature which is not present for lower Reynolds numbers. This might indicate the onset of branching behaviour similar to that observed in vortex-induced vibration purely transverse to the free stream (see Leontini et al. 2006). In comparison to previous experimental studies, which typically correspond to Reynolds numbers above 10 3 , we did not observe another excitation region associated with symmetrical vortex shedding. In contrast, we observed only the alternating mode of vortex shedding in the present simulations corresponding to the laminar wake regime. Figure 3 shows vorticity distributions in the wake at U * values corresponding to peak amplitudes for different Reynolds numbers. In all cases, the vorticity distributions display the familiar von Kármán vortex street similar to the wake of a stationary cylinder. As the Reynolds number is increased, the contours of individual vortices become more concentrated and peak vorticity values within them increase. This might be partly attributable to the increase in the amplitude of cylinder oscillation with Reynolds number, which accrues the generation of vorticity on the cylinder surface (Konstantinidis & Bouris 2016). In addition, the streamwise spacing between the centres of subsequent vortices decreases due to the increase of the normalized frequency of cylinder oscillation f * with Reynolds number. The absence of the other excitation region associated with the symmetrical vortex shedding may be attributable to the fact that, as has been shown in several previous studies where the cylinder is forced to oscillate in the streamwise direction at correspondingly low Reynolds numbers, the onset of this mode occurs at relatively high amplitudes above 0.1 diameters (Al-Mdallal et al.

2007; Marzouk & Nayfeh 2009; Kim & Choi 2019
). Since streamwise amplitudes of free vibration are much lower than that threshold, it is not surprising that the mode of symmetrical shedding and the corresponding excitation region were not observed in the present study.

Effect of mass ratio on in-line response at a fixed Reynolds number
Next, we concentrate on the effect of mass ratio on the in-line response at a fixed Reynolds number of Re = 180, at which the flow is expected to remain laminar and strictly two dimensional. The structural damping was set to zero (ζ = 0) to allow the highest possible amplitudes.  Figure 4 shows the variations of A * and f * with U * for four m * values. It can be seen that A * displays a single excitation region with peak amplitudes of approximately 1% of the cylinder diameter, irrespectively of m * . At high reduced velocities, i.e. U * > 4, the response amplitude gradually drops off down to a level that depends on the mass ratio, with A * becoming lower as m * increases. The response frequency f * initially decreases within the excitation region reaching a minimum value of approximately 0.372 for all mass ratios. As U * is increased beyond the point of minimum, f * increases asymptotically towards the value corresponding to twice the Strouhal number for a stationary cylinder, i.e. f * ≈ 2S = 0.384. It is interesting to note that f * < 2S over the entire U * range for all m * , i.e. the cylinder always oscillates at a frequency slightly lower than twice the frequency of vortex shedding from a stationary cylinder. This is consistent with forced harmonic vibration studies, which show that vortex-induced vibration due to alternating vortex shedding occurs for f < 2f v0 (Nishihara et al. 2005;Konstantinidis & Liang 2011).

Cylinder response
The variations of A * and f * shown in figure 4 appear to be strongly correlated, i.e. A * increases as f * decreases and vice versa. However, it should be noted that peak amplitudes occur at marginally lower U * values than those that correspond to the minimum in f * . The U * values at which peak amplitudes occur depend on m * quite substantially. In fact, peak amplitudes occur approximately at the point where the frequency of cylinder oscillation approaches the natural frequency of the system in still fluid, i.e. f ≈ f n,a . This condition can also be expressed in dimensionless form as: Normalized frequency at point of peak amplitude : The above condition can be verified in Table 5, which summarizes the response charac- teristics at peak amplitudes for different mass ratios. Effectively, we see that the reduced velocity at which the peak amplitude occurs primarily depends on the mass ratio but the peak response amplitude does not depend on the mass ratio. The drop in normalized frequency f * within the excitation region in fact illustrates the tendency of the oscillation to 'lock-in' at the natural frequency of the structure in still fluid, i.e. f ≈ f n,a , over a range of reduced velocities. It is important to distinguish this 'lock-in' tendency, which only occurs in free vibration, from 'vortex lock-in' (a.k.a. 'vortex lock-on'), which regards the synchronization of the vortex shedding and the cylinder oscillation and may occur in both forced and free vibration (Konstantinidis 2014). In forced vibration, there is no natural frequency of the structure so the lock-in relationship f = f n,a is meaningless. It should also be noted that for all simulations reported in figure 4 (i.e. for all U * and m * values), the vortex shedding was synchronized with the cylinder oscillation so that f = 2f vs , where f vs is the frequency of vortex shedding from the freely vibrating cylinder. This is tantamount to the sub-harmonic vortex lock-on in the context of forced oscillations where alternating vortex shedding synchronizes at half the frequency of cylinder oscillation, i.e. f vs = 1 2 f (see, e.g., Kim & Choi 2019). However, the conventional lock-in f ≈ f n,a occurs over a narrow range of reduced velocities in the region of peak-amplitude response.

Magnitude and phase of fluid forces
In figure 5 we present the variation of the unsteady force coefficients in-line with and transversely to the free stream, respectively denoted as C x and C y , as functions of U * . The dashed lines indicate constant values corresponding to a stationary cylinder at Re = 180 (taken from Qu et al. 2013). C x exhibits similar trends with U * for all mass ratios. Initially, C x increases gradually reaching a peak at approximately the point of peak amplitude response. A maximum value of C x = 0.106 is attained for all mass ratios, which is nearly three times the value corresponding to a stationary cylinder; a substantial increase despite the low amplitude of peak oscillation of only 1% of the cylinder diameter. Following the peak there is a steep decrease within a narrow U * range at the end of which C x tends to zero. Interestingly, for all mass ratios this occurs at U * = 2.625, a point at which the response frequency coincides with the natural frequency of the structure in vacuum, i.e. f = f n (or in non-dimensional parameters f * = 1/U * ). Hereafter, this special operating point will be referred to as the 'coincidence point' and its ramifications will be discussed in more detail in Sect. 3.3. Beyond the coincidence point, C x gradually reaches to a plateau at a value that is proportional to m * .
It can be seen in the bottom plot in figure 5 that the variation of C y as a function of U * is very similar to that of C x (cf. top plot). However, peak and trough C y values are only 6% higher and 5% lower, respectively, than the value corresponding to a stationary cylinder. For all m * , a maximum value of C y = 0.631 exactly is attained at the point of peak response amplitude. In addition, at the coincidence point, i.e. U * = 2.625, C y attains the same value of 0.591 exactly for all m * , which is just slightly below the value corresponding to a stationary cylinder. Since there is no body acceleration in the transverse direction, changes of C y can be related to changes in the vortex dynamics around the oscillating cylinder (Leontini et al. 2013). Then, the small variation in the C y magnitude illustrates that the process of vortex formation and shedding does not substantially change over the entire range of reduced velocities, even though the cylinder can be oscillating with different but generally small amplitudes. The phase angle φ x between the driving force F x (t) and cylinder displacement x c (t), as defined by harmonic approximations in equations (1.1) and (1.2), is often considered to be a useful parameter in studies of vortex-induced vibration. In addition to φ x , we also compute here the phase angle φ y between the unsteady force acting in the transverse direction F y (t) and the displacement, assuming that F y (t) can also be approximated as a single-harmonic function of time. The calculation method of the phase angle is described in appendix B. It should be noted that the harmonic approximations work very well at low Reynolds numbers considered here except for the time history of F x (t) near the coincidence point as will be discussed in detail further below. Figure 6 shows the variations of the phase angles φ x and φ y as functions of U * .
In figure 6, it can be seen that φ x jumps suddenly from 0 • to 180 • across the coincidence point at U * = 2.625 for all m * . This behaviour can be predicted from the steady-state harmonic solution, which is given in appendix A, as follows. Equation (A 6) requires that the condition sin φ x = 0 be satisfied when the structural damping is null (ζ = 0). This constrains φ x to be either 0 • or 180 • . This taken in tandem with the requirement per equation (A 7) that cos φ x must change from positive to negative as U * increases through f * U * = 1 translates to the jump in φ x from 0 • to 180 • appearing at the coincidence point. We would like to stress that the equation of cylinder motion constraints φ x , which makes it impossible to infer any changes in the flow from the variation of φ x as the reduced velocity is varied. In contrast to φ x , there is no analogous constraint on φ y , which instead of a jump displays a smooth variation with U * , as can be seen in figure 6. For all mass ratios, φ y increases from an initial value of approximately 20 • at low reduced velocities to a terminal value of approximately 105 • at high reduced velocities, with precise values being slightly depended on m * . The U * range over which φ y changes rapidly is consistent with the range of relatively high-amplitude response, which broadens as m * decreases (see figure 4). The variation of φ y clearly suggests a gradual change in the vortex dynamics as U * is varied over the prescribed range. Based on evidence from previous studies (see Konstantinidis et al. 2005;Konstantinidis & Liang 2011), our main hypothesis is that the variation of the phase angle φ y can be linked to a gradual shift in the timing of vortex shedding with respect to the cylinder oscillation as U * is increased. This hypothesis will be verified in the following subsection by inspection of vorticity distributions at different phases of the cylinder oscillation.

Vorticity distributions
We have selected three U * values, which are listed in table 6, for presenting vorticity distributions in the wake. In purpose, the normalized frequency of cylinder oscillation f * is nearly the same in all three cases so that the vorticity patterns are easier to describe since the streamwise spacing of the shed vortices scales with f * (Griffin 1978). However, the frequency ratio f /f n = U * f * and the phase angle φ y both increase with U * , at different degrees, as also shown in table 6. The two higher U * values are just before and after the 'coincidence point', U * f * = 1. For each U * , figure 7 shows instantaneous vorticity distributions at two phases corresponding to the maximum and the subsequent zero displacement of the cylinder as indicated in the time traces (see bottom plots in   figure 7). The significant change of the phase angle φ y as a function of U * can be readily inferred from the time traces of the displacement and the transverse force, which have been appropriately normalized with their maximum amplitudes for better visualisation of their waveforms. As U * is increased, the position of individual vortices at the same instants appears to have shifted slightly downstream. The shift in the streamwise position of individual vortices fits very well with the change in φ y as U * is varied; when U * changes from 2.35 to 2.6, φ y more than doubles and a large shift is observed; when U * changes from 2.6 to 2.7, φ y changes by few degrees and the shift is hardly perceptible. Making an estimation of the phase change from the position of the vortex centres when U * changes from 2.35 to 2.6 yields a shift of 54 • , which is very close to the computed change of 55.4 • in φ y . Similar inferences can be made by looking into the stage of formation of vortices just behind the cylinder. For instance, a stripe of negative vorticity (in blue colour) connecting the second vortex to the third vortex taken from right to left, i.e. as time progresses, can be observed at maximum displacement for U * = 2.35. However, the corresponding vortices are no longer connected by a stripe of vorticity for U * = 2.6 and 2.7, which illustrates that this instant corresponds to a later stage in the vortex shedding process. These observations provide firm evidence that the variation of the phase angle φ y as a function of U * is directly linked to changes in the timing of vortex shedding from the cylinder.

Variable added mass
Changes in the frequency of cylinder response are often correlated with variations in the inertial force due to the added mass . This follows from the relationship Here, C a is the ideal added mass coefficient from potential flow theory and C EA is a variable added mass coefficient (Aguirre 1977). The latter coefficient has become known as the 'effective added mass' in the context of transverse free vibration (Khalak & Williamson 1996;Williamson & Govardhan 2004). Equation (3.2) is equivalent to (A 7) of the harmonic approximation solution where the component of the force in-phase with displacement has been substituted by and the ratio of the natural structural frequency in vacuum to that in still fluid is given by In figure 8, we present the variation of C EA as a function of U * for each m * investigated. C EA values were computed via equation (3.2) and the known frequency ratio; note that f /f n,a = f * U * a = f * U * (f n /f n,a ). It can be seen that C EA decreases continuously from positive to negative values with U * ; C EA decreases from as high as 40.8 to as low as −13.3 for m * = 20. At some point, C EA takes the theoretical C a value of unity, which is indicated by the dashed line in the inset of figure 8. Interestingly, this occurs at approximately the reduced velocity of peak amplitude response, which is given by equation (3.1). Furthermore, for m * = 20 a gap separating operating points with C EA values above and below unity appears at U * = 2.614, which is not present at lower m * values. The discontinuous variation remained in place even though an extremely fine step of ∆U * = 0.002 was employed around this point. When the frequency of cylinder oscillation approaches the natural frequency of the structure in vacuum, equation (3.2) yields C EA = 0, which occurs at the coincidence point U * = 2.625 for all mass ratios.
The very wide variation of C EA values as a function of U * is improbable to represent some physical change, e.g. due to inertial effects, as we have already seen that flow physics remain fairly robust over the entire range of reduced velocities. We point out this simply because we would like to decipher the true effect due to the added mass, which seems impossible to do through the effective added mass concept. On the other hand, it is shown further below with the aid of the new theory that a constant value of the added mass coefficient aligns very well the observations from the simulations. Nonetheless, the empirical values of C EA as defined above are valid.

Super-harmonic fluid forcing at the coincidence point
As already pointed out, the magnitude of the unsteady in-line force coefficient tends to zero level as the frequency of cylinder vibration approaches the natural frequency of the structure in vacuum. This occurs at U * ≈ 2.625 for all mass ratios. This can be predicted from the harmonic solution as follows. Equations (A 6) and (A 7) given in appendix A can be combined so as to eliminate the phase angle φ x . In the case of zero structural damping (ζ = 0), the following relationship is obtained (3.5) The above relationship shows that if the vibration frequency approaches the natural frequency of the structure in vacuum, which can be expressed in non-dimensional form as f * U * → 1, then the only feasible solution is that C x1 → 0. In the present simulations we have found that for all mass ratios f approaches f n within less than 0.1%, which was made possible by using very fine steps of ∆U * = 0.01 or even less. Figure 9 shows time series of the cylinder displacement and the fluctuating part of the in-line force at different reduced velocities. The displacement and the force were normalized with their corresponding maxima in order to reveal their relative waveforms because absolute values of the force are extremely small. Although the displacement remains remarkably harmonic, the force starts to deviate from pure harmonic as U * is varied with a fine step around the coincidence point. Spectra of the force show that a secondary peak appears at the first super-harmonic of the vibration frequency. At U * = 2.63, the main spectral peak appears at the first super-harmonic of the vibration frequency, whereas a very small peak remains at the main frequency of vibration; the latter corresponds to an extremely small magnitude of C x1 , which is however sufficient to sustain the vibration at the main harmonic of the vibration.

Development of new linear theory
In this section, we develop a theoretical framework for fluid mechanics of vortexinduced in-line vibration. We begin with linearising the quasi-steady drag term in equation (1.6) and substituting the harmonic approximations in equations (1.1) and (1.7) into the linearised force, which we express explicitly as a function of time t as where ω = 2πf . The first term on the right-hand side represents the steady part of the in-line force whereas the remaining terms represent the unsteady part, which comprises the combination of cosine and sine functions at the frequency of cylinder oscillation. Equating the steady parts in equations (1.2) and (4.1) yields From equation (4.2) follows that C d is the mean drag coefficient. Equating the cosine and sine terms in equations (1.2) and (4.1) yields the following relationships: The set of equations (4.2-4.4) establishes relationships between fluid forcing components as functions of A * and f * alone, i.e. the structural parameters do not appear. Therefore, the same values of the fluid forcing components also hold when the cylinder is forced to oscillate at the same operating points in the A * : f * parameter space as the free vibration. We use the above set of relationships in order to determine C dv and φ dv appearing in the new theoretical model from the cylinder response (A * , f * ) and the fluid forcing (C x1 , φ x ) obtained from the simulations. It should be remembered here that C a = 1 is the known ideal added-mass coefficient of unity. U $ a f $ ? dv Figure 10. The variation of the vortex drag magnitude C dv and phase φ dv as functions of U * a f * for different mass ratios (symbol legend as in figure 5). Figure 10 shows the variation of C dv and φ dv as functions of the parameter U * a f * . We have employed this parameter because peak amplitudes occur at U * a f * ≈ 1. As can be seen in the top plot, C dv also displays a peak at the same point with a maximum value of 0.067 for all m * . The mutual amplification of the response amplitude and the vortex drag are representative of resonance. Hence, we refer to the condition U * a f * = 1 as the 'resonance point'. Away from resonance C dv tends asymptotically to 0.036, which corresponds to the fluctuating drag coefficient for a non-vibrating cylinder at Re = 180.
In the bottom plot in figure 10 can be seen that φ dv follows an S-type increase as a function of U * a f * from approximately 0 • to 180 • . At the resonance point, U * a f * = 1, φ dv passes through 90 • for all m * . Overall, the variation of φ dv is very similar to that of φ y ; in fact there is a direct relationship between them, which is illustrated in figure 11. Since the variation of φ y is invariably linked to the vortex dynamics, in particular to the timing of vortex shedding as shown earlier, it follows that φ dv appropriately captures related changes as the reduced velocity is varied.

Peak response at resonance point
At the resonance point, we can use the fact that φ dv = 90 • irrespectively of the m * value. In this case, equations (4.3) and (4.4) can be simplified to obtain (4.6) The above set of equations establish relationships between the vortex, mean, and total drag coefficients at resonance point, which involve A * and f * . Considering that all three force coefficients, C dv , C d , and C x1 , can be 'uniquely' specified in the parameter space  Figure 11. Relationship between phase angles φ dv and φy for different mass ratios (symbol legend as in figure 5).
resonance point, equation (4.5) illustrates that the component of the vortex force inphase with the cylinder velocity counterbalances the quasi-steady drag force, whereas equation (4.6) illustrates that the only contribution of the fluid force in-phase with the cylinder acceleration is the inviscid added-mass force. This is commensurate to the situation where an elastically mounted cylinder oscillates freely within a fluid medium with no net viscous forces; in such a situation, it would be anticipated that the frequency of cylinder oscillation is equal to the natural frequency of the structure in an inviscid fluid, i.e. f = f n,a . This is as if the fluid was phenomenologically interacting with the cylinder motion only through its ideal inertia. It should be emphasized that in the case of a viscous fluid, the zero net viscous force results from the cancelling out of the vortex and quasi-steady drag forces, which is brought about by a gradual change in the phasing of the vortex shedding, thereby in φ y and φ dv , as the the reduced velocity is varied. The above changes can be viewed from another perspective, which is more insightful. As U * increases, the phasing of vortex shedding gradually shifts to follow the changes in the frequency of oscillation with the flow velocity. Since these variations occur in a continuous manner for the low-Re cases considered here, a point is reached where the timing of vortex shedding induces a vortex drag exactly in-phase with the cylinder velocity, φ dv = 90 • . At this point, the net viscous force must become null (Eq. 4.5), whereas only the added mass contributes to the force in-phase with the cylinder acceleration (Eq. 4.6). It is important to note again that these changes are brought in entirely by the fluid dynamics.
We can also combine the new hydrodynamic model with the equation of cylinder motion in order to illustrate the interaction between the fluid and the structural dynamics. The steady-state solution can be obtained following a similar procedure as described in appendix A, which results in the following relationships: (4.8) When the structural damping is zero (ζ = 0), equation (4.7) becomes identical to (4.5), and its solution apparently does not depend on m * . It should be noted however that the feasible operating points in the A * : f * space are limited to those that satisfy equation (4.7), which correspond to the contour of zero energy transfer (sin φ x = 0). In addition, equation ( This result shows that the steady-state solution does not depend on m * ; rather the operating point is solely determined by the fluid dynamics as explained above. When the structural damping is finite positive (ζ > 0), equation (4.9) may be used to estimate approximately the peak amplitude for very low values of the mass and the damping such that π 2 m * ζ/U * C d , on the condition that the resonance point is attained. The early experiments of Aguirre (1977) support the above inference; he observed little variation of the peak amplitude in the second excitation region, i.e. the one associated with alternating vortex shedding, which occurred at U * a ≈ 1.2/(2S) for m * = 1.2 and 4.3 (see figure 43 of his thesis). The mass-damping was reported in terms of a stability parameter k s0 to be approximately 0.5, which was estimated from free-decay tests in still water, i.e. it includes contributions from hydrodynamic damping from the cylinder as well as its mountings. From his data, we estimate here that the corresponding k s value in air was significantly lower than the one measured in water so that π 2 m * ζ/U * < 0.05; this is small compared to the mean drag coefficient C d . Thus, we argue that the present analysis is consistent with the experimental facts at Reynolds numbers higher than considered in this study.

A *
max as a function of Re in the laminar regime The variation of the peak amplitude as a function of the Reynolds number can be estimated using the following assumptions. The mean drag coefficient C d and the inverse of the Strouhal number 1/S for stationary bluff bodies display similar variations as functions of Re so that the product SC d remains almost constant (Alam & Zhou 2008). For vibrating cylinders, there is also evidence that the corresponding product (f * /2)C d remains constant (Griffin 1978). Alam & Zhou (2008) showed that SC d ≈ 0.25 for stationary bluff bodies of various cross-sectional shapes throughout the sub-critical regime. However, SC d appears to increase slightly with Re in the laminar regime as seen in their data compilation (see their figure 5). We have confirmed this dependency on Re in the laminar regime for the stationary as well as the freely-vibrating circular cylinder as shown in table 7. Nevertheless, the variation of SC d and (f * /2)C d with Re is small, and in order to keep some generality in the result, we assume that f * C d remains approximately constant. In this case, equation (4.9) shows that A * max ∝ C dv . In addition, we assume C dv is magnified by a constant factor at resonance so that A * max ∝ C dv0 , where C dv0 is the fluctuating drag coefficient for a stationary cylinder. Then, peak amplitudes at different Re can be estimated from data at one particular Reynolds number Re 0 , i.e.
(4.10) Figure 12 shows predictions of A * max using C dv0 data at different Re (taken from Qu et al. 2013) and the known response at Re = 180. It can be seen that A * max predictions fit well with the amplitude from the simulations at Re = 100 while they also extrapolate to the result at Re = 250. The enhancement of A * max with Re is close to quadratic, which is quite marked compared the almost constant amplitude of purely transverse vortex-induced vibration in the laminar regime.

Re
SC d (f * /2)C d 100 0.217 0.216 180 0.252 0.244 Table 7. Compilation of data for stationary and vibrating circular cylinders in the laminar regime. Data for the stationary cylinder are from the study of Qu et al. (2013) and for the freely vibrating cylinder are from the present simulations at peak amplitude.

Conclusions
In this study we have developed a theoretical model for the fluid force acting on a circular cylinder vibrating in-line with a free stream. The streamwise fluid force comprises the sum of three components: an inviscid added-mass force opposing the inertia of the cylinder, a quasi-steady drag force opposing the velocity of the cylinder, and a vortex drag force. The key element here is that the viscous force is split into a quasi-steady drag and a vortex drag, which can be thought of as separate contributions from the vorticity bounded to the cylinder and that shed in the wake.
The theory enabled us to decipher three important fluid-dynamical aspects. First, the phase angle of the vortex force with respect to the cylinder displacement φ dv increases smoothly with U * due to a gradual shift in the timing of vortex shedding. This contrasts the variation of the phase angle between the total force and the cylinder displacement, which must remain fixed at φ x = 0 • for U * f * < 1 and at φ x = 180 • for U * f * > 1, with a sudden jump at the coincidence point (U * f * = 1); this constraint is imposed by the equation of cylinder motion when the structural is null and has no bearing on the fluid dynamics. Second, the added mass coefficient for inviscid flow is also applicable for viscous flow about a cylinder oscillating in-line with a free stream. Third, the vortex drag is amplified in the excitation region of relatively high-amplitude response, which supports the classification of vortex-induced in-line vibration as a resonance phenomenon.
The above findings could be confirmed because, unlike previous works dealing with vortex-induced vibration that primarily consider the fluid force in the direction of cylinder motion, here we calculated the phase angle between the force transverse to the direction of motion and the cylinder displacement, φ y . This allowed us to illustrate that the smooth variation of φ y as a function of U * , in contrast to φ x that displays a sudden jump by 180 • , is directly linked to a gradual shift in the timing of vortex shedding.
The theory developed in this investigation is linear, which means that all force components can be well approximated by single-harmonic functions of time, although the fluid dynamics is non linear by default. Departures from linearity can arise at much higher amplitudes and/or frequencies of cylinder vibration because of the quadratic relative velocity term. In addition, other non linearities may arise due to competing requirements posed by the structural dynamics and the fluid dynamics. In fact, we have pinpointed in the present study that such non-linear effects become apparent at the 'coincidence point' where the vibration frequency becomes equal to the natural frequency of the structure in vacuum. At this point, the equation of cylinder motion requires that the main component of the in-line force at the main harmonic of the vibration must become null. In this case, the component at the first super-harmonic of the vibration frequency dominates the driving in-line force.
An important result predicted from the theory is that the response does not depend on m * at the point where φ dv = 90 • . The simulations show that this occurs at the 'resonance point' where the same peak amplitude is attained for all m * values, in accord with the theory. It should be noted however that the theory cannot predict that the maximum amplitude occurs when φ dv = 90 • , or that this operating point does materialize, since the solution is given in open form. On the other hand, the flow physics suggest that at some point φ dv = 90 • due to the gradual shift in the timing of vortex shedding, which accompanies the continuous variation of the vibration frequency as the velocity is varied. Nevertheless, we have also observed that discontinuities may arise by changing some parameter, e.g. Re or m * . Further study over wider ranges of these parameters than considered here may worth undertaking in the future.
The simulations show that Reynolds number has a remarkable effect on the maximum amplitudes A * max attained over the entire reduced velocity range; increasing Re from 100 to 250 results in a 12-fold increase of A * max . This stands in sharp contrast to the free vibration transversely to a free stream, in which case peak amplitudes remain fairly constant in the laminar regime. This may be attributable to variations of the added mass coefficient in the latter configuration (Konstantinidis 2013). In view of this complexity, we consider that free in-line vibration offers a more convenient test case to uncover the fluid dynamics. The fluid dynamics could be elucidated because in-line response amplitudes remain very small for the low Reynolds numbers investigated in the present study. As a consequence, the fluid excitation comes solely from the primary wake instability associated with alternating vortex shedding, which remains robust and similar as in the wake of a non-vibrating cylinder. It has been well established in the published literature that other instabilities can be excited, even at small amplitudes, with increasing the Reynolds number, such as the symmetrical mode of vortex shedding. There may also exist competition between different modes. Further instabilities due to gradual transition to turbulence will inevitably perplex the phenomenology of vortexinduced vibration. However, we maintain that the theory developed here remains valid and can be used to analyse the more complex phenomena at higher Reynolds numbers, possibly with some adjustments for different fluid excitation mechanisms.  Table 8. List of non-dimensional parameters employed in the present study.
where t = t − t 0 is the time with origin at a particular instant t 0 during the simulation where the displacement oscillation is at peak after oscillations have stabilised to a steady state. It should be noted that in the present simulations the frequency of alternating vortex shedding and frequency of the unsteady transverse force are both equal to half the frequency of cylinder oscillation (sub-harmonic synchronization). In the following, the procedure to calculate the phase angle φ y of the transverse force with respect to displacement is shown in detail. Expansion of the cosine term in (B 3) yieldsC y (t) = C y cos φ y cos (πf t ) − C y sin φ y sin (πf t ). MultiplyingC y (t) by cos (πf t ) and then taking the time average yields C y (t) cos (πf t ) = C y cos φ y cos 2 (πf t ) − C y sin φ y sin (πf t ) cos (πf t ), (B 4) where overlines denote time averaging. Similarly, multiplyingC y (t) by sin (πf t ) and then taking the time average yields C y (t) sin (πf t ) = C y cos φ y cos (πf t ) sin (πf t ) − C y sin φ y sin 2 (πf t ).

(B 5)
The time average of sin (πf t ) cos (πf t ) over an even integer number of cycles is zero by default while cos 2 (πf t ) = sin 2 (πf t ). Thus, substitution of these values into (B 4) and (B 5) and then combining them yields φ y = tan −1 −C y (t) sin (πf t ) C y (t) cos (πf t ) .

(B 6)
A similar procedure can be employed to compute the phase angle φ x of the in-line force by multiplyingC x (t) in Eq. (B 3) with cos (2πf t ) and independently with sin (2πf t ), averaging both equations and then combining them, which leads to the following result φ x = tan −1 −C x (t) sin (2πf t ) C x (t) cos (2πf t ) .
(B 7) The accuracy of the above procedure is limited by the time step employed in the simulations, which is very small and thus leads to negligible errors in the calculation of the phase angles.