1 Introduction

Vortex shedding is often observed in engineering,such as circular cylinders in cross-flow^[1],oscillating airfoils and turbomachine blades^[2],and micro air vehicles^[3]. Lai and Platzer^[2] studied different vortex patterns shedding from the trailing edge of an NACA 0012 airfoil,and found that the vortex patterns oscillated sinusoidally in the plunge were captured in water tunnel tests with a Reynolds number range from 500 to 21 000 based on the airfoil chord. Gostelow et al.^[4] showed a strong similarity between the vortex wakes shed from cylinders and airfoils with the sinusoidal plunge motion in the low-speed flow and the wakes shed from a turbine nozzle cascade in the transonic flow,and presented the thermo-acoustic effect associated with the vortex shedding in the rotating machines. Lawaczek and Heinemann^[5-6] found a turbine blade with a blunt trailing edge to shed a von Kármán vortex street. Aeroelastic problems may occur easily in operating the conditions owing to the vortex shedding,such as the flutter for the aircraft wings,helicopter,and turbine blades.

For a fluid-structure interaction (FSI) problem,the vortex-induced vibration is difficult in modelling the coupled loads. Because both experimental^[7] and numerical^[8-11] approaches are expensive to obtain robust results,simplified models were developed to reveal the nonlinear effects of the system parameters on the dynamic responses of the system. Generally,turbine blades are simplified as the elastic clamped-free beams in both experimental^[12] and numerical analyses^[13-16]. Cao et al.^[17] and Chu et al.^[18] analyzed the two-dimensional friction contact problem and impact vibration characteristics of a rotational blade in the centrifugal force field,and simplified the blade as a continuous cantilever.

The van der Pol model for the FSI problem have been widely studied. The results show that the self-excited oscillation for the structure owns to the characteristic of the time-varying^[19-21]. Lee et al.^[22] used a van der Pol oscillation to model an aeroelastic system possessing limit cycle oscillations. Specifically,Hartlen and Currie^[23] introduced a van der Pol-based model,which captured many of the features seen in experimental results. Skop and Griffin^[24] subsequently modified and improved the van der Pol-based model. The vortex-induced vibrations for a cylinder^[25],an offshore riser^[26],and a turbine blade^[27] have also been analyzed,where the time-varying characteristics of the vortices are modelled by the van der Pol oscillation,and the effects of the structural motion on the fluid are studied. Barron and Sen^[15] added a van der Pol damping term to the equations of four coupled elastic beams to represent the self-excitation and fluid-structure interaction of the system. Barron^[16] simulated the FSI problem for turbine blades by adding a van der Pol self-exciting term to a single partial differential equation of a linear beam. Wang et al.^[28] used the van der Pol oscillator to simulate the time-varying characteristics of the lift coefficient for the fluid-structure interactions of turbine blades.

Moreover,it has been demonstrated that the vortex shedding from oscillating airfoils and cylinders can be significantly affected by the body oscillation^{[1, 25-26, 29-33]}. Therefore,to analyze the vibrations of the vortices further,it is necessary to consider the action of the structural motion on the fluid.

The motivation of this paper is to investigate the dynamic response and bifurcation characteristics for the vortex-induced vibrations of a rotating blades. The blade is modelled as a cantilever beam,while a van der Pol equation is used to mimic the time-varying characteristics for the vortices. The reaction of the structural motion on the fluid is considered and represented by a linear inertial coupling. The 1:1 internal resonance analysis is carried out with the multiple scale method. The two-parameter bifurcation diagram is derived by the singularity theory,and the frequency-response curves in different parameter regions are obtained. The time histories,phase portraits,and Lyapunov exponents are calculated by the Runge-Kutta method for the original system to verify the validity of the multiple scale method.

2 Mathematical modelling 2.1 Cantilever beam model for blade

In engineering practice,the rotating Euler-Bernoulli beam^[34-35] is often used to analyze the dynamic characteristics of engineering systems,such as turbomachinery,wind turbines,robotic manipulators,and rotorcraft blades. To investigate the interaction mechanism of the structure and fluid,the blade is assumed to be a continuous uniform straight cantilever beam based on the Euler-Bernoulli formulation in the centrifugal force field (see Fig. 1).

The transverse motion of the beam will be mainly studied,and the effects of the shear deformation and rotatory inertia are neglected. özgür and Gökhan^[35] studied the in-plane vibrations of a rotating Euler-Bernoulli beam with the Coriolis force,and revealed that all the coefficients related to the Coriolis force would vanish during the reduction processing,where a single degree-of-freedom model for the ith mode vibration was extracted by the Galerkin method. Therefore,there is no Coriolis term in the single degree-of-freedom equation.Xu et al.^[36] demonstrated that the Coriolis force had no influence on the axial and torsional vibrations of the rotating blades,and the effect on the tangential vibration was very small and could be neglected. Therefore,the effect of the Coriolis force on the rotating blade is neglected in this study.

The transverse displacement of the cantilever beam is assumed as w(x,t),and the equation of the blade motion can be derived by considering the equilibrium of the forces and moments acting on the differential segment of the blade with the length of dx (see Fig. 2). For more details,please refer to Ref.^[17].

The first dynamic equilibrium relationship with respect to the transverse displacement w(x,t) can be obtained by summing all the forces in the transverse direction for the segment,i.e.,

where m=(p+p_f)A is the total mass consisting of the structure and the added mass induced by the fluid,p and p_f are the densities of the structure and air flow,respectively,and A is the area of the cross-section for the cantilever beam. p and A are constants according to the beam assumption. Q is the shear force acting on the cross-section of the blade. c is the viscous damping coefficient. f is the axial load expressed as follows (see Fig. 1(b)):

where L is the length of the blade,and Ω is the rotating speed of the blade. Ff is the fluid force acting on the blade induced by the vortices,i.e.,

where C_L is the lift coefficient representing the time-varying characteristics of the vortices,and U is the total velocity defined by in which V is the freestream velocity and D_0 is the characteristic length of the cross-section of the cantilever beam.

Neglecting the shear deformation and the rotation of the cross-section according to the beam assumption,when the condition of moment equilibrium is satisfied,we can introduce the basic moment-curvature relationship as follows:

where M is the moment acting on the cross-section of the beam,and EI is the flexural rigidity of the structure. With Eq.(2),we can simplify Eq.(1) as follows:

Substitute the following non-dimensional variables:

into Eq.(3). Then,we have

The boundary conditions of the cantilever beam should meet the following conditions:

(i) When z=0,the displacement and rotation of the beam should be equal to zero,i.e.,

(ii) When z=1,the moment and shear force should be equal to zero,i.e.,

2.2 Van der Pol oscillator

The van der Pol oscillator is often used to model the time-varying and self-sustained characteristics of the flow in the analytical/experimental research of vortex-induced vibrations. The lift coefficient can represent the flow characteristics^[25-26],which can affect the lift force and then the structural vibrations. Herein,the van der Pol oscillator is introduced to simulate the time-varying characteristics of the lift coefficient,i.e.,

where q(z,t)=2 C_L/C_L0,in which C_L0 is the reference lift coefficient,w_f=2π StV/D₀ is the shedding frequency of the vortex,St is the Strouhal number,and λ is the van der Pol damping coefficient. Moreover,the action of the structural vibration on the fluid motion F_S is considered to be outlined in the introduction. It has been shown in Refs.^[25] and ^[26] that the inertial coupling is the ideal form to describe the reaction of the blades to the flows. Therefore,the linear inertial coupling is studied,and the force F_S induced by the structural motion can be assumed as follows:

where N~ is the linear coupling parameter.

To analyze in the same time scale,let T =Ω t. Then,Eq.(6) becomes

Therefore,the coupled equations (5) and (7) model the interactions of the blade and vortices.

2.3 Reduced model with Galerkin discretization

An arbitrary oscillation of the structure v(z,T ) can be accurately described with a sum of all of its modal responses,which can be mathematically expressed with a Taylor series composed of the individual modal components^[37] as follows:

where v~_i (z) and v_i (T) (i=1,2,...) represent the mode shapes and the modal coordinates of the beam,respectively. When the shapes of the deformation are known,the modal coordinates v_i (T) (i=1,2,...) define how the amplitudes of the associated deformations v~_i (z)(i=1,2,...) change with time. Thus,the continuous system can be discretized by defining it on the modal subspaces.

For the approximate solution,the orthonormal sets of the eigen amplitude functions for a cantilever beam as the modal functions for the blade are introduced,i.e.,

where β_i(i=1,2,...) are the roots of the transcendental equation

which is derived from the boundary conditions of the cantilever beam.

Herein,the cross-section of the beam is assumed to be uniform. Therefore,the orthogonality of different modal functions can be written as follows:

It means that the spatial distributions of the structure are perpendicular to each other in the modal space if they are not of the same order.

When the modal functions are in the same order,the integral can be calculated as follows:

It can be seen from Eq.(11) that the value of the integral can be determined by the corresponding coefficient β_i.

The span-wise distribution of the shedding variable has been widely studied^[38-43],in which the spatial component of q(z,T ) is described with the shape of the respective normal mode and separated from the temporal part of the response. Hence,the form of the solution q(z,T ) can be assumed to be expressed with a Taylor series of modal components as analyzed in Refs.^[39]-^[43],i.e.,

where the relations between the mode-shapes and the modal coordinates q_i (T )(i=1,2,...) are the same as those in Eq.(7).

The modal functions for the van der Pol oscillator are introduced according to Ref.^[43],i.e.,

which satisfy the following orthogonality condition:

With the relation

substituting Eq.(8) into Eq.(5) and orthogonalizing the latter with respect to the set v~_j (z),we have

Because the fluid-structure interaction problem considered in this study is an ideal one,the fluid force can induce the structural oscillation in only one mode. Therefore,the van der Pol oscillator can be expected to have the same and single modal distribution in the space,and the first modal approximations of the cantilever beam and the van der Pol oscillator can be investigated.

The first mode equation of the structure can be derived from Eq.(15) when i=j=1,i.e.,

Thus,the first-order equation for the cantilever beam can be simplified as follows:

where w_R =w f/Ω is the non-dimensional frequency of the fluid, is the non-dimensional frequency of the structure,and . In the above equations,

The values of the above definite integrals can be obtained when the coefficient β₁ is determined.

Similarly,substituting Eq.(12) into Eq.(7) and orthogonalizing the obtained results with respect to q~_j (z),we have

Furthermore,it can be rewritten as follows:

where

When m=n=p=i=1,Eq.(20) can be simplified as follows:

When i=j=1,we introduce

as the inertia coupling parameter,which denotes the action for the first mode motion of the structure on the first mode motion of the van der Pol oscillator. Thus,the first mode motion of the van der Pol oscillator can be derived as follows:

Equations (17) and (22) model the interactions between the structural vibration and the van der Pol oscillation,and will be investigated in the following sections.

3 Bifurcation analysis 3.1 1:1 internal resonance analysis with multiple scale method

When the resonance occurs^[44-45],such as the primary resonance,the internal resonance,and the superharmonic/subharmnic resonance,the nonlinear systems display the rich dynamic characteristics. The multiple scale method is a good way for understanding the qualitative characteristics of the system presenting the internal resonance^[46-48]. The phenomenon of lock-in or synchronization is often analyzed in the vortex-induced vibration for the long cylinders (e.g.,the offshore risers),which means that the vortex shedding frequency fV tends to the natural frequency of the structure f_S ^[29]. Similarly,the phenomenon for the frequency approximation can also be encountered in the coupled system proposed in this paper. For Eqs. (17) and (22),the non-dimensional frequencies w_S and w_R of the first-order coupled system (see Fig. 3) are approximate to each other around the rotating speed Ω≈ 450 rad·s⁻¹.

Hence,the 1:1 internal resonance of the coupled system is interesting,and will be investigated by the multiple scale method. The relation between the non-dimensional frequencies w_S and w_R can be written as follows:

where σ is the detuning parameter. Introducing the scaling parameters

into Eqs.(17) and (22),we have

Assume the approximate form of the solutions as follows:

Then,substituting Eq.(24) into Eq.(23) and equating the coefficient of like powers of ε,we have the results of the order ε ⁰ as follows:

and the order ε ^1 as follows:

where

and is the partial differential operator. Then,we can obtain the solutions of Eqs.(25) and (26) as follows:

Substitute Eqs.(29) and (30) into Eqs.(27) and (28). Consider the internal resonance condition. Then,we have

where c.c. stands for the complex conjugate of the proceeding terms.

The solvable conditions of Eqs.(31) and (32) are derived by equating the coefficients of the secular terms to be zero,i.e.,

The derivatives of the amplitudes A and B with respect to T₁ can be obtained by Eqs.(33) and (34) as follows:

The functions A and B can be expressed in the polar coordinates as follows:

where a_j and θ_j(j=1,2) are the amplitudes and phase angles,respectively. Substituting Eq.(37) into Eqs.(35) and (36) yields the following first-order differential equations after separating the real and imaginary parts:

where (′) denotes the derivatives with respect to T₁,and φ = θ ₂ - σ T₁ - θ ₁. The derivative of φ with respect to T₁ can be obtained by eliminating θ₁ and θ_2 from Eqs.(39) and (41) as follows:

The equilibrium solutions of Eqs.(38),(40),and (42) correspond to the periodic motion of the coupled system. To obtain the equilibrium solutions,a′_j(j=1,2) and φ′ in Eqs.(38),(40),and (42) are assumed to be equal to zero. Then,we can obtain the frequency-response equations as follows:

Equations (43) and (44) can reveal the effects of the system parameters on the responses.

3.2 Singularity analysis for steady-state responses

The structural motion and the van der Pol damping have important effects on the vibration of the fluid. Then,the fluid motion in turn can affect the structural vibration. To investigate the bifurcation characteristics of the coupled system in a wider parameter space,an engineering unfolding analysis is carried out (see Refs.^[49]-^[51] for the details). Equations (43) and (44) can be rewritten as follows:

where σ is the bifurcation parameter,λ and N are the engineering unfolding parameters,and a₁ and a₂ are the state variables. The transition set is derived by the singularity method as follows^[52]:

where B,H,and D represent the bifurcation set,the hysteresis set,and the double limit set,respectively. The partial derivatives of the bifurcation functions G₁ and G_2 with respect to σ,a₁,and a₂ are calculated to obtain the bifurcation sets (see Appendix A for details). Therefore,the transition set of the coupled system can be derived by solving the algebraic equations shown in Eqs.(48)-(50). According to the analysis in Refs.^[17],^[25],and ^[43],the other parameters are fixed as follows:

It can be seen from Fig. 4 that the two-parameter space is divided into twelve parts by the transition set,where the red dashed lines denote the bifurcation set,and the blue solid lines present the hysteresis set. Herein the double limit set is D=ø after calculation.

The transition set can be used to classify different kinds of responses of the system. The responses of the coupled system can display different dynamic characteristics when the van der Pol damping λ and the coupling parameter N are chosen in different parameter subspaces.

The representative frequency-response curves of the system in each region of the two-parameter space shown in Fig. 4 are calculated,and the stability is determined by examining the eigenvalues of the corresponding characteristic equation for Eqs.(38),(40),and (42). The characteristic equation is a triple polynomial as follows:

A solution is stable if all of the eigenvalues have negative real parts. It is decided by the Routh-Hurwitz criterion^[53] as follows:p₁ ＞0,p₃＞0,and p₁ p₂-p₃＞0.

Figure 5 shows that the amplitudes a₁ and a₂ of the steady-state solutions are unstable for λ=-0.0639 and N=0.1836 in Region I,and the varying trends of the amplitudes a₁ and a₂ with respect to the detuning parameter σ are opposite to each other.

Figure 6 depicts that the trivial solutions of the amplitudes a₁ and a₂ lose stabilities via a saddle-node bifurcation at S_N1,resulting in the occurrence of a two-mode solution. One of the two-mode solution is unstable for both the amplitudes a₁ and a₂,which has opposite varying trends for the amplitudes a₁ and a₂ with respect to the detuning parameter σ. The other one of the two-mode solution is unstable until a Hopf bifurcation occurs at H₁,resulting in a change of the solutions from unstable to stable. The stable solutions decrease until encounter another saddle-node bifurcation at S_N2,leading to a change of the solutions from stable to unstable. Then,the unstable solutions decrease until another Hopf bifurcation occurs at H₂,resulting in a change of the solutions from unstable to stable.

Figure 7 shows that,as the detuning parameter σ increases,the trivial solutions jump to the large unstable solutions via a Hopf bifurcation at σ=0.0183 (see H). Then,the unstable amplitude a₁ decreases while a₂ increases as σ increases. As σ decreases,the unstable solutions jump to the trivial ones via a saddle-node bifurcation at σ=0.0181 (see S_N),leading to a change of the solutions from unstable to stable.

It can be seen from Fig. 8 that as σ increases,the responses increase until encounter a saddle-node bifurcation at S_N1,resulting in that the responses jump to the large solutions. The large responses increase until reach the maximum values,then decrease until a saddle-node bifurcation occurs at S_N2,leading to that the solutions jump to the small solutions. Then,the responses decrease all the way as σ increases. Similarly,as σ decreases,the amplitudes a₁ and a₂ have the same varying trends,and encounter two saddle-node bifurcations at S_N3 and S_N4,

Figure 9 shows that the amplitudes of the steady-state solutions of the beam and van der Pol oscillator are stable for λ=0.0400 and N=0.0519. Both Fig. 8 and Fig. 9 show that the varying trend of the amplitude a₁ with respect to the detuning parameter is the same as that of the amplitude a₂.

It can be seen from Fig. 10 that the stable trivial solutions become unstable via a Hopf bifurcation for σ=0.0090.

Figures 11-13 show that the solutions in Regions VII-IX are unstable solutions expect the trivial solutions. While the varying trends of the amplitudes a₁ and a₂ with respect to the detuning parameter σ are opposite in Fig. 11 and the same in Figs. 12 and 13.

It can be seen from Fig. 14 that the trivial solutions jump to the large stable solutions via a Hopf bifurcation at σ= -0.0015 (see H₁) as σ increases. Then,the amplitude a₁ decreases and the amplitude a₂ increases to the trivial solutions as σ increases. While as σ decreases,the amplitude a₁ increases and the amplitude a₂ decreases until a Hopf bifurcation occurs at σ=-0.0052 (see H₂) ,leading to a change of the solutions from unstable to stable. When σ decreases beyond H₂,the amplitudes a₁ and a₂ jump to the trivial solutions.

It can be seen from Fig. 15 that as σ increases,the amplitude a₁ grows while the amplitude a₂ decreases all the way until a Hopf bifurcation occurs at σ=0.0045 (see H₁),after that,the responses encounter a saddle-node bifurcation at σ=0.0064 (see SN₁),resulting in a change of the solutions from unstable to stable. Beyond S_{N₁,the amplitude a₁ jumps to the small solution,while a₂ jumps to the large solution. Then,the amplitude a₁ decreases while a₂ increases all the way as σ increases. The amplitudes a₁ and a₂ have symmetrical varying trends as σ decreases.

Figure 16 shows that the trivial solutions jump to the stable solutions via a Hopf bifurcation for σ=0.0003. The amplitude a₁ jump to the maximum values rapidly. Then,as σ increases,the amplitude a₁ decreases while the amplitude a₂ increases gradually until to the trivial solutions.

The results show that when λ and N have the same signs,i.e.,λ N ＞0,the varying trends for the structural and van der Pol vibrations are the same with respect to the detuning parameter σ. Specifically,when λ＞0 and N＞0,the solutions are stable or can encounter the saddle-node bifurcation as σ varies; when λ ＜0 and N＜0,the solutions are unstable as σ varies. When λ and N have opposite signs,i.e.,λ N＜0,the varying trends of thetwo coupled motions are reverse to each other,which means the energy transfer between the two modes. The Hopf bifurcation can occur for certain parameter valueswhen λ N＜0.

4 Numerical results

The response bifurcations of the coupled system including the Hopf and saddle-node bifurcations can occur under certain parameter values as analyzed in the previous section. To verify the validity of the multiple scale method,the characteristics of the dynamic responses for the original systems (17) and (22) are investigated. A representative point (σ=0.0285) after the Hopf bifurcation shown in Fig. 6 is chosen with the unfolding parameters λ=-0.0372 and N=0.3720. The other parameters are fixed as follows:

The time histories,phase portraits,and Lyapunov exponents for the original structural and van der Pol motions with the above parameter values are computed by the Runge-Kutta method (see Figs. 17 and 18). It shows that the largest Lyapunov exponent of the original system with the initial condition values v₁(0)=0.05,y(0)=0.00,q₁(0)=0.50,and p(0)=0.00 (assuming v₁(T)=y(T) and q₁(T)=p(T)) tends to zero as time T increases,which means that the quasi-periodic solutions occurs under this parameter conditions. Therefore,the Hopf bifurcation for the modulated solutions indicates the occurrence of the quasi-periodic solutions for the original system,which can result in important physical consequences^[54],such as a behavior from transition to chaos.

The saddle-node bifurcation is another bifurcation type analyzed in this paper,which indicates the stability changing of the responses,resulting in the occurrences for the jump phenomenon and multiple solutions. A typical point (σ=0.016) between the two saddle-node bifurcation points when λ=0.0629 and N=0.2178 in Region IV (see Fig. 8) is chosen to obtain two stable solutions by the Runge-Kutta method. The other parameters are fixed as follows:

The two sets of the time histories and phase portraits for the structural and van der Pol motions are obtained with theinitial condition values v₁(0)=0.1,y(0)=0.0,q₁(0)=2.0,and p(0)=0.0 and v₁(0)=0.40,y(0)=0.00,q₁(0)=2.65,and p(0)=0.00,respectively,which are shown in Fig. 19.

Moreover,the time histories of the two-degree-of-freedom motions can be predicted by Eqs.(24),(29),and (38)-(41) analytically. The comparisons of the time histories obtained by the analytical and numerical methods for the two different solutions shown in Fig. 19 are carried out with the same system parameter values as follows:

Figures 20 and 21 show that the analytical results for the time histories agree with the numerical results obtained by the Runge-Kutta method.

5 Conclusions

The vortex-induced vibrations of a rotating blade are investigated. The blade is modelled as a uniform and straight cantilever beam,and the van der Pol oscillator is used to simulate the time-varying of the vortex. The reaction for the motion of the blade on the fluid is represented by a linear inertial coupling. The multiple scale method is used to analyze the 1:1 internal resonance of the coupled system. The bifurcation equations are derived,and a two-parameter bifurcation diagram for the van der Pol damping λ and the coupling parameter N is obtained with the singularity theory of two state variables. The bifurcation characteristics for the frequency-responses in different bifurcation regions are investigated. The phenomena including the saddle-node and Hopf bifurcations are found to occur under certain parameter regions. The bifurcation analysis shows that when the parameters λ and N have the same signs,i.e.,λ N>0,the varying trends of the structural and van der Pol motions are similar to each other. Specifically,when both λ and N are positive,i.e.,λ>0 and N>0,two kinds of dynamic characteristics for the responses can occur,i.e.,(i) the responses are totally stable as σ varies,e.g.,the responses in Region V; (ii) the solutions can encounter the saddle-node bifurcation at a certain value for the detuning parameter σ,e.g.,the responses in Region IV. When both λ and N are negative,i.e.,λ ＜0 and N＜0,the solutions are unstable as σ varies,e.g.,the responses in Region VIII. When λ and N have opposite signs,i.e.,λ N＜0,the varying trends of the two motions are opposite to each other,which means the energy transfer between the two modes. The Hopf bifurcation can be encountered for certain parameter values when λ and N have opposite signs. The time histories,phase portraits,and Lyapunov exponents are calculted with the Runge-Kutta method for the original system at a point after the Hopf bifurcation in Region III. The results reveal that the occurrence of the Hopf bifurcation for the modulation equations indicates the quasi-periodic motion for the original system. Additionally,the coexisting multiple solutions generated for the saddle-node bifurcation in the parametric Region IV are illustrated,whilst their corresponding numerical and analytical time histories are compared subsequently,both of which are in good agreement with each other. These results indicate the validity of the analytical solutions obtained by the multiple scale method.

Acknowledge The first author thanks for the hospitality of the University of Aberdeen.