A full-discretization method for prediction of milling stability

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

Contents lists available at ScienceDirect

International Journal of Machine Tools & Manufacture journal homepage: www.elsevier.com/locate/ijmactool

Short Communication

A full-discretization method for prediction of milling stability Ye Ding a, LiMin Zhu a, XiaoJian Zhang b, Han Ding a,n a

State Key Laboratory of Mechanical System and Vibration, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China State Key Laboratory of Digital Manufacturing Equipment and Technology, School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan 430074, China b

a r t i c l e in f o

a b s t r a c t

Article history: Received 13 December 2009 Received in revised form 10 January 2010 Accepted 12 January 2010 Available online 2 February 2010

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 efficiency without loss of any numerical precision. The code of the algorithm is also attached in the appendix. Crown Copyright & 2010 Published by Elsevier Ltd. All rights reserved.

Keywords: Milling stability Full-discretization method Floquet theory Time delay

1. Introduction Machining stability prediction is important for optimal selection of spindle speed and cutting depth to avoid chatter and improve production efficiency. Many methods including analytical, 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 coefficients. The ZOA method is efficient 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

n

Corresponding author. Tel.: + 86 21 34206086; fax: + 86 21 34206086. E-mail address: [email protected] (H. Ding).

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 finite 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 semidiscretization (SD) method, developed by Insperger and Ste´pa´n [9,10] is an efficient 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 identified 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].

0890-6955/$ - see front matter Crown Copyright & 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmachtools.2010.01.003

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

 The TFEA method is efficient 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 than 10% of the cutter diameter [16]. However, their computational efficiencies are not high. Time domain numerical simulation methods are quite powerful. They take true kinematics of the milling, mechanics of cutting, the influence 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 efficiency. The remainder of this paper is organized as follows. In Section 2, the mathematical model is introduced and the numerical schemes are described in detail. In Section 3, two benchmark examples for one and two degrees of freedom milling models are given to illustrate the efficiency of the novel approach. In Section 4, 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 _ ¼ A0 xðtÞ þAðtÞxðtÞ þ BðtÞxðtTÞ xðtÞ

ð1Þ

where A0 is a constant matrix representing the time-invariant nature of the system, A(t) and B(t) are two periodic-coefficient 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 first 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=mt, where m is an integer. On each time interval kt rtr(k+1)t, (k=0, y, m), the response of Eq.(1) with the initial condition xk =x(kt) can be obtained via the direct integration scheme as follows: xðtÞ ¼ eA0 ðtktÞ xðktÞ þ

Zt

feA0 ðtxÞ ½AðxÞxðxÞ þBðxÞxðxTÞgdx

ð2Þ

kt

Eq.(2) can be equivalently expressed as " #) Z t( Aðkt þ txÞxðkt þtxÞ þ eA0 x xðkt þ tÞ ¼ eA0 t xðktÞ þ dx Bðkt þtxÞxðkt þtxTÞ 0

503

where 0rt r t. Then, xk + 1, i.e. x(kt + t) can be got from Eq.(3): xk þ 1 ¼ eA0 t xðktÞ þ

" #) Z t( Aðkt þ txÞxðkt þ txÞ þ eA0 x dx Bðkt þ txÞxðkt þ txTÞ 0

ð4Þ

The next step is to handle the Duhamel term of Eq. (4), i.e. the integral item. Firstly, the time delay item x(kt + t  x T) is approximated linearly using xk + 1  m and xk  m, i.e. the two boundary values at the time interval [(k m)t,(k+ 1 m)t], resulting in xðkt þ txTÞ6xk þ 1m þ xðxkm xk þ 1m Þ=t

ð5Þ

The state item x(kt + t  x) in Eq. (4) can also be approximated linearly using xk and xk + 1, i.e. the two boundary values at the time interval [kt,(k+ 1)t], resulting in xðkt þ txÞ6xk þ 1 þ xðxk xk þ 1 Þ=t

ð6Þ

Similarly, the time-periodic items A(kt + t  x) and B(kt + t  x) in Eq. (4) could also be approximated by linearly interpolating the two boundary values at the time interval [kt,(k+ 1)t], resulting in ðkÞ Aðkt þ txÞ6AðkÞ 0 þ A1 x

ð7Þ

  ðkÞ where AðkÞ 0 ¼ Ak þ 1 , A1 ¼ Ak Ak þ 1 =t, and Ak denotes the value of A(t) sampled at the time tk =kt, and ðkÞ Bðkt þ txÞ6BðkÞ 0 þB1 x

ð8Þ

  ðkÞ where BðkÞ 0 ¼ Bk þ 1 , B1 ¼ Bk Bk þ 1 =t and Bk denotes the value of B(t) sampled at the time tk = kt. Substituting Eqs. (6)–(8) into Eq. (4) leads to xk þ 1 ¼ ðF0 þ F0;1 Þxk þFk þ 1 xk þ 1 þ Fm1 xk þ 1m þFm xkm

ð9Þ

where F0 ¼ U0

ð10Þ

ðkÞ F0;1 ¼ ðU2 =tÞAðkÞ 0 þ ðU3 =tÞA1

ð11Þ

ðkÞ Fk þ 1 ¼ ðU1 U2 =tÞAðkÞ 0 þ ðU2 U3 =tÞA1

ð12Þ

ðkÞ Fm1 ¼ ðU1 U2 =tÞBðkÞ 0 þ ðU2 U3 =tÞB1

ð13Þ

ðkÞ Fm ¼ ðU2 =tÞBðkÞ 0 þðU3 =tÞB1

U0 ¼ eA0 t ; U1 ¼

Z t

eA0 x dx; U2 ¼

0

ð14Þ Z t

xeA0 x dx; U3 ¼

0

Z t

x2 eA0 x dx

0

ð15Þ Clearly, U1, U2 and U3 can be expressed explicitly in terms of matrices U0 and A1 0 , i.e.

U1 ¼ A1 0 ðU0 IÞ

ð16Þ

U2 ¼ A1 0 ðtU0 U1 Þ

ð17Þ

2 U3 ¼ A1 0 ðt U0 2U2 Þ

ð18Þ

Alternatively, the matrices U0–U3 in Eq. (15) can be numerically calculated using the precise time-integration (PTI) method without calculating the inverse matrix A1 0 [20,21]. From Eq. (9), it is obvious that if matrix [I Fk + 1] is nonsingular, xk + 1 can be expressed in the following explicit form: xk þ 1 ¼ ½IFk þ 1 1 ðF0 þ F0;1 Þxk þ ½IFk þ 1 1 Fm1 xk þ 1m

ð3Þ

þ ½IFk þ 1 1 Fm xkm

ð19Þ

ARTICLE IN PRESS 504

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

In the case that matrix [I Fk + 1] is singular, the state item x(kt + t  x) in Eq. (6) could be replaced by the zero-order hold expression, i.e. xðkt þ txÞ ¼ xk :

ð20Þ

Then, xk + 1 becomes xk þ 1 ¼ ðF0 þF0;2 Þxk þFm1 xk þ 1m þ Fm xkm

ð21Þ

sized grid, the matrix exponentials must be calculated Ns  Na  m times to obtain the stability chart using the semi-discretization method, while only Ns times are needed to calculate the matrix exponentials by using the proposed method. This is the main reason that the computational efficiency 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.

where ðkÞ F0;2 ¼ U1 AðkÞ 0 þ U 2 A1

ð22Þ

According to Eq. (19), a discrete map can be defined as y k þ 1 ¼ Dk y k

ð23Þ

where the n(m+ 1) dimensional vector yk denotes yk ¼ colðxk xk1 . . .xk þ 1m xkm Þ

ð24Þ

and Dk is defined as 2 6 6 6 6 6 6 Dk ¼ 6 6 6 6 6 4

½IFk þ 1 1 ðF0 þ F0;1 Þ I

0 0

0 0

 

0 0

½IFk þ 1 1 Fm1 0

0

I

0



0

0

^

^

^

&

^

^

0

0

0



0

0

0

0

0



I

0

0

0

0



0

I

3 ½IFk þ 1 1 Fm 7 0 7 7 7 0 7 7 ^ 7 7 7 0 7 7 5 0 0

3. Verification 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 utilized. Computer programs of the proposed approach are all written in MATLABs 7.4 and implemented on a personal computer [Intel Core (TM) 2 Duo Processor, 2.1 GHz, 1 GB]. 3.1. Single DOF milling model The dynamic equation of a single DOF milling model is [7,10] € þ 2zon xðtÞ _ þ o2n xðtÞ ¼  xðtÞ

ð25Þ Now, the transition matrix U over one periodic time interval can be constructed by using the sequence of discrete maps Dk, (k=0, y, m 1), i.e. ym ¼ Uy0

ð26Þ

ð27Þ

Now, 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(kt + t  x) and time-delay term x(kt + t  x  T) in Eq. (4) are approximated by linearly interpolating the two boundary values at the corresponding time intervals [kt,(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-coefficient matrices of the system are approximated by their averaging values over the small time interval; while, in this paper, the time-periodic terms A(kt + t  x) and B(kt + 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, U2 and U3 in 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 Ns  Na

ð28Þ

where z is the relative damping, on is the angular natural frequency, w is the depth of cut, and mt is the modal mass of the tool. The time delay T is equal to the tool passing period 60/(NO), where N is the number of the cutter teeth and O is the spindle speed in rpm. h(t) is the cutting force coefficient which is defined as

where U is defined by

U ¼ Dm1 Dm2 . . .D1 D0

whðtÞ ðxðtÞxðtTÞÞ mt

N X

hðtÞ ¼

gðfj ðtÞÞsinðfj ðtÞÞ½Kt cosðfj ðtÞÞ þKn sinðfj ðtÞÞ

ð29Þ

j¼1

where Kt and Kn are the tangential and the normal linearized cutting force coefficients, respectively, and fj(t) is the angular position of the jth tooth defined by

fj ðtÞ ¼ ð2pO=60Þt þ ðj1Þ 2p=N

ð30Þ

The function g(fj(t)) is defined as  1 if fst o fj ðtÞ o fex gðfj ðtÞÞ ¼ 0 otherwise

ð31Þ

where fst and fex are the start and exit angles of the jth cutter tooth. For down-milling, fst = arcos(2a/D 1) and fex = p; for upmilling, fst = 0 and fex =arcos(1 2a/D), where a/D is the radial immersion ratio.  T _ þ mt zon xðtÞand xðtÞ ¼ xðtÞ yðtÞ . Through Let yðtÞ ¼ mt xðtÞ some simple transformations, the state-space form of the single DOF milling model can be represented as _ ¼ A0 xðtÞ þ AðtÞxðtÞ þBðtÞxðtTÞ xðtÞ where 2 6 A0 ¼ 4 BðtÞ ¼

zon

mt ðzon Þ2 mt o2n " # 0 0 whðtÞ

0

3 1 mt 7 5; zon

ð32Þ

" AðtÞ ¼

0

0

whðtÞ

0

# ;

ð33Þ

The system parameters are: a two fluted cutter, the natural frequency fn = on/(2p)= 922 Hz, the relative damping is z = 0.011,

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

the modal mass is mt = 0.03993 kg, the cutting force coefficients _ are Kt =6  108 N/m2 and Kn =2  108 N/m2. Since the item xðtTÞ does not appear in Eq. (28), the dimension of the transition matrix can be reduced [10]. The parameter m is chosen as 40, the same as the one employed in the semi-discretization method [10]. The stability charts are calculated using both methods over a 400  200 sized grid of parameters. The computational time is summarized in Table 1 for 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 calculated 401  201  40 = 3 224 040 times by using the semi-discretization method, while only 401 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,

Table 1 Comparison of methods for single DOF milling model.

0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

1

Computational

w (m)

0.01 0.009 0.008 0.007 0.006

w (m)

1

1.5 Ω (rpm)

2

The proposed method

0.005 0.004 0.003 0.002 0.001 1

1.5 Ω (rpm)

2

0.005 0.004 0.003 0.002 0.001 0 0.5

2.5 x 104

1

1.5 Ω (rpm)

2

293.1

Semi-discretization method

The proposed method 0.01 0.009

0.008

0.008

0.007

0.007

0.006

0.006

w (m)

w (m)

0.01 0.009

0.005

0.004

0.003

0.003

0.002

0.002

0.001 0 0.5

0.001 0 0.5

1.5 Ω (rpm)

1328.4

2

2.5 x 104

2.5 x 104

0.005

0.004

1

2.5 4 x 10

0.01 0.009 0.008 0.007 0.006

1339.5

time (s)

time (s)

2.5 x 104

Semi-discretization method

Computational

Computational

2

358.0

0 0.5

a/D=0.05

1.5 Ω (rpm)

0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

1393.5

time (s)

a/D=0.1

The proposed method

w (m)

a/D=1

w (m)

Semi-discretization method

505

1

1.5 Ω (rpm)

291.8

2

2.5 x 104

ARTICLE IN PRESS 506

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

the numerical precision is still high enough. Compared with the semi-discretization method, the computation time of the proposed method can be reduced by nearly 75%.

hyx ðtÞ ¼

3.2. Two DOF milling model

hyy ðtÞ ¼

"

h i gðfj ðtÞÞcosðfj ðtÞÞ Kt sinðfj ðtÞÞ þ Kn cosðfj ðtÞÞ

ð38Þ

gðfj ðtÞÞsinðfj ðtÞÞ½Kt cosðfj ðtÞÞ þ Kn sinðfj ðtÞÞ

0

0

mt

;

2mt zon

0

0

2mt zon

" A0 ¼

M1 C=2

ð36Þ

w (m)

w (m)

a/D=0.1

1

Computational

Computational time (s)

2

2.5 x 104

1

1.5 Ω (rpm)

2

954.1

Semi-discretization method

The proposed method

0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

w (m)

w (m)

1.5 Ω (rpm)

2309.8

time (s)

a/D=0.05

The proposed method 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

1

1.5 Ω (rpm)

2318.2

2

2.5 x 104

0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

M1

CM1 C=4K CM1 =2 2 0 0 0 6 0 0 0 6 AðtÞ ¼ 6 6 whxx ðtÞ whxy ðtÞ 0 4 whyx ðtÞ whyy ðtÞ 0

Table 2 Comparison of methods for two DOF milling model.

0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0.5

;

mt o2n 0

0 mt o2n

#

" and

xðtÞ

#

yðtÞ

ð39Þ

where

j¼1

Semi-discretization method

# "

_ ¼ A0 xðtÞ þ AðtÞxðtÞ þBðtÞxðtTÞ xðtÞ

ð35Þ

gðfj ðtÞÞcosðfj ðtÞÞ½Kt cosðfj ðtÞÞ þ Kn sinðfj ðtÞÞ

# "

mt

in Eq. (34), respectively. Then, let pðtÞ ¼ Mq_ þ Cq=2 and x(t)  T denotes qðtÞ pðtÞ [23]. Finally, the two DOF milling model can be represented as

j¼1

hxy ðtÞ ¼

N X

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

where z is the relative damping, on is the angular natural frequency, and mt is the modal mass of the cutter, and they are assumed to be equal at x and y directions. hxx(t), hxy(t), hyx(t) and hyy(t) are the cutting force coefficients defined as

N X

ð37Þ

j¼1

ð34Þ

hxx ðtÞ ¼

gðfj ðtÞÞsinðfj ðtÞÞ½Kt sinðfj ðtÞÞ þ Kn cosðfj ðtÞÞ

j¼1

The dynamic equation of a two DOF milling model with a symmetric tool is [7,10] " #" #" # # " _ € xðtÞ xðtÞ 0 mt 0 2mt zon þ _ € yðtÞ yðtÞ 0 mt 0 2mt zon " #" # 2 xðtÞ 0 mt on þ yðtÞ 0 mt o2n " #" # " #" # whxx ðtÞ whxy ðtÞ whxx ðtÞ whxy ðtÞ xðtÞ xðtTÞ ¼ þ whyx ðtÞ whyy ðtÞ whyx ðtÞ whyy ðtÞ yðtÞ yðtTÞ

N X

N X

1

1.5 Ω (rpm)

960.6

2

2.5 x 104

2.5 x 104

# ; 0

3

07 7 7; 07 5 0

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

2

0 6 0 6 BðtÞ ¼ 6 6 whxx ðtÞ 4 whyx ðtÞ

3

0

0

0

0

0

whxy ðtÞ

0

whyy ðtÞ

0

07 7 7 07 5 0

ð40Þ

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  200 sized grid of parameters, and the computational time are listed in Table 2 for up-milling with radial depth of cut ratios a/D=0.1, 0.05, respectively. Again, it is shown that the computational efficiency of the proposed method is much higher without loss of any numerical precision. Compared 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

507

efficient 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 semidiscretization 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 first one is to test the potential application of the fulldiscretization 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

ARTICLE IN PRESS 508

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

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

References [1] S. Smith, J. Tlusty, Efficient Simulation Programs for Chatter in Milling, CIRP Annals 42 (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 Industry 90 (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 Industry 115 (1) (1993) 1–8. [4] Y. Altintas, E. Budak, Analytical prediction of stability lobes in milling, CIRP Annals—Manufacturing Technology 44 (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 ASME 120 (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 ASME 126 (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 endmilling, 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 finite 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 Engineering 55 (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 Engineering 61 (1) (2004) 117–141. [11] F. Ismail, E. Soliman, A new method for the identification of stability lobes in machining, International Journal of Machine Tools and Manufacture 37 (6) (1997) 763–774.

509

[12] E. Solis, C.R. Peres, J.E. Jime´nez, J.R. Alique, J.C. Monje, A new analyticalexperimental method for the identification of stability lobes in high-speed milling,, International Journal of Machine Tools and Manufacture 44 (15) (2004) 1591–1597. [13] G. Quintana, J. Ciurana, I. Ferrer, C.A. Rodrı´guez, Sound mapping for identification of stability lobe diagrams in milling processes, International Journal of Machine Tools and Manufacture 49 (3-4) (2009) 203–211. [14] J. Gradisek, M. Kalveram, T. Insperger, K. Weinert, G. Stepan, E. Govekar, I. Grabec, On stability prediction for milling, International Journal of Machine Tools and Manufacture 45 (7–8) (2005) 769–781. [15] P.V. Bayly, J.E. Halley, B.P. Mann, M.A. Davies, Stability of interrupted cutting by temporal finite element analysis, Journal of Manufacturing Science and Engineering, Transactions of the ASME 125 (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 Technology 53 (2) (2004) 619–642. [18] W., Zhong, S., Tan, and Z., Wu, The optimal control of time-delay system. in: Proceedings of the 24th 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 Mechanics 39 (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 Science 208 (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 Journal 39 (12) (2001) 2394–2399.