当前位置：文档库 › A full-discretization method for prediction of milling stability

A full-discretization method for prediction of milling stability

Short Communication

Ye Ding a ,LiMin Zhu a ,XiaoJian Zhang b ,Han Ding a,n

State Key Laboratory of Mechanical System and Vibration,School of Mechanical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China

b State Key Laboratory of Digital Manufacturing Equipment and Technology,School of Mechanical Science and Engineering,Huazhong University of Science and Technology,Wuhan 430074,China

a r t i c l e i n f o

Article history:

Received 13December 2009Received in revised form 10January 2010

Accepted 12January 2010

Available online 2February 2010Keywords:

Milling stability

Full-discretization method Floquet theory Time delay

a b s t r a c t

This paper presents a full-discretization method based on the direct integration scheme for prediction of milling stability.The fundamental mathematical model of the dynamic milling process considering the regenerative effect is expressed as a linear time periodic system with a single discrete time delay,and the response of the system is calculated via the direct integration scheme with the help of discretizing the time period.Then,the Duhamel term of the response is solved using the full-discretization method.In each small time interval,the involved system state,time-periodic and time delay items are simultaneously approximated by means of linear interpolation.After obtaining the discrete map of the state transition on one time interval,a closed form expression for the transition matrix of the system is constructed.The milling stability is then predicted based on Floquet theory.The effectiveness of the algorithm is demonstrated by using the benchmark examples for one and two degrees of freedom milling models.It is shown that the proposed method has high computational ef?ciency without loss of any numerical precision.The code of the algorithm is also attached in the appendix.

1.Introduction

Machining stability prediction is important for optimal selec-tion of spindle speed and cutting depth to avoid chatter and improve production ef?ciency.Many methods including analy-tical,numerical and experimental ones [1–17]have been proposed.Smith and Tlusty [1]presented a method to generate stability lobes by time domain simulations of the chatter vibrations in milling process.Sridhar et al.[2]developed a mathematical model to describe the dynamic milling process and solved it by a numerical method.Minis and Yanushevsky [3]proposed a comprehensive analytic method and solved the two dimensional milling problem by introducing the theory of periodic differential equations.Altintas and Budak [4]presented an analytical method (ZOA method)for predicting milling stability lobes based on the mean of the Fourier series of the dynamic milling coef?cients.

The ZOA method is ef?cient and fast,but it cannot predict the existence of the additional stability regions and period doubling bifurcations in the case of low radial depth cut for milling.To overcome this problem,several methods have been developed.The multi-frequency (MF)solution of chatter stability was explored by Budak and Altintas [5]and then extended by Merdol

and Altintas [6].It considers harmonics of the tooth spacing angle and spread of the transfer function with the harmonics of the tooth passing frequencies.The temporal ?nite element analysis (TFEA)for milling process simulation was presented by Bayly et al.[7].It is an extension of the method developed by Bayly et al.[8]for simulation of an interrupted turning process.The semi-discretization (SD)method,developed by Insperger and Ste

′pa ′n [9,10]is an ef?cient numerical method for stability analysis of linear delayed systems.It can be applied to predict milling stability.

Also,some experimental methods are utilized to get the stability boundaries in milling.Ismail and Soliman [11]quickly identi?ed the stability lobes in milling by ramping the spindle speed while monitoring the behaviour of a chatter indicator.Solis et al.[12]combined a chatter’s analytical prediction method with the experimental multi-degree of freedom system modal analysis to obtain the stability lobes information.More recently,Quintana et al.[13]determined the stability charts of a milling process by applying sound mapping methodology.

The methods reviewed above have their advantages and disadvantages,respectively:

The ZOA method provides an analytic solution.It is the fastest

approach that solves for the chatter free cutting condition so far.However,it is not very suitable for the low radial immersion problem [14].

Contents lists available at ScienceDirect

journal homepage:https://www.wendangku.net/doc/1114397247.html,/locate/ijmactool

International Journal of Machine Tools &Manufacture

Corresponding author.Tel.:+862134206086;fax:+862134206086.E-mail address:hding@https://www.wendangku.net/doc/1114397247.html, (H.Ding).

International Journal of Machine Tools &Manufacture 50(2010)502–509

The TFEA method is ef?cient and accurate for small times in the cut,but not quite suitable in full and near-full immersion cases[15].

The semi-discretization and multi-frequency methods take into account the effect of higher harmonics mostly due to multiple mode excitation or highly interrupted cutting,which happens when the radial immersion is lower than10%of the cutter diameter[16].However,their computational ef?cien-cies are not high.

Time domain numerical simulation methods are quite power-ful.They take true kinematics of the milling,mechanics of cutting,the in?uence of inner and outer modulation,cutter geometry,runout and other non-linearities into consideration, but their computational costs are too expensive[17].

The experimental methods need expensive apparatus and it is time consuming to conduct experiments to obtain the milling stability lobes information.

The dynamic milling process considering the regenerative effect is generally modeled as a linear time periodic system with a single discrete time delay,which can be approximately solved by using the analytical or numerical methods reviewed above.In[9],the semi-discretization method was introduced since it was found to be more effective than the full-discretization method which discretizes all the actual time domain terms.Actually,from the viewpoint of discretization of a continuous time delay system,the full-discretization method reported in[9]is based on the direct difference scheme.Another kind of full-discretization methods is based on the direct integration scheme,which was suggested by Zhong et al.[18,19].This scheme has more robust numerical stability than the former one.In this paper,a full-discretization method based on the direct integration scheme is proposed for prediction of milling stability.It has high computational ef?ciency. The remainder of this paper is organized as follows.In Section2, the mathematical model is introduced and the numerical schemes are described in detail.In Section3,two benchmark examples for one and two degrees of freedom milling models are given to illustrate the ef?ciency of the novel approach.In Section4, conclusions with a brief discussion on future works are presented.

2.Mathematical model and algorithm

Without loss of generality,the dynamics of the milling process considering the regenerative effect can be described by a n-dimensional linear time periodic system with a single discrete time delay expressed in the following state-space form

_xetT?A

xetTtAetTxetTtBetTxetàTTe1Twhere A0is a constant matrix representing the time-invariant nature of the system,A(t)and B(t)are two periodic-coef?cient matrices satisfying A(t+T)=A(t+T)and B(t+T)=B(t+T),and T is the time period which equals to the time delay.

The?rst step to solve Eq.(1)numerically is to discretize the time period T,i.e.,equally divide T into m small time intervals such that T=m t,where m is an integer.On each time interval k t r t r(k+1)t, (k=0,y,m),the response of Eq.(1)with the initial condition x k=x(k t) can be obtained via the direct integration scheme as follows:

xetT?e A0etàk tTxek tTt

Z t

k t

f e A0etàxT?AexTxexTtBexTxexàTT

g d xe2T

Eq.(2)can be equivalently expressed as

xek tttT?e A0t xek tTtZ t

e A0x

Aek tttàxTxek tttàxTt

Bek tttàxTxek tttàxàTT

()

d x

e3T

where0r t r t.Then,x k+1,i.e.x(k t+t)can be got from Eq.(3):

x kt1?e A0t xek tTt

Z t

e A0x

Aek tttàxTxek tttàxTt

Bek tttàxTxek tttàxàTT

()

d xe4T

The next step is to handle the Duhamel term of Eq.(4),i.e.the

integral item.Firstly,the time delay item x(k t+tàxàT)is

approximated linearly using x k+1àm and x kàm,i.e.the two

boundary values at the time interval[(kàm)t,(k+1àm)t],

resulting in

xek tttàxàTT6x kt1àmtxex kàmàx kt1àmT=te5T

The state item x(k t+tàx)in Eq.(4)can also be approximated

linearly using x k and x k+1,i.e.the two boundary values at the time

interval[k t,(k+1)t],resulting in

xek tttàxT6x kt1txex kàx kt1T=te6T

Similarly,the time-periodic items A(k t+tàx)and B(k t+tàx)

in Eq.(4)could also be approximated by linearly interpolating

the two boundary values at the time interval[k t,(k+1)t],

resulting in

Aek tttàxT6AekT0tAekT1xe7T

where AekT

?A kt1,AekT

?A kàA kt1

àá

=t,and A k denotes the value

of A(t)sampled at the time t k=k t,and

Bek tttàxT6BekT0tBekT1xe8T

where BekT

?B kt1,BekT

?B kàB kt1

àá

=t and B k denotes the value of

B(t)sampled at the time t k=k t.

Substituting Eqs.(6)–(8)into Eq.(4)leads to

x kt1?eF0tF0;1Tx ktF kt1x kt1tF mà1x kt1àmtF m x kàme9T

where

F0?U0e10T

F0;1?eU2=tTAekT0teU3=tTAekT1e11T

F kt1?eU1àU2=tTAekT0teU2àU3=tTAekT1e12T

F mà1?eU1àU2=tTBekT0teU2àU3=tTBekT1e13T

F m?eU2=tTBekT0teU3=tTBekT1e14T

U0?e A0t;U1?

Z t

e A0x d x;U2?

Z t

x e A0x d x;U3?

Z t

x2e A0x d x

e15T

Clearly,U1,U2and U3can be expressed explicitly in terms of

matrices U0and Aà10,i.e.

U1?Aà1

eU0àITe16T

U2?Aà1

et U0àU1Te17T

U3?Aà1

et2U0à2U2Te18T

Alternatively,the matrices U0–U3in Eq.(15)can be numeri-

cally calculated using the precise time-integration(PTI)method

without calculating the inverse matrix Aà1

[20,21].

From Eq.(9),it is obvious that if matrix[IàF k+1]is nonsingular,

x k+1can be expressed in the following explicit form:

x kt1??IàF kt1 à1eF0tF0;1Tx kt?IàF kt1 à1F mà1x kt1àm

t?IàF kt1 à1F m x kàme19TY.Ding et al./International Journal of Machine Tools&Manufacture50(2010)502–509503

In the case that matrix[IàF k+1]is singular,the state item x(k t+tàx)in Eq.(6)could be replaced by the zero-order hold expression,i.e.

xek tttàxT?x k:e20TThen,x k+1becomes

x kt1?eF0tF0;2Tx ktF mà1x kt1àmtF m x kàme21Twhere

F0;2?U1AekT

0tU2AekT

e22T

According to Eq.(19),a discrete map can be de?ned as

y kt1?D k y ke23Twhere the n(m+1)dimensional vector y k denotes

y k?colex k x kà1...x kt1àm x kàmTe24Tand D k is de?ned as

D k?

?IàF kt1 à1eF0tF0;1T00ááá0?IàF kt1 à1F mà1?IàF kt1 à1F m I00ááá000

0I0ááá000

^^^&^^^

000ááá000

000áááI00

000ááá0I0

e25T

Now,the transition matrix U over one periodic time interval can be constructed by using the sequence of discrete maps D k, (k=0,y,mà1),i.e.

y m?U y0e26T

where U is de?ned by

U?D mà1D mà2...D1D0e27TNow,according to Floquet theory[22],the stability of the system can be determined:if the moduluses of all the eigenvalues of the transition matrix U are less than unity,the system is stable, otherwise,unstable.

Remark.The main differences between the proposed method and the well known semi-discretization method[10]lie in three aspects.Firstly,the key point of the semi-discretization method is that only the time-delay term is discretized while the time domain terms are all unchanged.However,in this paper,the state item x(k t+tàx)and time-delay term x(k t+tàxàT)in Eq.(4) are approximated by linearly interpolating the two boundary values at the corresponding time intervals[k t,(k+1)t]and [(kàm)t,(k+1àm)t],respectively,which means that it is a full-discretization method.Secondly,in the semi-discretization method,the periodic-coef?cient matrices of the system are approximated by their averaging values over the small time interval;while,in this paper,the time-periodic terms A(k t+tàx) and B(k t+tàx)in Eqs.(7)and(8)are all approximated by means of linear interpolation of boundary values.Finally and the most importantly,when calculating for the milling stability prediction, the matrices U0,U1,U2and U3in Eq.(15)depend only on the spindle speed,and do not affected by the depth of cut.This implies that during the process of sweeping the range of the depth of cut to determine the transition matrix U,no additional calculation is needed to compute these matrix exponentials repeatedly.However,computation of the exponential of a matrix is necessary while sweeping the range of the depth of cut in the semi-discretization method.Supposing the parameter plane of the spindle speed–depth of cut is divided into a N s?N a sized grid,the matrix exponentials must be calculated N s?N a?m times to obtain the stability chart using the semi-discretization method,while only N s times are needed to calculate the matrix exponentials by using the proposed method.This is the main reason that the computational ef?ciency of the proposed method is higher than that of the semi-discretization method.In the next section,this point will be emphasized in detail via examples.

3.Veri?cation

In this section,two benchmark examples for one and two degrees of freedom(DOF)milling models are illustrated.To compare with the semi-discretization method,the same model parameters from Insperger and Ste′pa′n[10]are https://www.wendangku.net/doc/1114397247.html,puter programs of the proposed approach are all written in MATLAB s 7.4and implemented on a personal computer[Intel Core(TM)2 Duo Processor,2.1GHz,1GB].

3.1.Single DOF milling model

The dynamic equation of a single DOF milling model is [7,10]

€xetTt2z o n_xetTto2

xetT?à

whetT

m t

exetTàxetàTTTe28Twhere z is the relative damping,o n is the angular natural frequency,w is the depth of cut,and m t is the modal mass of the tool.The time delay T is equal to the tool passing period60/(N O), where N is the number of the cutter teeth and O is the spindle speed in rpm.h(t)is the cutting force coef?cient which is de?ned as

hetT?

X N

j?1

gef jetTTsinef jetTT?K t cosef jetTTtK n sinef jetTT e29T

where K t and K n are the tangential and the normal linearized cutting force coef?cients,respectively,and f j(t)is the angular position of the j th tooth de?ned by

etT?e2p O=60Tttejà1T2p=Ne30TThe function g(f j(t))is de?ned as

gef jetTT?

1if f st o f jetTo f ex

0otherwise

e31T

where f st and f ex are the start and exit angles of the j th cutter tooth.

For down-milling,f st=arcos(2a/Dà1)and f ex=p;for up-milling,f st=0and f ex=arcos(1à2a/D),where a/D is the radial immersion ratio.

Let yetT?m t_xetTtm t z o n xetTand xetT?xetTyetT

??T

.Through some simple transformations,the state-space form of the single DOF milling model can be represented as

_xetT?A

xetTtAetTxetTtBetTxetàTTe32Twhere

A0?

àz o n

m tez o nT2àm t o2nàz o n

75;AetT?00

àwhetT0

;

BetT?

whetT0

e33T

The system parameters are:a two?uted cutter,the natural frequency f n=o n/(2p)=922Hz,the relative damping is z=0.011,

Y.Ding et al./International Journal of Machine Tools&Manufacture50(2010)502–509 504

the modal mass is m t=0.03993kg,the cutting force coef?cients are K t=6?108N/m2and K n=2?108N/m2.Since the item_xetàTTdoes not appear in Eq.(28),the dimension of the transition matrix can be reduced[10].

The parameter m is chosen as40,the same as the one employed in the semi-discretization method[10].The stability charts are calculated using both methods over a400?200 sized grid of parameters.The computational time is summarized in Table1for down-milling with radial depth of cut ratios a/D=1,0.1,0.05,respectively.To demonstrate the difference between the presented method and the semi-discretization method more clearly,the MATLAB code of the proposed algorithm for the single DOF milling stability analysis,which follows the same names of variables and the program structure in [10],is given in Appendix A.For this example,the matrix exponentials must be calculated401?201?40=3224040 times by using the semi-discretization method,while only401 times are needed to calculate the matrix exponentials by using the presented method.Since there is no need to compute the matrix exponential in the loop of sweeping the range of the depth of cut,the computational burden is reduced remarkably in contrast to the semi-discretization method;meanwhile,

Table1

Y.Ding et al./International Journal of Machine Tools&Manufacture50(2010)502–509505

the numerical precision is still high https://www.wendangku.net/doc/1114397247.html,pared with the

semi-discretization method,the computation time of the proposed method can be reduced by nearly 75%.3.2.Two DOF milling model

The dynamic equation of a two DOF milling model with a symmetric tool is [7,10]m t 00m t "#€x et T€y et T"#t

2m t z o n

002m t z o n

#_x et T_y et T"#t

m t o 2n

00m t o 2n

"#x et Ty et T"#

?àwh xx et Tàwh xy et Tàwh yx et Tàwh yy et T"#x et Ty et T"#t

wh xx et Twh xy et Twh yx et Twh yy et T"#x et àT Ty et àT T

"#e34T

where z is the relative damping,o n is the angular natural

frequency,and m t is the modal mass of the cutter,and they are assumed to be equal at x and y directions.h xx (t ),h xy (t ),h yx (t )and h yy (t )are the cutting force coef?cients de?ned as h xx et T?

X N j ?1g ef j et TTsin ef j et TT?K t cos ef j et TTtK n sin ef j et TT

e35T

h xy et T?X N j ?1

g ef j et TTcos ef j et TT?K t cos ef j et TTtK n sin ef j et TT

e36T

h yx et T?X N j ?1g ef j et TTsin ef j et TT?àK t sin ef j et TTtK n cos ef j et TT e37T

h yy et T?X N j ?1

g ef j et TTcos ef j et TTàK t sin ef j et TTtK n cos ef j et TT

h i

e38T

All the parameters here have the same meanings as those employed in the single DOF model.

To represent the original system Eq.(34)in the state-space form,let M ,C ,K and q (t )denote the matrices m t 00

m t

"#;

2m t z o n

"#;m t o 2n 00

m t o 2n

"#and

x et Ty et T

in Eq.(34),respectively.Then,let p et T?M _q tCq =2and x (t )denotes q et Tp et T??T

[23].Finally,the two DOF milling model can be represented as

_x et T?A 0x et TtA et Tx et TtB et Tx et àT Te39T

where

A 0?

àM à1C =2

M à1

CM à1C =4àK àCM à1=2

"#;A et T?00000000àwh xx et Tàwh xy et T0

0àwh yx et Tàwh yy et T0

664

77775;Table 2

Y.Ding et al./International Journal of Machine Tools &Manufacture 50(2010)502–509

506

B et T?0000

0000wh xx et T

wh xy et T00wh yx et T

wh yy et T

66664

77775e40T

The system parameters selected are the same as those employed in the single DOF milling model case.The parameter m is still chosen as 40.The stability charts are calculated by using the proposed method and the semi-discretization method over a 400?200sized grid of parameters,and the computational time are listed in Table 2for up-milling with radial depth of cut ratios a /D =0.1,0.05,respectively.Again,it is shown that the computational ef?ciency of the proposed method is much higher without loss of any numerical https://www.wendangku.net/doc/1114397247.html,pared with the semi-discretization method,the computation time of the proposed method can be reduced by about 60%.

4.Conclusions and future works

In this work,a full-discretization algorithm is proposed for prediction of milling stability.The dynamics of the milling process considering the regenerative effect is described by a linear time periodic system with a single discrete time delay,and the response of the system is calculated by a direct integration scheme with the help of discretizing the time period.The algorithm is easy to implement because on each small time interval,the involved system state,time-periodic and time delay items are simultaneously approximated by means of linear interpolation.More importantly,the algorithm is computationally

ef?cient because the involved matrix exponentials are calculated only in the outer loop for sweeping the range of the spindle speed,and not needed to be updated in the inner loop for sweeping the range of the depth of cut.Two benchmark examples for one and two degrees of freedom milling models,which have been validated by experiments [7],are utilized to verify the proposed algorithm.When using the same sized grid for the spindle speed–depth of cut plane and the same sampling period for the tooth passing period,the computational times are reduced by about 60–75%in comparison with those elapsed by the semi-discretization method [10]without loss of any numerical precision.

The present work focuses on the topic of prediction of milling stability using the full-discretization method.Some other topics are worth considering.Two of them are of the most interest.The ?rst one is to test the potential application of the full-discretization method for prediction of surface location error.The second one is to test the feasibility of the full-discretization method to handle more complex dynamic milling systems,e.g.considering spindle speed variations,cutter runout and the process damping,etc .

Acknowledgements

This work was partially supported by the National Key Basic Research Program under grant 2005CB724103,and the Science &Technology Commission of Shanghai Municipality under grant 09QH1401500.

Appendix A.MATLAB code for the single DOF milling stability

analysis

Y.Ding et al./International Journal of Machine Tools &Manufacture 50(2010)502–509

507

Y.Ding et al./International Journal of Machine Tools&Manufacture50(2010)502–509 508

References

[1]S.Smith,J.Tlusty,Ef?cient Simulation Programs for Chatter in Milling,CIRP

Annals42(1)(1993)463–466.

[2]R.Sridhar,R.E.Hohn,G.W.Long,A stability algorithm for the general milling

process,ASME Journal of Engineering for Industry90(2)(1968)330–334. [3]I.Minis,R.Yanushevsky,A new theoretical approach for the prediction of machine

tool chatter in milling,,Journal of Engineering for Industry115(1)(1993)1–8.

[4]Y.Altintas,E.Budak,Analytical prediction of stability lobes in milling,CIRP

Annals—Manufacturing Technology44(1)(1995)357–362.

[5]E.Budak,Y.Altintas,Analytical prediction of chatter stability in milling—Part

I:general formulation,Journal of Dynamic Systems,Measurement and Control,Transactions of the ASME120(1)(1998)22–30.

[6]S.D.Merdol,Y.Altintas,Multi frequency solution of chatter stability for low

immersion milling,,Journal of Manufacturing Science and Engineering, Transactions of the ASME126(3)(2004)459–466.

[7]P.V.Bayly,B.P.Mann,T.L.Schmitz,D.A.Peters,G.Stepan,T.Insperger,Effects

of radial immersion and cutting direction on chatter instability in end-milling,American Society of Mechanical Engineers,13,Manufacturing Engineering Division,MED,2002,pp.351–363.

[8]P.V.,Bayly,J.E.,Halley,B.P.,Mann,M.A.,Davies,,Stability of interrupted

cutting by temporal?nite element analysis,in:Proceedings of the ASME Design Engineering Technical Conference.Vol.6C,2001,pp.2361–2370. [9]T.Insperger,G.Stepan,Semi-discretization method for delayed systems,,

International Journal for Numerical Methods in Engineering55(5)(2002) 503–518.

[10]T.Insperger,G.Ste′pa′n,Updated semi-discretization method for periodic

delay-differential equations with discrete delay,,International Journal for Numerical Methods in Engineering61(1)(2004)117–141.

[11]F.Ismail,E.Soliman,A new method for the identi?cation of stability lobes in

machining,International Journal of Machine Tools and Manufacture37(6) (1997)763–774.[12]E.Solis,C.R.Peres,J.E.Jime′nez,J.R.Alique,J.C.Monje,A new analytical-

experimental method for the identi?cation of stability lobes in high-speed milling,,International Journal of Machine Tools and Manufacture44(15) (2004)1591–1597.

[13]G.Quintana,J.Ciurana,I.Ferrer, C.A.Rodr?′guez,Sound mapping for

identi?cation of stability lobe diagrams in milling processes,International Journal of Machine Tools and Manufacture49(3-4)(2009)203–211. [14]J.Gradisek,M.Kalveram,T.Insperger,K.Weinert,G.Stepan,https://www.wendangku.net/doc/1114397247.html,ekar,I.

Grabec,On stability prediction for milling,International Journal of Machine Tools and Manufacture45(7–8)(2005)769–781.

[15]P.V.Bayly,J.E.Halley,B.P.Mann,M.A.Davies,Stability of interrupted cutting

by temporal?nite element analysis,Journal of Manufacturing Science and Engineering,Transactions of the ASME125(2)(2003)220–225.

[16]D.S.Merdol,Virtual three-axis milling process simulation and optimization,

Ph.D.Thesis,University of British Columbia,Vancouver,2008.

[17]Y.Altintas,M.Weck,Chatter stability of metal cutting and grinding,CIRP

Annals—Manufacturing Technology53(2)(2004)619–642.

[18]W.,Zhong,S.,Tan,and Z.,Wu,The optimal control of time-delay system.in:

Proceedings of the24th Chinese Control Conference,Guangzhou,P.R.China, 2005,pp.394–398.

[19]W.Zhong,Z.Wu,S.Tan,in:Theory and Computation of State-Space Control,

Science Press,Beijing,2007.

[20]S.Tan,W.Zhong,Precise integration method for duhamel terms arising from

non-homogenous dynamic systems,Chinese Journal of Theoretical and Applied Mechanics39(3)(2007)374–381.

[21]W.X.Zhong, F.W.Williams,A precise time step integration method,

Proceedings of the Institution of Mechanical Engineers.Part C:Mechanical Engineering Science208(6)(1994)427–430.

[22]M.Farkas,in:Periodic Motions,Springer-Verlag,New York,1994.

[23]Y.Gu,B.Chen,H.Zhang,Z.Guan,Precise time-integration method with

dimensional expanding for structural dynamic equations,AIAA Journal39

(12)(2001)2394–2399.

Y.Ding et al./International Journal of Machine Tools&Manufacture50(2010)502–509509