Magneto-elastic dynamics and bifurcation of rotating annular plate?

Yu-Da Hu(胡宇達(dá)),Jiang-Min Piao(樸江民),and Wen-Qiang Li(李文強(qiáng))

1 School of Civil Engineering and Mechanics,Yanshan University,Qinhuangdao 066004,China

2 Key Laboratory of Mechanical Reliability for Heavy Equipment and Large Structures of Hebei Provincial, Yanshan University,Qinhuangdao 066004,China

Magneto-elastic dynamics and bifurcation of rotating annular plate?

Yu-Da Hu(胡宇達(dá))1,2,?,Jiang-Min Piao(樸江民)1,2,and Wen-Qiang Li(李文強(qiáng))1,2

1 School of Civil Engineering and Mechanics,Yanshan University,Qinhuangdao 066004,China

2 Key Laboratory of Mechanical Reliability for Heavy Equipment and Large Structures of Hebei Provincial, Yanshan University,Qinhuangdao 066004,China

In this paper,magneto-elastic dynamic behavior,bifurcation,and chaos of a rotating annular thin plate with various boundary conditions are investigated.Based on the thin plate theory and the Maxwell equations,the magneto-elastic dynamic equations of rotating annular plate are derived by means of Hamilton’s principle.Bessel function as a mode shape function and the Galerkin method are used to achieve the transverse vibration differential equation of the rotating annular plate with different boundary conditions.By numerical analysis,the bifurcation diagrams with magnetic induction,amplitude and frequency of transverse excitation force as the control parameters are respectively plotted under different boundary conditions such as clamped supported sides,simply supported sides,and clamped-one-side combined with simply-another-side.Poincaré maps,time history charts,power spectrum charts,and phase diagrams are obtained under certain conditions, and the influence of the bifurcation parameters on the bifurcation and chaos of the system is discussed.The results show that the motion of the system is a complicated and repeated process from multi-periodic motion to quasi-period motion to chaotic motion,which is accompanied by intermittent chaos,when the bifurcation parameters change.If the amplitude of transverse excitation force is bigger or magnetic induction intensity is smaller or boundary constraints level is lower,the system can be more prone to chaos.

magneto-elastic,rotating annular plate,Bessel function,bifurcation and chaos

1.Introduction

In the fields of aerospace,nuclear industry,large generator,computer storage devices,turbine engines,and electromagnetic sensor device,the rotating annular plate as the basic component is widely utilized.Especially transverse vibration of a high-speed rotating circular plate in conductive electromagnetic environment will be unstable when subjected to small disturbance,and even lead to huge harm to electromechanical system.Thus,it is most significant to study magnetoelastic dynamics of rotating conductive annular plate.

Leo et al.[1]presented approximately linear dynamic differential equation of a circular plate with multiple circular holes by means of Rayleigh method,and used Bessel functions to get the vibration response equation of the system.Based on shear deformation theory of circular plate,Malekzadeh et al.[2]obtained vibration equation of the annular plate in thermal environment,and analyzed the effects of temperature,inner diameter,outer diameter,and other parameters on natural frequency of the annular plate.According to the linear deformation theory,Saidi et al.[3]derived and solved the vibration control equations in different boundary conditions. Allahverdizadeh et al.[4]got the vibration equation of rotating circular plate by using Galerkin method,and obtained more accurate solutions by means of semi-analytic differential perturbation method.Ratko[5]utilized calculus to study the transverse vibration and stability of a rotating circular plate. Hashemi et al.[6]have made the research on the transverse vibration and the stability of rotating circular plates with the finite element method.

For vibration of a conductor in electromagnetic field, many researchers have got results.Zheng et al.[7]studied bend and vibration problems of a cantilever conductive plate in magnetic field,and found that there is unstable vibration of the circular conductive plate under specific transverse excitation load.Gao et al.[8]obtained the analytical solution of a circular plate combined with Maxwell’s equations,and analyzed the resonance characteristics.Hasanyan and Librescu[9]investigated vibration of a beam-plate in magnetic field,and discussed the influence of magnetic induction intensity and boundary conditions on the free vibration by using numerical analysis.Hu and Wang[10]studied magneto-elastic free vibration of a rotating circular plate in static magnetic field. Li et al.[11]obtained the three-dimensional analytical solution for functionally graded the magneto-electro-elastic circular plates subjected to uniform load through a step-by-step integration scheme,with five integral constants determinable from the boundary conditions at the cylindrical surface in the Saint Venant’s sense.Alaimo et al.[12]and Razavi and Shooshtari[13]respectively accomplished static and free vibration analyses of the magneto-electro-elastic multilayeredplates based on the first-order shear-deformation plate theory.

Bifurcation and chaos are very common issues in research and application of many fields.Lu et al.[14]investigated dynamic stability and bifurcation of an alternating load and magnetic field excited magnetoelastic beam,and obtained the solutions.Hu and Zhang[15]deliberated the bifurcation characteristics of axially moving plate in magnetic field,derived bifurcation equation and transition set corresponding to the universal unfolding,and discussed the impacts of bifurcation control parameters on bifurcation and chaos.For a functionally graded circular plate,Hu et al.[16]investigated unfolding problems of bifurcation equation,and plotted bifurcation diagrams.Hu and Zhang[17]analyzed bifurcation of the circular functionally graded plate with combination resonances.In response to geometrically nonlinear problem of a circular plate, Touzé et al.[18]derived ordinary differential equations by using Galerkin method,and plotted the bifurcation diagrams and Poincaré maps.Coman[19]analyzed the influence of initial tension on the non-axisymmetric bifurcation of the circular plate system.Shahverdi and Khalaf i[20]investigated the bifurcation problems of functionally graded plate under hypersonic aerodynamic loads.Zhao and Zhang[22]analyzed bifurcation and chaos for aeroelastic airfoil with freeplay structural nonlinearity in pitch by using the Poincaré map method.

In this paper,Bessel functions are used as vibration mode functions to solve magneto-elastic vibration equation.The bifurcations diagrams,response charts,power spectrum charts, phase diagrams,Poincaré maps,and time history diagrams of the system are explored in different bifurcation control parameters,such as magnetic induction intensity,excitation amplitude and frequency.The influence of different control parameters and boundary conditions on the bifurcation and chaos are discussed.

2.Fundamental equations

Consider a isotropic conductive annular plate rotating in magnetic field with magnetic field intensity B0z,which rotates at a constant angular speed ?.The plate has a uniform thickness h,outer radius b,and inner radius a.The Young’s modulus,the Poisson ratio,and the density of the disk are E,μ,and ρ,respectively,as shown in Fig.1.

Fig.1.Rotating annular plate in magnetic field.

2.1.Kinetic energy and deformation potential

A fixed cylindrical coordinate system(r,θ,z)is used in the model.The displacement vector of any point in the plate is assumed as

where urand uθare the radial and circumferential displacement in the middle plane,respectively;w is the transverse displacement;u1=??w/?r and v1=??w/(?θr)are angular displacements;r,θ,and z are radial,circumferential,and normal coordinates,respectively;and t is the time variable.

The kinetic energy of the rotating annular plate can be obtained as

The deformation potential caused by bending deformation of the plate can be given as

where Mrand Mθare bending moments,Mrθis torque,κrand κθare curvatures,and κrθis torsion.

The potential energy of strain in the middle plane of the plate can be given as

where Nθand Nrθare internal forces in the middle plane,εr, εθ,and εrθare strains in the middle plane.

2.2.External virtual work

It is assumed that the plate is loaded in a transverse excitation P.Thus,the expression of external virtual work δUPis written as

where δw is the transverse virtual displacement.

2.3.Virtual work of the electromagnetic force

Current density vector of the conductive plate in magnetic if eld is written as

where σ0,E,V,and B are conductivity,electric field intensity vector,absolute speed vector,and magnetic induction intensity vector,respectively.

Lorentz force loading on the annular plate is written as

In transverse magnetic field with magnetic induction B0z, electromagnetic force can be obtained as

where Fr,Fθ,and Fzare radial,circumferential,and normal electromagnetic force on unit area,respectively;mrand mθare radial and circumferential electromagnetic torque on unit area,respectively.

Thus,the virtual work of the electromagnetic force is written as

2.4.Magneto-elastic dynamic equation

Based on Hamilton principle

and combined with Eqs.(3)—(9),considering geometric nonlinear condition and neglecting the effect of the longitudinal displacement on the transverse vibration of the plate,the magneto elasticity vibration equation of the annular plate can be obtained as

3.Mode shape function and Galerkin method

3.1.Mode shape function

In order to solve the annular plate transverse nonlinear vibration equation,a solution of the liner transverse free vibration equation is employed to assume the transverse displacement w(r,θ,t)of the system(11).

Here,the liner transverse free vibration equation is introduced as

where DM=E h3/[12(1?μ2)]is the flexural rigidity.

Obviously,the solution of Eq.(12)can be put in the following form by the method of separation of variables:

where Wn(r)are displacement functions,which are given in terms of Bessel functions and modified Bessel functions of order n.Considering symmetric vibration problem(n=0), Wn(r)can be rewritten as

where K,A,B,C,and D are determined by boundary conditions and a normalizing condition.

Based on different boundary conditions,the mode shape functions of Eq.(14)are given respectively.

(i)Clamped boundary conditions

The boundary conditions for the annular plate,which is clamped by the collar(r=a)and rim(r=b)where the displacement and angular vanish,are given by

Substituting Eq.(14)into boundary conditions Eq.(15), the equation of boundary conditions is rewritten as

Since equation(16)has obviously linear correlation, which have non-trivial solutions,the determinant composed of the indeterminate coefficient is equal to zero,that is

Here,equation(17)is transcendental equation,which contains only an unknown variable K,and thus it has infinitely many values K.Since equation(16)is a linear system,whose coefficients(A,B,C,and D)are uncertain,a normalizing method can be used to determine these coefficients.It assumed that the maximum value of W0(r)is 1,and then it is can be normalized.

(ii)Simplified boundary conductions

The boundary conditions for the annular plate simplified at the collar(r=a)and rim(r=b)are expressed as

Introducing Eq.(14)into Eq.(18),similar to the clamped boundary condition,equations must satisfy

where

Hence,the value K and the relationship of coefficients can be obtained by solving Eq.(19),and then the value of the coefficients can be determined by normalizing condition too.

(iii)Other boundary conditions

When the boundary conditions of the annular plate are clamped-inner and simply-outer or clamped-outer and simply-inner,the normalizing method can also be utilized to determine the value of K.

3.2.Galerkin method

When there is only first order mode shape in the system, modal function may be in the form

where W(r)is the mode shape function in different boundary conditions,whose first-order can be used in following calculation.

The sinusoidal excitation force is

where P0and ω are the amplitude and frequency of the excitation force,respectively.

For axisymmetric problems,equation(11)can be reduced to

where DN=E h/(1?μ2)is the tensile stiffness.

Substituting Eq.(20)into Eq.(22),the nonlinear differential equation of a rotating annular plate in magnetic field is presented by means of Galerkin method,as

where

4.Numerical analysis of bifurcation and chaos

Here,we take some result comparison which contains a result of a free vibration mode of a ring-shape thin plate on the elastic base and a result of the finite dynamic element of rotating circular disks.

It is assumed that we do not consider the impacts of magnetic induction intensity,excitation,rotational speed and nonlinearity.We use the same parameters as Ref.[22]:the mass density ρ=7800 kg/m3,Young’s modulus E=200 GPa Poisson ratioμ=0.3,conductivity σ=3.63×107(?·m)?1,the thickness h=0.002 m,the outer radius b=1.2 m,and the ratio of outer radius and inner radius k=b/a.The natural frequency of the system varying with the ratio of outer radius and inner radius k is plotted in Fig.2.When k=3,the natural frequency is ω0=106.4960 rad/s.This result coincides with the result(ω0=106.7019 rad/s)of Ref.[22],which does not consider the impact of the base modulus.

Fig.2.The natural frequency of the system varying with the ratio of outer radius and inner radius k.

Here we do not consider the impacts of magnetic induction intensity and excitation,and use the same parameters as Ref.[23]:the mass density ρ=7800 kg/m3,Young’s modulus E=1.961 MPa,Poisson ratioμ=0.3,the thickness h=0.004 m,the inner radius a=0.1 m,and outer radius b=1.0 m.The natural frequency of this system varying with rotational speed is plotted in Fig.3,where the natural frequency of the system increases extremely slowly with the increase of rotational speed.This change of the natural frequency is similar to the results with rotational speed 0–2500 r/min in Ref.[23].

Fig.3.The natural frequency of the system varying with the rotational speed.

Consider the rotating annular plate under the sinusoidal excitation load,which is made of aluminum with the mass density ρ=2670 kg/m3,Young’s modulus E=71 GPa,Poisson ratioμ=0.34,conductivity σ=3.63×107(?·m)?1,the thickness h=0.001 m,the inner radius a=0.050 m,and outer radius b=0.150 m.

Parameter K and the corresponding mode function in different boundary conditions are separately as follows:

Introducing K=47.15895 with simply supported boundary conditions into Eqs.(18)and(19),the corresponding mode function is obtained as

Introducing K=32.195 with clamped supported boundary conditions into Eqs.(16)and(17),the corresponding mode function is obtained as

Similarly,when K=38.055,the corresponding mode function with clamped-outer and simply-inner boundary conditions is obtained as

When K=41.026,the corresponding mode function with clamped-inner and simply-outer boundary conditions is shown as

Sequentially,considering magnetic induction intensity, excitation amplitude,and excitation frequency as the bifurcation control parameters,the analysis of nonlinear dynamic behavior of the annular plate is presented.

4.1.Magnetic induction intensity as a bifurcation control parameter

Bifurcation diagram with simply supported inner–outer sides is plotted in Fig.4,where magnetic induction intensity is a bifurcation control parameter.Here,fundamental parameters are the rotation speed ?=10000 r/min and excitation amplitude P0=40 kN/m2.

Fig.4.Bifurcation diagram of the system(23)for ?=10000 r/min and P0=40 kN/m2 with simply supported inner–outer sides,whose control parameter is the magnetic induction intensity B0z.Here 4 kinds of typical regions are selected to analyze.

When the magnetic induction intensity increases,the vibration amplitude of the system is reduced,and the motion of the system is a kind of alternating-periodic motion between period doubling motion and chaotic motion.Here,within some specific regions,there are some unique phenomena,as shown in Fig.4.

As shown in Fig.5,when the magnetic induction intensity is small(Region 1 in Fig.4),there is a closed curve in the Poincaré map lattice section.Moreover,it is found that vibrational energy is mainly distributed in about 3.5 times frequency from the power spectrum chart(Fig.5(c)).Thus,there is a quasi-period motion in this system.

When the magnetic induction intensity B0z=0.6 T(Region 2 in Fig.4),as shown in Fig.6,the motion of this system converts from quasi-period motion to double-period motion, and vibration energy is mainly concentrated in the 3.5 times frequency(Fig.6(c)).

Continuing to increase magnetic induction intensity to 0.7 T＜Boz＜2.0 T(Region 3 in Fig.4),the motion of the system repeatedly converts from multi-period motion to chaos motion.The amplitude of the system is almost a constant,and the response of the system is similar,shown in Fig.7 when the magnetic induction intensity B0z=0.7 T.

Fig.5.In region 1 of Fig.4,the motion of the system is a quasi-period motion when B0z=0.5 T,which is reflected in(a)Poincaré map that is a mapping circle formed by a number of points,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.6.In region 2 of Fig.4,the motion of the system is a double-period motion when B0z=0.6 T,which is reflected in(a)Poincaré map that obtains two mapping points,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.7.In region 3 of Fig.4,the motion of the system is a chaotic motion when B0z=0.7 T,which is reflected in(a)Poincaré map that is a lattice,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.8.In region 4 of Fig.4,the motion of the system becomes a quasi-period motion again when B0z=3 T,which is reflected in(a)Poincaré map that is mapping circle formed by a number of points,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

As shown in Fig.8,when the magnetic induction intensity B0z=3 T(Region 4 in Fig.4),the motion of this system is quasi-period motion,where vibration of the system has been extremely closed to single-period motion as shown in Fig.8(b).The amplitude of the system significantly decreases, which is compared with the phase diagram and time history diagram when B0z=0.6 T(Figs.6(b)and 6(d)).

When magnetic induction intensity is a bifurcation control parameter,bifurcation diagrams with clamped supported sides(Fig.9(a))and clamped-outer and simply-inner supported sides(Fig.9(b))are also plotted.The motions of two systems are both the repeated-period phenomena from chaotic motion to multi-period motion and to chaotic motion.The chaotic region of the two systems is smaller and multi-period region is larger,compared with Fig.4.

Fig.9.Bifurcation diagram of the system(23)varying with control parameter B0z when ?=10000 r/min and P0=40 kN/m2(a)in bothclamped supported sides and(b)in clamped-outer and simply-inner supported sides.

4.2.Excitation amplitude as the bifurcation control parameter

In order to reflect the effect of excitation amplitude on the system,global bifurcation diagram with simply supported sides is obtained,when B0z=2.0 T,?=10000 r/min,and ω=1.02ω0,where ω0is natural frequency of undamped linear system.Here,the amplitude of the excitation force is a bifurcation control variable,as shown in Fig.10.Amplitude of the system increases with the increase of excitation force and the vibration form of the system is basically a kind of periodic motion when P0is small(0＜P0＜50 kN/m2).Nevertheless,when the amplitude of the excitation force is large (50 kN/m2＜P0＜170 kN/m2),amplitude of the system increases with the increase of excitation force,whose vibration form is chaotic motion mixed with multi-period motion.Then, when the amplitude of the excitation force continues to increase(P0＞170 kN/m2),amplitude of the system contains a constant,and motion of the system mainly is a chaotic motion with a small number of doubling-period motion.

Fig.10.Bifurcation diagram varying with excitation amplitude P with both simply supported sides when B0z=2.0 T,?=10000 r/min,and ω=1.02ω0,where ω0 is natural frequency of undamped linear system.

The motion of system is single-period motion when P0= 12 kN/m2,as shown in Fig.11.There are different frequencies shown in Fig.11(c),where vibrational energy is mainly concentrated in the single-frequency and successively decreases with increase of the frequency multiple.

When P0continuously increases,the motion of system is a motion from single-period motion to multi-period motion shown in Fig.12.When P0=15 kN/m2,the motion of system is quadruple-period motion.Frequency of the system is composed of multiple and 1/3 times frequency shown in Fig.12(c).

When P0=20 kN/m2,there is chaotic motion in this system,shown in Fig.13.Poincaré map,existing as a mask form, converts from several discrete points into a lattice where there is a certain geometric feature.The motion of the system begins to be a kind of irregular reciprocating motion,as shown in time history diagram(Fig.13(b)).At this point,the proportion of the single-frequency power is greater than other frequency powers,and there are some continuous impulse response nearly 1/3 times and double frequency(Fig.13(c)). The phase locus of the system is a complex and messy curve (Fig.13(d)).

When excitation amplitude is large,the motion of the system is mainly chaotic motion mixed with some periodic motions,shown in Fig.10.For example,the motion of system is chaotic motion when P0=50 kN/m2(Fig.14),but singleperiod motion when P0=64 kN/m2(Fig.15).Moreover,as shown in Figs.14(c)and 15(c),the power of the system are both centralized at the single and triple frequencies.

Fig.11.The motion of the system is a single-period motion when P0=12 kN/m2,which is reflected in(a)Poincaré map that obtains a mapping point,(b) time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.12.The motion of the system is a quadruple-period motion when P0=15 kN/m2,which is reflected in(a)Poincaré map that obtains four mapping points,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.13.The motion of the system is a chaotic motion when P0=20 kN/m2,which is reflected in(a)Poincaré map that is a mapping lattice existing as a mask form,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.14.The motion of the system is still a chaotic motion when P0=50 kN/m2,which is reflected in(a)Poincaré map that is a mapping lattice,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.15.The motion of the system is a single-period motion again when P0=64 kN/m2,which is reflected in(a)Poincaré map that obtains a mapping point, (b)time history diagram,(c)power spectrum,and(d)phase diagram.

In order to reflect the effect of magnetic induction intensity on bifurcation motion,we plot the bifurcation diagrams in different cases,where bifurcation control parameter is still the excitation amplitude P0(P0is less than that in Fig.10), shown in Fig.16.We obtain that the motions of the system are mainly multi-period motion and quasi-period motion,when magnetic induction intensity is relatively small(B0z=0 T and B0z=1.0 T).Moreover at this point the system is more likely to develop the quasi-period motion and chaotic motion by changing the excitation amplitude compared with Fig.10.

Fig.16.Brifurcation diagram varying with excitation amplitude P0 for(a) B0z=0 T and(b)B0z=1 T with ?=10000 r/min and ω=1.02ω0,where ω0 is natural frequency of undamped linear system.

In order to analyze the effect of force on bifurcation with the other boundary conditions,the bifurcation diagrams varying with excitation amplitude are plotted with clamped supported conditions(Fig.17(a))and with clamped-outer and simply-inner boundary conditions(Fig.17(b)).In comparison with the diagrams in Figs.4 and 11(b),the system subjected to nonlinear effect will be smaller and less prone to show the chaotic motion,as there is a clamped side in the system.

4.3.Excitation frequency as the bifurcation control parameter

Global bifurcation diagram(?=10000 r/min,B0z= 2.0 T,and P0=40 kN/m2)with simply supported sides is plotted shown in Fig.18.Here,ˉω as a dimensionless parameter is the ratio of excitation frequency and natural frequency of undamped linear system.

Fig.17.Brifurcation diagram varying with excitation amplitude P0(a)in both-clamped supported sides and(b)in clamped-outer and simply-inner sides.Here,B0z=2 T,?=10000 r/min,and ω=1.02ω0,where ω0 is natural frequency of undamped linear system.

Fig.18.Bifurcation diagram of the system varying with control parameter ˉω with both-simply supported sides.Here,?=10000 r/min,B0z=2.0 T, and P0=40 kN/m2.

Figure 18 embodies the motions of the system varying with dimensionless parameterˉω.For example,whenˉω=0.2, the motion of the system is single-period motion shown in Figs.19(a)and 19(d),and the energy of the system mainly focuses on the single-frequency.However,whenˉω=0.3 (Fig.20),the motion of the system is chaotic motion.At this point,Poincaré map is an array of lattice with the special geometry(Fig.20(a)),and there are some noise signals in power spectrum(Fig.20(c)),but the power of the system is still centralized at single frequency compared with that in Fig.19(c).

Fig.19.The motion of the system is a single-period motion when=0.2,which is reflected in(a)Poincaré map that obtains one mapping point,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

In order to reflect the effect of the boundary conditions on the system,the bifurcation diagram is plotted with the selection of two kind of boundary conditions and keeping other parameter sun changed,shown in Fig.21,whereˉω is still bifurcation control parameter.The motions of the system are mainly quasi-period motion and multi-period motion,and there are slight chaotic motions only in the 1/3 subharmonics region (=3)For instance,the response diagrams of the system in primary resonance region(=1)and subharmonics resonance region(=3)are plotted in Figs.22 and 23,respectively.When=1,the motion of the system is quasi-period motion,but single period when=3.

Fig.21.Bifurcation diagram varying with control parameter with bothclamped supported sides and other things being the same as Fig.18.

Fig.22.The motion of the system is a quasi-period motion when=1,which is reflected in(a)Poincaré map that is a mapping special shapes,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

Fig.23.The motion of the system is a single-period motion when=3,which is reflected in(a)Poincaré map that is a mapping point,(b)time history diagram,(c)power spectrum,and(d)phase diagram.

5.Conclusion

In this paper,Bessel functions as mode shape functions are utilized to solve the magneto-elasticity nonlinear vibration equation of the annular plate.The bifurcations diagrams,response charts power spectrum charts,phase diagrams,Poincaré maps,and time history diagram of the systems are obtained with different bifurcation control parameters by numerical calculating.When magnetic induction,excitation amplitude,or frequency respectively reaches a specific value,there may be bifurcations and chaos in the system, which are very sensitive to the control parameters.In addition,there are multi-period motions between the chaos fields in the system.In most cases,with the increases of the bifurcation control parameters,the motions of the system are repeatedly from multi-period motion to chaotic motion and to multi-period motion.Moreover,when magnetic induction intensity increases to some degree,the motion of the system may be multi-period motion or quasi-period motion rather chaotic motion.

[1]Leo A,Saeedi K and Rama B 2012 J.Mech.Sci.Technol.26 1439

[2]Malekzadeh P,Haghighi M R G and Atashi M M 2011 Acta Mech.46 893

[3]Saidi A R,Baferani A and Jomehzadeh E 2011 Acta Mech.219 309

[4]Allahverdizadeh A,Naei M H and Bahrami M N 2008 J.Sound Vib. 310 966

[5]Ratko M 2005 J.Sound Vib.280 467

[6]Hashemi S H,Farhadi S and Carra S 2009 J.Sound Vib.323 366

[7]Zheng X J,Zhang J P and Zhou Y H 2005 Int.J.Solid.Struct.42 2417

[8]Gao Y,Xu B and Huh H 2010 Acta Mech.210 99

[9]Hasanyan D J and Librescu L 2005 Comput.Struct.83 1205

[10]Hu Y D and Wang T 2016 Nonlinear Dyn.85 1825

[11]Li X Y,Ding H J and Chen W Q 2008 Compos.Struct.83 381

[12]Alaimo A,Benedetti I and Milazzo 2014 Compos.Struct.107 643

[13]Razavi S and Shooshtari A 2015 Compos.Struct.119 377

[14]Lu Q S,To C W S and Huang K L 1995 J.Sound Vib.181 873

[15]Hu Y D and Zhang Z Q 2011 Chaos Soliton.Fract.44 739

[16]Hu Y D,Hu P and Zhang J Z 2015 J.Comput.Nonlinear Dyn.10 021010

[17]Hu Y D and Zhang Z Q 2012 Nonlinear Dyn.67 1779

[18]Touzé C,Thomas O and Amabili M 2011 Int.J.Non-Linear Mech.46 234

[19]Coman C D 2013 Mech.Res.Commun.47 11

[20]Shahverdi H and KhalafiV 2016 Compos.Struct.146 84

[21]Zhao D M and Zhang Q C 2010 Chin.Phys.B 19 030518

[22]Li S and Zhang J Y 1997 Metal Forming Technol.15 34(in Chinese)

[23]Zhang Y S,Gao D P and Yin L Y 1989 J.Nanjing Aerount.Inst.21 18 (in Chinese)

22 November 2016;revised manuscript

26 April 2017;published online 11 August 2017)

10.1088/1674-1056/26/9/094302

?Project supported by the National Natural Science Foundation of China(GrantNo.11472239),the Hebei Provincial Natural Science Foundation of China(Grant No.A2015203023),and the Key Project of Science and Technology Research of Higher Education of Hebei Province of China(Grant No.ZD20131055).

?Corresponding author.E-mail:huyuda03@163.com

?2017 Chinese Physical Society and IOP Publishing Ltd http://iopscience.iop.org/cpb http://cpb.iphy.ac.cn