Bayesian uncertainty quantification and propagation for prediction of milling stability lobe

Mechanical Systems and Signal Processing 138 (2020) 106532

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Bayesian uncertainty quantification and propagation for prediction of milling stability lobe Kai Li a, Songping He a,⇑, Hongqi Liu b, Xinyong Mao b, Bin Li a,b, Bo Luo c a

State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan 430074, China National NC System Engineering Research Center, Huazhong University of Science and Technology, Wuhan 430074, China c Advanced Manufacturing Research Center, University of Sheffield, Western Bank, Sheffield S10 2TN, England, United Kingdom b

a r t i c l e

i n f o

Article history: Received 20 June 2019 Received in revised form 31 October 2019 Accepted 18 November 2019

Keywords: Milling stability Uncertainty Chatter prediction Bayesian inference

a b s t r a c t To predict the stability boundary of cutting operation, it is necessary to establish process force model and dynamic mechanical behavior model. Uncertainty of model parameters may lead to significant differences between predicted and measured stability boundary. In this paper, a novel method for determining boundary of stability in milling, which performs the stability analysis on an uncertain dynamic milling model, is presented. Firstly, the input parameters of milling stability prediction model are considered as random variables. Subsequently, sensitivity analysis is used to determine the most influential parameters, the posterior distributions of all influential parameters are estimated using a Bayesian inference framework, which is built based on an improved Ensemble Markov Chain Monte Carlo (IGWMCMC) algorithm. Compared with the traditional method, IGWMCMC is not affected by affine transformations of space, which can achieve much better convergence on badly scaled problems and the computational costs are cheaper. Finally, the posterior distributions of uncertainties calculated using Bayesian inference are propagated through the computational model. The experimental results verify the reliability of the calculated boundaries of stability lobe diagrams. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Chatter is a self-excited vibration, which is accompanied by instability, chaotic behavior and abnormal fluctuation. It is a major problem that causes low material removal rate, increased tool wear, machine tool failure, poor surface finish, excessive noise and increased cost in machining applications [1,2]. To avoid chatter in the production process, the stability lobe diagram (SLD) is a very effective tool, which shows a distinct boundary that separates unstable and stable machining regimes, to facilitate determination of the optimal cutting parameters in advance. Conventional chatter models typically assume that parameters are constant. However, system dynamics, including natural frequencies, damping ratios and cutting coefficients, are the main parameters that can be changed in the actual machining [3–6,9]. The predicted results often do not match the experimental ones during the cutting tests. Thus, these uncertainties must be taken into account for a robust chatter prediction.

⇑ Corresponding author at: State Key Laboratory of Digital Manufacturing Equipment and Technology School of Mechanical Science and Engineering, Huazhong University of Science and Technology (HUST), 1037 Luoyu Road, Hongshan District, Wuhan 430074, Hubei Province, China. E-mail address: [email protected] (S. He). https://doi.org/10.1016/j.ymssp.2019.106532 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

In recent years, many researchers have chosen different methods to deal with the uncertainty analysis of SLDs. A robust chatter prediction method (RCPM) was proposed by G. Totis [7] to analyze the stability of an uncertain dynamic milling model. It emerged that the RCPM approach is more accurate and reliable than the quick chatter prediction method (QCPM). Hajdu et al [8] combined the extended multi-frequency solution and the method of l-analysis to determine robust stability boundaries for milling operations with uncertain dynamics. Probability-based reliability has been intensively studied and widely applied to describe the uncertainty. Random structural system reliability analysis can provide a quantitative index to represent the influence of uncertainty. In Ref. [9], reliability analysis of a dynamic structural system was introduced into structural analysis of a turning process system. Chatter reliability was defined to represent the probability of stability (no chatter occurs). The Monte Carlo simulation (MCS) method and the first-order second moment (FOSM) method were adopted to predict the reliability of the stability in turning [10,11]. Loukil et al. [12] employed Monte Carlo and other reliability methods to study the variable impact on system instability risk. Unfortunately, it is very difficult for the above probability method to quantify the distributions of input parameters and probabilistic analysis techniques because nondeterministic properties cannot always be exactly represented using probabilistic concepts. In addition, sufficient data is often needed to obtain accurate results. Graham et al. [3,13] developed a robust chatter stability model utilizing the Edge theorem and the Zero Exclusion condition. Uncertainty was introduced into the natural frequency and cutting coefficient, which was determined by experiments. The limitation was that it did not fully evaluate the effects of each model parameter, and several assumptions were necessary, such as small uncertainty or time-invariant systems. Fuzzy algorithm as a non-model intelligent algorithms have been used for various problems in structural dynamics [14,15]. Fuzzy arithmetic techniques have been used for the chatter stability problem in Ref. [16]. The results showed that the fuzzy algorithm can be used to solve the process design problem with robustness to uncertain parameters. Hamann et al. [17] employed fuzzy arithmetic to analyze and investigate the influence of uncertainties on time-delayed stability charts and quantify the effect of given uncertain parameters. However, fuzzy arithmetic involves performing mathematical operations on fuzzy numbers, the calculation procedure is still limited to few parameters due to the high numerical complexity. Zhang et al. [18] developed a robust spindle speed optimization formulation and sensitivity analysis was used to obtain the upper and lower bounds and mean value of SLDs. However, this method calculated SLDs based on the modal parameters of the tool tip under idle conditions, so the upper and lower bounds obtained are relatively large, which still has limitations in choosing cutting parameters under machining conditions. Duncan et al. [19] used Monte Carlo simulation techniques to determine the associated uncertainties in the predicted stability limit at each spindle speed. However, this method not only converges slowly in evaluating the uncertainty interval, but also takes the tool tip FRFs as input parameters under idle condition; so the upper and lower bound intervals obtained are too large to provide useful advice for the parameter selection in an actual milling process. The goal of this paper is to introduce the ideas of Bayesian Inference into cutting stability analysis of a milling process system. The Zero-order approximation (ZOA) method proposed by Altintas and Budak [20] is a very common method for calculating SLDs in the frequency domain, ZOA is used as the computational model in our work, the input parameters of the ZOA model include radial cutting coefficients; tangential cutting coefficients; natural frequency and damping ratio of the tool tip in the x and y directions. In this paper, a novel perspective of uncertainty quantification based on Bayesian inference for SLDs is proposed, which has not been reported in other literatures. In addition to the uncertainties caused by the above parameters, there are some epistemic uncertainties caused by simplified assumptions in the model [21,22]. To address the deviation caused by this assumption, a model bias term,r, is introduced. This can be viewed as a model error in predicting blim and xc for the given input parameters. Bayesian inference can be used to quantify the uncertainty of the influential model parameters (along with the hyper-parameter r). Due to different input parameters having different degree effects on the uncertainty of model outputs, an advanced Morris sensitivity analysis is applied to determine the most influential input parameters. The parameter estimation problem in this study is to infer the distribution of the most influential parameters, which can be achieved by an improved Ensemble Markov Chain Monte Carlo (IGWMCMC) algorithm. Compared with the traditional MCMC and GWMCMC methods [23], IGWMCMC is not affected by affine transformations of space, which can achieve much better convergence on badly scaled problems, and the computationally expensive inner loop can be run in parallel, thus the computational costs are cheaper. A series of experiments with different spindle speeds and axial depths of cut have proved the effectiveness of the proposed method. The whole flowchart of the proposed method is summarized as shown in Fig. 1. This paper is organized as follows. The theoretical derivation of the analytical chatter stability model is first presented in Section 2. Advanced Morris Sensitivity Analysis Method is outlined in Section 3. The Bayesian inference and an improved GWMCMC algorithm can be seen in Section 4. In Section 5, the experimental analysis is carried out. The experimental results are discussed and the proposed method is verified in Section 6. Finally, the work is summarized in the conclusions section.

2. Dynamics milling model The chatter prediction method is based on the joint analysis of the dynamics of the machining system, the external force characteristics and the dynamic interaction between tool and workpiece. The pseudo Single-Degree-of-Freedom (SDOF) chatter stability model of the Multi-Degree-of-Freedom (MDOF) system proposed by Tlusty includes many basic assumptions [24], which often results in large stability prediction errors. By solving

3

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

Amplitude (N)

d=0.2 mm

1.Experiments of different spindle speeds and axial depths of cut

2.The measured surface topography

Amplitude (N)

Frequency (Hz)

d=0.2*n mm

2.The measured surface topography

Chatter

Frequency (Hz)

x1 ⎫ x2 ⎪⎪ ⎪ x3 ⎬ ⎪ ⎪ xn ⎪⎭

Sensitivity analysis(advanced Morris method)

All input model parameters

⎧x1 Prior ⎪x Prior ⎪⎪ 2 ⎨x3 Prior ⎪x Prior ⎪4 ⎪⎩x5 Prior

3.Uncertainty quantification

x2

x3

x4

4.Uncertainty Propagation The zero-order approximation model for stability lobe diagram prediction

x5 Posterior distribution of influential parameters

blim ( mm)

The most influential parameters

x1

× ×

×

× ×

× × ×

×

×

Ω (rpm)

Fig. 1. Uncertainty quantification and propagation framework for the SLDs computational model.

the eigenvalue problem proposed by Altintas and Budak [20,25,26], chatter stability can be presented more accurately. Considering the general equation of motion, a linear analysis model of chatter stability can be derived. An orthogonal milling scheme is shown in Fig. 2. The tool has N t teeth and rotates at a constant angular speed of X rad/s. The general equations of motion for the 2DOF milling system are:

mx €x þ cx x_ þ kx x ¼ € þ cy y_ þ ky y ¼ my y

PN

i¼1 F x;i

PN

i¼1 F y;i

ð1Þ

where mx ; my ; kx ; ky ; cx ; and cy represent the mass, stiffness and damping coefficients in the x and y directions, respectively. Radial F ri and tangential cutting forces F ti excite the structure on tooth i, resulting in dynamic displacement h of the cutter. This variation of chip thickness leads to self-excited vibration or chatter. The dynamic chip thickness is given by the following formula:

hð/i Þ ¼ ½Dxsinð/i Þ þ Dycosð/i Þgð/i Þ

ð2Þ

where Dx ¼ xðtÞ  xðt  TÞ and Dy ¼ yðtÞ  yðt  TÞ represent the dynamic displacements of the cutter structure at the present and previous tooth path, respectively. T is cutter rotation period. The immersion angle varies with time as /i ðtÞ ¼ Xt. The function g is a unit step function that has a value of unity when tooth i is engaged in the workpiece:



g ð/i ðtÞÞ ¼

1

if /st < /i ðtÞ < /ex

0

if /st > /i ðtÞ or /i ðtÞ > /ex

ð3Þ

A chip is produced only when the immersion angle, /st < /i ðtÞ < /ex . /st and /ex represent the entry and exit angles, respectively. Thus,g ð/i ðtÞÞ ¼ 1,hð/i Þ ¼ Dxsinð/i Þ þ Dycosð/i Þ. The tangential (F t;i ) and radial (F r;i ) cutting forces acting on the tooth i are proportional to the axial depth of cut (b) and chip thickness h (dynamic displacement h of cutter) and are defined as:

F t;i ¼ K tc bh þ bK te ¼ K tc bhð/i Þ F r;i ¼ K rc bh þ bK re ¼ K r F t;i

ð4Þ

4

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

Fig. 2. Dynamic model of milling with two degrees of freedom.

where K tc and K rc are the tangential cutting coefficient, and the radial cutting coefficient, respectively, K r ¼ K rc =K tc is the ratio of the coefficients. K te and K re represent tangential and radial edge coefficients, respectively, which have no effect on chatter and are dropped from Eq. (4), leaving only the dynamic part of the cutting force. By analyzing the cutting force of the given material and cutter, the cutting coefficients can be identified by experiment. According to general milling theory in [6], cutting forces in x and y directions can be decomposed into:

Fx ¼ Fy ¼

PN i¼1

PN

 F t;i cosð/i Þ  F r;i sinð/i Þ

ð5Þ

i¼1 F t;i sinð/i Þ  F r;i cosð/i Þ

The expressions of chip thickness Eq. (2) and tooth force Eq. (4) are substituted into Eq. (5) and rearranged in matrix form, yielding



Fx Fy



¼

 axx 1 bK tc 2 ayx

axy ayy



Dx Dy



ð6Þ

The relative displacement between tool and workpiece depends on four periodic milling force coefficients, which map chip thickness to cutting force. After taking the average component of the Fourier series expansion, the integrated functions of dynamic milling force coefficients in the time-varying direction are given by

axx ¼ 12 ½cos2/  2K r þ K r sin2///exst axy ¼ 12 ½sin2/  2/ þ K r cos2///exst

ð7Þ

ayx ¼ 12 ½sin2/ þ 2/ þ K r cos2///exst ayy ¼ 12 ½cos2/  2K r  K r sin2///exst

The system dynamics in Eq. (1) can be expressed in matrix form in the Laplace domain. The direct transfer function matrix at the tool-workpiece contact zone is given by

½GðixÞ ¼



Gxx ðixÞ Gxy ðixÞ Gyx ðixÞ Gyy ðixÞ



ð8Þ

Because only two orthogonal degrees of freedom are considered, this means that the cross-transfer function zero (Gxy ¼ Gyx ¼ 0).Gxx and Gyy are the tool point Frequency Response Functions (FRFs) in x and y directions respectively, which can be expressed as functions of modal parameters. This paper uses the general formula of complex eigenvalues, which is also used in [27]. k X Q r fwgi fwgTi Q r fwgi fwgT i Gxx;yy ðxÞ ¼ þ jx  ki jx  ki i¼1

!

¼

k  X i¼1

Ai Ai þ jx  ki jx  ki

 ð9Þ

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ki ; ki ¼ ni xi  j 1  n2i xi

5

ð10Þ

where k is the number of modes; fwgi is modal vector, Q r is proportional conversion factor,  denotes conjugation. ni is the damping ratio, xi is the natural frequency of the ith mode and x is the excitation frequency. Ai and Ai are residues that depend on the mode shape and can be assumed to be unaffected by the operating conditions [2]. Thus, they can be measured directly under the idle tool condition. The characteristic equation of the system can then be written as

a0 K2 þ a1 K þ 1 ¼ 0

ð11Þ

where the coefficients a0 and a1 are



a0 ¼ Gxx Gyy axx ayy  axy ayx a1 ¼ axx Gxx þ ayy Gyy

ð12Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  K ¼ ðKR þ KI iÞ ¼  2a10 a1  a21  4a0

The complex-valued eigenvalue K is related to the limiting depth of cut blim , and the final expression for chatter-free axial depth and chatter frequency can be determined as follows:

blim ¼ 

2 p KR 1 þ j2 xc T ¼ cos1 NK tc





j2  1 K j¼ I j2 þ 1 KR

ð13Þ

The corresponding spindle speeds can be found as 60x X ¼ Nt ðp2wþ2k pÞ

w ¼ tan1 j

ð14Þ

From the above theoretical derivation, the axial limit depths of cut are related to the radial coefficient K rc , tangential coefficient K tc and FRFs of the tool tip in the x and y directions. However, since the residue can be measured directly, the axial limit depths of cut thus only depend on the radial coefficient K rc , tangential coefficient K tc , and natural frequencies xx , xy and damping ratios nx , ny in the x and y directions. 3. Selection of parameters using sensitivity analysis (SA) Generally, in the computational model, different input parameters have different effects on the uncertainty of the model output. In addition, the higher the parameter dimension is, the higher the calculation cost of uncertainty analysis is. The method developed in this section mainly concerns the factors that have less influence on calculation. Sensitivity analysis first requires a screening method to select non-important factors that can remain unchanged. Morris method is one of the most popular screening methods [28]. The Morris method is generally used in its original form based on the random sampling of so-called trajectories. In this work, an advanced Morris sensitivity analysis is used to isolate the most influential parameters in the model, which use radial points instead of trajectories, and a Latin hypercube sampling instead of a random sampling, improves the efficiency of the Morris method. What’s more, this algorithm reduces the risk to underestimate and fix nonnegligible factors. SA focuses on the model parameters. In other words, the hyper parameter r is not considered because the limit range of the r value is given; it can conceal the sensitivity of the model parameters. We consider r in the uncertainty quantification problem along with other influential model parameters, and thus only the influence of model parameters on output uncertainty is considered. We can calculate the sensitivities of the limit depths of cut blim to the input parameters: the radial coefficient K rc , tangential coefficient K tc , natural frequencies xx and damping ratios nx in the x direction of the tool tip and the natural frequencies xy and damping ratios ny in the y direction of the tool tip. Fig. 3 represents the elementary effects obtained for each input parameter. The larger the variance is, the more significant the effect of the input parameter is on the output. It can be seen from Fig. 3 that the influence of the left side of the black dashed line on the output results is small, that is, the influence of cutting coefficients on milling stability is very small. In other words, the influence of cutting force coefficients on milling stability uncertainty is significantly less than that of the tool tip modal parameters. Therefore, it can be considered that the natural frequencies xx and damping ratios nx in the tool tip in the x direction and the natural frequencies xy and damping ratios ny in the tool tip in the y direction are the input parameters that have significant effects on the SLDs. 4. Uncertainty quantification using Bayesian inference The most influential model parameters along with the hyper parameter r are presented using a vector, h ¼ fxx ; nx ; xy ; ny ; rg. Let yðxi Þ represent the limit depths of cut blim and chatter frequency xc obtained from experiments; yðxi ; hÞ represent the limit depths of cut blim and chatter frequency xc the predictions from the model. They can be related by the following simple equation:

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

Limit between the key influential factors

6

Elementary effects

8

6

4

2

0

Kt

Kr

c

x y x Factors ordered by ascending maximum

cy

Fig. 3. Limit between factors of low and high influence.

yðxi Þ ¼ yðxi ; hÞ þ e þ rðxi Þ; i ¼ 1; 2    n

ð15Þ

where xi is the experimental condition that can be controlled; e is the experimental error term, r represents the model discrepancy (or model bias) term. 4.1. Bayesian method Bayesian inference is a statistical method, it regards parameter estimation as an inference problem, in which probability is used to measure the relative rationality of the results given a model of the system and measured data is applied to dynamic parameter identification [29]. Compared with non-Bayesian methods, one of the main advantages of Bayesian identification methods is that, in addition to the most probable estimates, they also provide a rigorous quantitative assessment of the uncertainties associated with model parameters for the given available data. When new data are obtained, Bayesian rules are used to update the probability estimates of the hypothesis [30]. Let h be the true set of uncertain parameters; suppose we get a set of observation data points xi ,i ¼ 1; 2;    ; n, which is called evidence D. The aim of Bayesian model updating technique is to find all possible parameters of the system that are consistent with the measured data. By following the Bayesian theorem, the updated probability distributions (posterior distributions),pðhjDÞ, can be obtained

pðhjDÞ ¼

pðDjhÞpðhÞ pðDÞ

ð16Þ

where pðhÞ is the prior Probability Density Function (PDF) that represents the available knowledge or user’s judgment about h before the evidence D is observed. The likelihood function pðDjhÞ denotes the probability of observing the evidence D assuming that the true parameters of such a model are h. 1=pðDÞ is the normalizing constant such that the integration of the posterior PDF over the parameter space equals unity. It is given by

Z

pðDÞ ¼

pðDjhÞpðhÞdh

ð17Þ

Then, posterior and likelihood distributions are as follows:

pðhjDÞ / pðDjhÞpðhÞ

ð18Þ

The prior distributions pðhÞ of all parameters will be discussed in detail in Section 5. Statistical models are used as shown in Eq. (15), where the model error r is assumed to be independent and identically distributed [20]. The likelihood pðDjhÞ of observing the data follows a normal distribution [31–33]:

# n 1 X 2 pðyjh; rÞ ¼ exp  2 ðyðxi Þ  yðxi ; hÞÞ ; i ¼ 1; 2    n n=2 2r i¼1 ð2pr2 Þ 1

"

ð19Þ

where r is a hyper parameter in the likelihood distribution, which can be inferred using the limit depths of cut and chatter frequency data fy1 ; y2 ;    yn g collected independently from n different cutting tests.

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

7

4.2. An improved ensemble Markov chain Monte Carlo method for evaluating posterior distributions of model parameters MCMC uses posterior density attributes to specify parameter values to fully explore the geometry of the distribution [34]. In this work, an improved ensemble Markov chain Monte Carlo method, which is called IGWMCMC, for evaluating posterior densities of model parameters is proposed. Compared with the traditional GWMCMC method [23], the computationally expensive inner loop can be run in parallel, thus the computational costs are lower. Let us consider an ensemble of walkers X ¼ ðx1 ; x2 ;    ; xL Þ and a Markov chain that walks on a product space, the distribution is

PðXÞ ¼ Pðx1 ;    ; xL Þ ¼ pðx1 Þpðx2 Þ    pðxL Þ

ð20Þ

Starting with Xð1Þ, a sequence XðtÞ is produced. The ensemble Markov chain can preserve the product density Eq. (20) without the individual walker sequences X k ðtÞ being independent. This method involves evolving an ensemble of K walkers W ¼ fX k g simultaneously where the proposal distribution for one walker K  1 is based on the current positions of the walkers in the complementary ensemble W k ¼ fX i ; 8i–kg. In [23], to update the position of a walker at position X k , the simplest move of this type is the stretch move. In a stretch move, a walker X i is drawn randomly from the remaining walkers W k . However, this violates the detailed balance in some cases and has low computing efficiency [44]. As a result, we divide the ensemble into two subsets W ð0Þ ¼ fX k ; 8k ¼ 1;    ; K=2g and W ð1Þ ¼ fX k ; 8k ¼ K=2 þ 1;    ; K g. Update all the walkers in W ð0Þ simultaneously based only on the positions of the walkers in the other set W ð1Þ . Then, using the new positions W ð0Þ , we can update W ð1Þ . In this case, the outcome is a valid step for all of the walkers. An improved algorithm is provided in the Appendix. A new position is obtained:

X k ðtÞ ! Y ¼ X i þ V ½X k ðtÞ  X i 

ð21Þ

where V is a random variable drawn from a distribution gðV ¼ vÞ. Obviously, if g satisfies

g



v 1



¼ v gðv Þ

ð22Þ

The particular distribution we use is the one suggested in [35]:

( gðv Þ /

1

 ;a

p1ffiffiffi

if v 2

0

otherwise

v

a

ð23Þ

where a is an adjustable scale parameter. Then, Eq. (21) is symmetric. In this case, the chain will satisfy the detailed balance if the proposal is accepted with probability

( q ¼ min 1;

k Y  X i kn1 PðYÞ

)

k X k ðtÞ  X i kn1 PðX k ðtÞÞ

  PðYÞ ¼ min 1; V n1 PðX k ðtÞÞ

ð24Þ

where n is the dimension of the parameter space. Then, this process is repeated for each walker in the ensemble following the procedure shown in Appendix. The convergence criteria for the posterior distributions of all model parameters are measured by autocorrelation [36]. Autocorrelation qxx ðtÞ for a sequence is the correlation between two points separated by a fixed distance t,

h  i   xi  x xiþt  x  qxx ðtÞ ¼   2 E xi  x E

ð25Þ

where EðÞ represents the expected value. MCMC will converge if the autocorrelation of each parameter gradually approaches zero during sampling. 4.3. Uncertainty propagation The posterior distributions of uncertainties calculated by Bayesian inference are propagated through the model. In this work, a Monte Carlo-based sampling technique is used because it is simple to implement and its efficiency is independent of the number of parameters in the model [22]. In each Monte Carlo iteration, the set of the most influential model parameters are sampled from their respective posterior distributions, as shown in Section 4.2. The schematic of the uncertainty propagation methodology is shown in Fig. 4. In every Monte Carlo iteration each of the influential parameters is assigned a value (depicted using a solid circle on the left-hand side of Fig. 4) by sampling from their posterior distributions and fed into the ZOA model along with the chatter frequency (xc ) and axial limit depths of cut (blim ). For a set of influential parameters, each Monte Carlo iteration generates a SLD. In our work, the number of iterations is 50, which proves to be reliable.

8

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

5. Experimental analysis 5.1. Test methodology All of the experiment tests were conducted on a 5-axis Computer Numerical Control (CNC) machining center to verify the proposed methodology, which was built by the National Engineering Research Center of Numerical Control System of China. The system dynamics under the idle condition were determined by experimental modal analysis (EMA), and the cutting coefficients were identified for a given workpiece. The workpiece material was a block of Al6061 aluminum of size 150 mm  100 mm  50 mm, which was clamped on the worktable; additional parameters of the workpiece are listed in Table 1. The cutter that was used in the milling experiment is an 8 mm diameter end mill; detailed geometric parameters are listed in Table 2. Then the chatter tests were carried out to simulate the chatter conditions and compared with the reliable predicted stability lobes. In milling chatter tests, a dynamometer (Kistler 9129A) was used to capture the milling forces with a sampling rate of 10 kHz. Macro-scale images of the machined surface were obtained for each chatter tests using a portable microscope.

5.2. Identification of tool tip FRF under the idle condition The model parameters as input variables to the milling stability prediction model are the structural dynamics parameters, i.e., the nature frequency and damping ratio of the tool tip, and the cutting coefficients. The FRFs of the tool tip under the idle condition were obtained by impact hammer test. The dynamic characteristics of the tool tip were acquired with an impact hammer (Dytran, sensitivity of 9.7 mV/N). A low-mass (1.0 g), wide-bandwidth accelerometer was mounted at the tool tip to monitor the vibrations in the x in the y direction. The experimental test bench and the measured FRFs of tool tip are shown in Fig. 5. Through the impact hammer tests, the dominant modal parameters for each set of data are identified by plotting the FRFs and fitting the data curve with Eq. (9). The actual measurements and fitting results are shown in Fig. 6a–d. Figs. 6e and 6f show the coherence functions of the input force and the vibration of the tool tip in the x direction and in the y direction, respectively. It can be seen from the figure that the measurements are in good agreement with the fitting results. To obtain a reasonable uncertainty interval, repeated modal impact tests are usually required. Measured FRFs obtained from 10 different measurements are presented in Fig. 5b. Each measurement is indicated by solid blue line, while the upper and lower

Fig. 4. Uncertainty propagation using Monte Carlo sampling.

9

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532 Table 1 Material properties of Al6061. Material type

Tensile strength (MPa)

Yield strength (MPa)

Fatigue strength (MPa)

Modulus of elasticity (Gpa)

Percentage elongation (%)

Al6061

124

55.2

62.1

68.9

25

Table 2 Tool parameters. Terminology

Values

Helix angle (degree) Tool diameter (mm) Tool blade length (mm) Number of teeth Tool length L1 (mm) Tool overhang

30 8 25 4 65 25

bounds are denoted by black lines (±1r). A 10-degree-of-freedom system is fitted onto the average of the measurements. As the modal parameters of each mode within the measurement bandwidth are identified, the FRF of each individual mode can be obtained, and they will be simply added to fit the direct FRF by frequency basis. It can be seen that neither the measurement nor the fitting process is perfect, so uncertainty analysis plays an important role in stability prediction. The identification results of the dominant modal parameters of the tool tip measured 10 times are listed in Table 3. The measured standard deviation shown in Table 3 would be a priori information for the related input parameters of the ZOA computational model. 5.3. Tool point FRF during cutting operation It is well known that, at high speeds gyroscopic moments cause separation of assembly modes into backward and forward ones [2,37]. FRFs of the tool tip during the cutting operation are different from those under the idle condition. Using FRFs under the idle condition to predict SLD at high speed may return significant errors: the differences in the natural frequency of dominant modes will be reflected in the inaccurate prediction of lobe locations, and the differences of modal damping will instead affect axial limit depths of cut [38]. In this paper, we only consider the uncertainties of input parameters measurement and uncertainties caused by simplified assumptions in the model. Therefore, to avoid the influence of spindle speeds on stability prediction, we need to calculate the tool tip FRFs during the cutting operation. Unfortunately, FRFs of the tool tip at high speeds are not obtainable via impact hammers. In recent studies, the inverse stability solution [2,39] has been presented as a novel and effective method and has been well applied. In this work, the tool-tip FRFs under the operational condition are calculated using:

1.2

-6 10 Measured FRF

Z

1

Fitted FRF

Amplitude(m/N)

X Y

0.8

PCB accelerometer

0.6

0.4 0.2 Hammer

0 200

(a)

400

600

800 1000 1200 Frequency(Hz)

1400

1600

1800

2000

(b)

Fig. 5. The impact tests: (a) Experimental test bench; (b) Measured FRFs, fitted FRFs and standard deviation via 10 tests.

10

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

-7 (a) 8 10 Measured Fitted

4 2

-6 -8

-2 200

400

600

-6 10

(c) 1

800 1000 Frequency(Hz)

1200

1400

0.5

-10 200

(d ) 0 I ma g ( m /N)

Re al ( m/ N)

Measured Fitted

0

-0.5

-1 200

400

(e) 1

600

800 1000 Frequency(Hz)

1200

1400

400

600

-7 10

800 1000 Frequency(Hz)

1200

1400

Measured Fitted -5

-10

-15 200

400

600

800 1000 Frequency(Hz)

1200

1400

(f) 1 Coherence

Coherence

Measured Fitted

-4

0

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

-7 10

-2

I ma g ( m /N)

Re al ( m/ N)

6

( b) 0

0

200 400 600 800 1000 1200 1400 1600 1800 2000 Frequency(Hz)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

200 400 600 800 1000 1200 1400 1600 1800 2000 Frequency(Hz)

Fig. 6. The actual measurement and fitting FRFs: (a) the real in x direction and (b) the imaginary in the x direction; (c) the real in the y direction and (d) the imaginary in the y direction; (e) coherence function in the x direction; (f) coherence function in the y direction.

Table 3 Measured boundaries of FRFs and standard deviation by 10 tests.

Lower boundary Upper boundary Std

xx (Hz)

nx (%)

xy (Hz)

ny (%)

897.13 906.24 7.32

2.02 2.46 0.61

932.71 940.94 7.66

3.75 4.54 0.56

 Experimental values of axial limit depths of cut and chatter frequency obtained by the Spindle Speed Ramp-up (SSR) tests [40].  Inverse calculation of input parameters by the ZOA computational model using Multi-objective Particle Swarm Optimization (MPSO) [41]. Compared with other heuristic algorithms such as genetic algorithm, MPSO is more efficient and achieves faster convergence. Modal parameters are calculated by minimizing the objective function f 0 , computed as:

f 1 ðiÞ ¼

blimexp ðiÞblimpre ðiÞ blimexp ðiÞ

f 2 ðiÞ ¼

xcexp ðiÞxcpre ðiÞ xcexp ðiÞ

ð26Þ

f0 ¼ k f1 k þ k f2 k where f 1 and f 2 are error vectors on axial limit depths of cut blim and chatter frequencies xc , respectively. i is the ith chatter limiting condition considered. pre represents the predicted value of analysis, and exp represents the experimental value. k k

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

11

is the 2-norm of the vector. The flow chart of the optimization algorithm is presented in Fig. 7. emin is a constraint tolerance, emin ¼ 1 e  6. In this paper, the spindle speeds of 7100 r/min and 7200 r/min were taken as examples to study. It is generally believed that the tool tip FRFs will hardly change if the spindle speed varies within 100 r/min [2]. The experimental results showed that when the spindle speed is 7100 r/min, the limit depth of cut is 0.7 mm, and the chatter frequency is 851 Hz; and that, when the spindle speed is 7200 r/min, the limit depth of cut is 0.6 mm and the chatter frequency is 864 Hz. According to the optimization objective function in Eq. (26), the process of optimization is as shown in Fig. 8; Fig. 8a shows the best fitness evolution for the input parameters, Generations = 400, MPSO algorithm converges and obtains the global optimal solution after 146 iterations. Fig. 8b shows the evolution of a 2-variable population over a known surface. This plotting function is useful for observing the performance of the swarm over a given 2-D search space. The errors between analytical and experimental results are 5.2% from Fig. 8. Table 4 shows the optimal input model parameters. The residues A and A in Eq. (9) depend on the mode shapes, and they can be assumed to be unaffected due to operating conditions [2]. Moreover, modal constants can be identified using experimentally obtained tool point FRFs at the idle state. According to the identified results, the dominant mode of the tool point FRFs for spindle speeds of 7100–7200 rpm was recalculated. The optimal input parameters shown in Table 4 were utilized to reorganize the corresponding FRFs based on the modal fitting technique. Fig. 9a and Fig. 9b show the tool point FRFs (500 Hz-1400 Hz) under the operation condition and the idle state in the x and y directions, respectively. As shown in Fig. 9, significant differences were observed in tool point dominant modes. The SLDs calculated from the tool tip FRFs under the idle condition were recalculated using the identified tool point FRFs under operational conditions (7100 r/min, 7200 r/min), as shown in Fig. 10. The recalculated stability diagram provided more accurate chatter prediction than the stability diagram obtained by using tool point FRFs at the idle state from Fig. 10. 5.4. Variations in cutting coefficients The cutting coefficient is a function of the workpiece material and tool performance. Any non-uniformity of the workpiece material can lead to uncertainties of cutting coefficient [6,12]. Because the parameters that were considered relatively noninfluential were still needed as input to the ZOA computational model, we ensured that the spindle speeds remain unchanged (7200 r/min) to estimate the cutting coefficients under the operational conditions. The cutting coefficients were calculated by the average cutting force method [42] with the feed rates of 0.01 mm/tooth, 0.015 mm/tooth, 0.02 mm/tooth,

SSR tests

limit depth of cut blim and chatter frequency c

ZOA computation model

Optimal x , x , y , y are calculated by MPSO

N

f0

min

?

Y Optimal input parameters Fig. 7. The flow chart of the MPSO optimization algorithm for input parameters.

12

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

Fig. 8. Optimization process of MPSO under the spindle speeds 7100 and 7200 rpm: (a) the best fitness evolution for the input parameters; (b) the performance of the swarm over a given 2-D search space.

Table 4 The optimal input parameters under the spindle speeds 7100–7200 rpm. K rc (N=m)

K tc (N=m)

xx (Hz)

nx (%)

xy (Hz)

ny (%)

6.213 E09

9.272 E09

865.69

2.68

875.26

2.16

0.025 mm/tooth and 0.03 mm/tooth. Linear expressions of average milling force as a function of feed per tooth can be derived from the data in Fig. 11; Figs. 11a and b show the milling forces in the x and y directions when the axial cutting depth is 0.5 mm, respectively. Figs. 11c and 11d show the milling forces in the x and y directions when the axial cutting depth is 1 mm, respectively. According to the calculation results of cutting coefficients in Fig. 11, there are obvious fluctuations. However, according to the conclusions in reference [42], the cutting coefficients are not theoretically affected by the change of axial depths of cut, so uncertain factors contribute to those fluctuations. At the same spindle speed, the mean values of radial and tangential cutting coefficients were identified by the 20 different axial cutting depths experiments shown in Table 5, which would be used as input parameters to the ZOA model. 6. Results and discussion In this section, we first discuss the prior distributions of the Bayesian method mentioned in Section 4.1. Taking the spindle speeds of 7100–7200 r/min as an example, Bayesian inference and IGWMCMC were used subsequently for uncertainty quantification of the most influential model parameters (along with the hyper-parameter r) under operational conditions, and 10-3

1.2

-3 10

1.5 Experimental tool tip FRF Indentified tool tip FRF

(a)

Experimental tool tip FRF Indentified tool tip FRF

(b)

Amplitude(mm/N)

Amplitude(mm/N)

1 0.8 0.6

0.4

1

0.5

0.2 0

0 0

200

400

600 800 Frequence(Hz)

1000

1200

1400

0

200

400

600 800 Frequence(Hz)

1000

1200

1400

Fig. 9. Tool point FRFs: (a) idle condition tool point FRF and calculated tool point FRF using the identified parameters for 7100–7200 rpm spindle speeds in the x direction; (b) idle condition tool point FRF and calculated tool point FRF using the identified parameters for 7100–7200 rpm spindle speeds in the y direction.

13

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532 5 Idle condition 7100rpm-7200rpm

4

Stable Chatter

blim (mm)

×

3

2

× × × ×

1

0 2000

3000

4000

5000

6000

×

7000

8000

9000

10000

(rpm) Fig. 10. SLDs obtained for the idle state and 7100–7200 rpm spindle speeds.

(b)

(a)

mean(F(y))

120 mean(F(x))

60

Fx Linear(Fx)

100 80

Fx=4160*f -8.8

Fy Linear(Fy)

50 40

Fy=2200*f -3.4

30

60 SSE: 44.21 R-square: 0.9723

SSE: 47.66 R-square: 0.974

20

40 0.01

0.014

(c) 320

0.018 0.022 f (mm/tooth)

0.026

0.01

0.03

0.014

(d)

Fx Linear(Fx)

120

0.018 0.022 f (mm/tooth)

0.026

0.03

Fy Linear(Fy)

280 260

Fx=4570*f +176.4

240

SSE: 36.79 R-square: 0.993

0.01

0.014

0.018 0.022 f (mm/tooth)

0.026

mean(F(y))

mean(F(x))

300 110 100 Fy=2734*f +40.52

90 80

SSE: 42.85 R-square: 0.9776

70 0.03

0.01

0.014

0.018 0.022 f (mm/tooth)

0.026

0.03

Fig. 11. Linear regression diagram showing the change of the x and y directions in average milling force with feed per tooth: (a) milling force in the x direction when the axial cutting depth is 0.5 mm and (b) milling force in the y direction when the axial cutting depth is 0.5 mm; (c) milling force in the x direction when the axial cutting depth is 1 mm and (b) milling force in the y direction when the axial cutting depth is 1 mm.

the propagation of uncertainties of posterior distributions through the model was discussed. Finally, utilizing the milling system in Fig. 15, different combinations of axial depths of cut and spindle speeds were tested to see if chatter occurred and the prediction results of the proposed method were verified.

6.1. Priori information of Bayesian inference The determination of Bayesian prior information plays an important role in the calculation and convergence of posterior distributions. For the natural frequencies and damping ratios of the input parameters in the x and y directions, the prior distributions are assumed to be normal distribution, the optimization results in Section 5.3 are taken as the mean value, and the standard deviation measured repeatedly in Section 5.2 are the standard deviations of the normal distribution. The parameters that are considered relatively non-influential (K rc , and K tc ) are still needed as the input to the model. The prior distributions of cutting coefficients can be assumed to be uniform distribution. These prior distributions of input parameters

14

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532 Table 5 The mean cutting coefficients identified by 20 cutting tests. Cutting coefficients

K rc (N=m)

K tc (N=m)

Mean

5.529 E09

8.767 E09

are determined by experience based on the results of many tests by manufacturing big data centers built by the National Engineering Research Center for Digital Manufacturing Equipment. Table 6 summarizes the distributions assigned to each of the input parameters. 6.2. Uncertainty quantification by Bayesian-IGWMCMC The input parameters (cutting coefficients, and modal parameters of the tool tip) of the milling stability prediction model are uncertain. Uncertainties in the computational model can be attributed to many factors, including but not limited to: (i) measurement errors, (ii) numerical errors, and (iii) simplification of boundary conditions, which lead to the uncertainties of the SLDs. For robust chatter predictions, these uncertainties must be considered. Bayesian inference was used to quantify uncertainties, and IGWMCMC was applied to estimate the posterior probability density distributions of the input parameters and hyper-parameter r; IGWMCMC algorithm was set to run for 50,000 samples, thinning by taking every 5th sample and a burn-in of 10,000 samples was used, as shown in Fig. 12, the acceptance rate/fraction is 37.6%. To diagnose convergence of Markov chain, we have applied the variance ratio of each parameter to diagnose convergence [43]; the results showed that the variance ratio of each parameter is close to 1 after 30,000 iterations. Fig. 12a shows the autocorrelation plot of the input parameters. When the autocorrelation of input parameters approaches zero, it is shown that the Markov chains converge [35]. The full posterior distributions of all input parameters are shown in Fig. 12b. Fig. 13 shows the prior and posterior densities of a set of influence parameters (along with parameter r). As seen from Fig. 13, the standard deviation of the posterior distributions of the damping ratios are obviously larger than that of the prior distributions, which indicates that the uncertainties of the tool tip damping ratios under the cutting condition are obviously greater than those under the idle condition. The mean and standard deviation of posterior distributions for the most influential model parameters are shown in Table 7. 6.3. Uncertainty propagation and experimental verification The posterior distributions of uncertainties calculated by Bayesian inference were propagated through the model to make reliable limit depths of cut blim predictions. Monte Carlo (MC) based sampling techniques were used to propagate uncertainties through the model. In every MC iteration, by sampling from their respective posterior distributions, a value (represented by a solid red circle) was assigned to each of the most influential model parameter, as shown in Fig. 4. The SLDs obtained after 50 iterations are shown in Fig. 14. In Fig. 14, red and blue solid lines represented the upper and lower bounds of SLDs, while green solid line represents the SLD with maximum probability, which was the best estimate of all of the most influential parameters for posterior distributions. To verify the SLDs in Fig. 14, we ensured that the feed rate (0.02 mm/tooth) and radial width of cut were constant (8 mm) during end milling; the spindle speeds range was from 2000 to 20000 rpm, and the spindle speeds increased by 100 rpm increments. For each spindle speed, the corresponding axial depths of cut varied by 0.1 mm each time until chatter occurred. Fig. 15a shows the milling test bench for different axial depths of cut, Fig. 15b represents the data acquisition system during the tests, Charge Amplifiers and Personal Computer (PC) is shown in Fig. 15c. Fig. 16 shows the cutting force signals collected by two different axial depths of cut when the spindle speed is 7800 r/min. Fig. 16a and Fig. 16d show the time domain signals, and Short Time Fourier Transform (STFT) of the time domain signals are shown in Fig. 16b and Fig. 16e. Fig. 16c and Fig. 16f show the Fast Fourier Transform (FFT) frequency spectrum when the axial depth of cut is 0.4 mm and 0.5 mm, respectively. As seen from Fig. 16c and Fig. 17, the limit axial depth of cut is 0.5 mm when the spindle speed is 7800 rpm. Table 6 Priori distributions hypothesis assigned to each of the parameter for Bayesian uncertainty quantification. Type

Parameter

Distribution

Input parameters

K rc



U ð5  10 ; 7  10 Þ

N=m

K tc

U ð8  109 ; 10  109 Þ

N=m

xx nx

xy ny Model discrepancy hyper-parameter

r

9



Units 9



Hz



–

N ð865:69; 7:32Þ N ð2:68; 0:61Þ

Hz



–



–

N ð875:26; 7:66Þ N ð2:16; 0:56Þ U ð0; 1Þ

15

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

(a)

Markov Chain Autocorrelation

1

K rc K tc x

0.8

x

y y

autocorrelation

0.6

0.4

0.2

0 Effective Sample Size (ESS): 689

-0.2 0

20

40

60

80

100

120

140

160

y

y

x

x

K tc

(b)

K rc

K tc

x

x

y

y

Fig. 12. Posterior distributions of model input parameters: (a) the autocorrelation plot of input parameters; (b) the full posterior distributions of all input parameters.

A macro-scale image (amplified 100) of the machined surface was obtained using a portable microscope, as shown in Fig. 17. It can be seen from Fig. 17a that the cutting was stable when the axial depth of cut was 0.4 mm, and chatter occurred when the axial depth of cut was 0.5 mm, as shown in Fig. 17b. The experimental chatter detection results using same method are shown in Fig. 14. From Fig. 14, it can be concluded that the predicted results have good reliability in the range of spindle speeds of 6000–9000 rpm. When the cutting depth is lower than the lower boundary, it must be stable, and when the cutting depth is above the upper boundary, chatter will occur. There is an uncertain region between the upper boundary and the lower boundary. In addition, from the experimental results, it seems that there is a maximum probability that the limit depths of cut fall near the green solid line. Moreover, for the spindle speeds less than 6000 rpm, it can be found that the predicted results are no longer consistent with the experimental detection results, which may be due to the spindle speeds having a significant effect on the FRFs of the tool tip at this moment.

6.4. Comparisons with the results in [18] To further prove the superiority of the proposed method, we compared the results of Bayesian inference with those method obtained by the method in reference [18] under the same experimental conditions. In reference [18], a robust spindle speed optimization formulation was developed, in which the upper bound of surface location error and lower bound of

16

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

(b) Prior distribution Posterior distribution

0.08 0.06 0.04

0

845 850 855 860 865 870 875 880 885 890 895 x (Hz)

(c)

Prior distribution Posterior distribution

110

0

-2

-1

0

1

(d)

2

3

cx (%)

4

0.06

0.04

0.02

5

6

7

Prior distribution Posterior distribution

80 Probability density

0.08 Probability density

115

105

0.02

0

Prior distribution Posterior distribution

120

Probability density

Probability density

(a) 0.1

75

70

65

850 855 860 865 870 875 880 885 890 895

0

-4

-2

y(Hz)

(e)

0

2

c y(%)

4

6

8

Probability density

4

3

2

1

0

-6

-4 -2 0 variance hyperparameter(

2

lo g (

)

4 )

6

Fig. 13. Prior and posterior densities of (a) the natural frequency of the tool tip in the x direction; (b) damping ratio of the tool tip in the x direction; (c) the natural frequency of the tool tip in the y direction; (d) the damping ratio of the tool tip in the y direction; (d) variance hyper-parameter (logðrÞ).

Table 7 Mean and standard deviation of the posterior distribution for each input parameter. Parameter

xx (Hz)

nx (%)

xy (Hz)

ny (%)

logðrÞ

Mean Std

868.32 5.02

2.72 2.36

877.11 5.13

2.12 1.91

3.862 1.282

SLDs were adopted as the optimization object and the constraint condition, respectively. The optimization problems were solved by an augmented Lagrange function method. The comparison results are shown in Fig. 18. The results show that our work provided a narrower uncertain region with accurate upper and lower boundaries, which provided more useful suggestions for parameter selection in the actual milling process at the workshop level. SLDs.

17

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532 5 chatter

4.5

no chatter

4

Upper bound Lower bound MC sample The Maximum probability

blim ( mm)

3.5 3 2.5 2

Probability density

1.5 1 0.5

0 2000

3000

4000

5000

6000

7000

8000

9000

10000

(rpm)

Fig. 14. SLDs calculated after uncertainties propagate by Monte Carlo.

a

b

Dewesoft

c Workpiece

dynamometer Charge amplifier

PC

Fig. 15. Milling test setup for different axial depths of cut: (a) Milling test bench; (b) Data acquisition system and (c) Charge Amplifier and PC.

7. Conclusions In this paper, a novel perspective of uncertainty quantification based on Bayesian inference for SLDs was proposed. The parameter estimation problem in this study was to infer the posterior distributions of all of the most influential input parameters, which was achieved by an improved Ensemble Markov Chain Monte Carlo (IGWMCMC) algorithm. The main contributions of this paper include: (1) In addition to the uncertainties caused by the input parameters K rc ; K tc ; xx ; nx ; xy ; ny , some epistemic uncertainties caused by simplified assumptions in the ZOA computational model proposed by Altintas and Budak were also considered. (2) Different input parameters have varying degrees of influence on the uncertainties of the model output. To reduce computational costs, an advanced Morris sensitivity analysis method was used to select the most influential parameters. (3) The fact that the FRFs of the tool tip during the cutting operation are different from those under the idle condition was considered. To obtain more robust SLDs, the MPSO optimization algorithm was used to obtain the prior distributions of modal parameters of the tool tip under operational conditions instead of the idle state. (4) The uncertain region obtained by the proposed method are narrower than those in reference [18], which can provide more useful suggestions for the optimization of cutting parameters in the actual milling process at the workshop level.

18

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

(a) 150

(d) 200 150 Amplitude (N)

50 0

-50

100 50

0 -50

-100

-100

Frequency (kHz)

(b) 5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (s) Spindle speed=7800 r/min

-200

1

0 20

4

0

3

-20 -40

2

-60 1

(e) 5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (s) Spindle speed=7800 r/min

1 20

4

0

3

-20 -40

2

-60 1 -80

-80 0

100 200 300 400 500 600 Time (ms)

(c) 70

700 800

0

900

900

f

60 Amplitude (N)

Amplitude (N)

60 50 7f

40 30

X 936 Y 3.97

3f

20 8f

10 0

100 200 300 400 500 600 700 800 Time (ms)

(f) 70

f

Power/frequency (dB/Hz)

0

Power/frequency (dB/Hz)

-150

-150

Frequency (kHz)

Amplitude (N)

100

0

1000

50 7f 40

X 935 Y 9.027

3f

30

8f

20 10

2000 3000 Frequency (Hz)

4000

0

5000

0

1000

2000 3000 Frequency (Hz)

4000

5000

Fig. 16. Cutting force signals of two different axial depths of cut; (a) time domain for an axial depth of cut of 0.4 mm; (b) the time–frequency spectrograph using STFT for an axial depth of cut of 0.4 mm; (c) the frequency domain FFT for an axial depth of cut of 0.4 mm; (d) the time domain for an axial depth of cut of 0.5 mm; (e) the time–frequency spectrograph using STFT for axial depth of cut 0.5 mm; (f) the frequency domain FFT for an axial depth of cut of 0.5 mm.

(a)

(b)

100

100

Fig. 17. The surface quality photographs of workpiece taken by portable microscope: (a) axial depth of cut of 0.4 mm; (b) axial depth of cut of 0.5 mm.

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

19

5 chatter

4.5

no chatter

4

The proposed method

blim ( mm)

3.5

The proposed in [18]

3 2.5 2 1.5

1 0.5 0 2000

3000

4000

5000

6000

7000

8000

9000

10000

(rpm) Fig. 18. Comparisons between the proposed method and that of reference [18] for calculating.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research is supported by the National Natural Science Foundation of China under Grant nos. 51875224, 51705174 and Major special projects in Jiangsu Province of China under Grant nos. SBE2017020146. Appendix

Algorithm—An improved Ensemble Markov Chain Monte Carlo (GWMCMC) [44] 1: for j 2 f0; 1g do 2: for k ¼ 1;    ; K=2 do 3: Draw a walker X i at random from the complementary ensemble W ð jÞ ðtÞ 4: 5: 6: 7: 8: 7: 9:

ðjÞ

Sk v V gðv Þ,Eq.(23) Y X i þ v ½X k ðtÞ  X i  v N1 pðYÞ=pðX k ðtÞÞ q r R ½0; 1 if r q,Eq.(24) then Y X k ðt þ 12Þ Xk

10: else X k ðt þ 12Þ 11: end if 12: end for13: t

X k ðtÞ t þ 1214: end for

References [1] M. Kayhan, E. Budak, An experimental investigation of chatter effects on tool life, Proc. Inst. Mech. Eng. Part B-J. Eng. Manuf. 223 (11) (2009) 1455– 1463. [2] O. Özsßahin, E. Budak, H.N. Özgüven, In-process tool point FRF identification under operational conditions using inverse stability solution, Int. J. Mach. Tools Manuf. 89 (2015) 64–73. [3] E. Graham, M. Mehrpouya, R. Nagamune, S.S. Park, Robust prediction of chatter stability in micro milling comparing edge theorem and LMI, CIRP Ann. 7 (1) (2014) 29–39. [4] D. Hajdu, T. Insperger, G. Stépán, Quantification of uncertainty in machining operations based on probabilistic and robust approaches, Procedia CIRP 77 (2018) 82–85. [5] M. Löser, A. Otto, S. Ihlenfeldt, G. Radons, Chatter prediction for uncertain parameters, Adv. Manuf. 6 (3) (2018) 319–333.

20

K. Li et al. / Mechanical Systems and Signal Processing 138 (2020) 106532

[6] E. Gözü, Y. Karpat, Uncertainty analysis of force coefficients during micromilling of titanium alloy, Int. J. Adv. Manuf. Technol. 93 (1-4) (2017) 839–855. [7] G. Totis, RCPM—A new method for robust chatter prediction in milling, Int. J. Mach. Tools Manuf. 49 (3-4) (2009) 273–284. [8] D. Hajdu, T. Insperger, D. Bachrathy, G. Stepan, Prediction of robust stability boundaries for milling operations with extended multi-frequency solution and structured singular values, J. Manuf. Process. 30 (2017) 281–289. [9] Y. Liu, T. Li, K. Liu, Y. Zhang, Chatter reliability prediction of turning process system with uncertainties, Mech. Syst. Signal Process. 66 (2016) 232–247. [10] Y. Liu, Z. Wang, K. Liu, Y. Zhang, Chatter stability prediction in milling using time-varying uncertainties, Int. J. Adv. Manuf. Technol. 89 (9-12) (2017) 2627–2636. [11] X. Huang, Y. Zhang, C. Lv, Probabilistic analysis of dynamic stability for milling process, Nonlinear Dyn. 86 (3) (2016) 2105–2114. [12] M.T. Loukil, V. Gagnol, T.P. Le, Reliability evaluation of machining stability prediction, Int. J. Adv. Manuf. Technol. 93 (1-4) (2017) 337–345. [13] E. Graham, M. Mehrpouya, S.S. Park, Robust prediction of chatter stability in milling based on the analytical chatter stability, J. Manuf. Process. 15 (4) (2013) 508–517. [14] H. De Gersem, D. Moens, W. Desme, DirkVandepitte, A fuzzy finite element procedure for the calculation of uncertain frequency response functions of damped structures: Part 2—Numerical case studies, J. Sound Vibr. 288 (3) (2005) 463–486. [15] M. Hanss, K. Willner, Fuzzy arithmetical modeling and simulation of vibrating structures with uncertain parameters, Shock. Vib. 38 (4) (2006) 374– 375. [16] N.D. Sims, G. Manson, B. Mann, Fuzzy stability analysis of regenerative chatter in milling, J. Sound Vibr. 329 (8) (2010) 1025–1041. [17] D. Hamann, N.P. Walz, A. Fischer, M. Hanss, P. Eberhard, Fuzzy arithmetical stability analysis of uncertain machining systems, Mech. Syst. Signal Process. 98 (2018) 534–547. [18] X. Zhang, L. Zhu, D. Zhang, H. Ding, Y. Xiong, Numerical robust optimization of spindle speed for milling process with uncertainties, Int. J. Mach. Tools Manuf. 61 (2012) 9–19. [19] G. Duncan, T.L. Schmitz, M.H. Kurdi, Uncertainty propagation for selected analytical milling stability limit analyses, Soci. Manuf, Eng, 2000. [20] Y. Altintasß, E. Budak, Analytical prediction of stability lobes in milling, CIRP Ann-Manuf. Technol. 44 (1) (1995) 357–362. [21] P. Angelikopoulos, C. Papadimitriou, P. Koumoutsakos, Bayesian uncertainty quantification and propagation in molecular dynamics simulations: a high performance computing framework, J. Chem. Phys. 137 (14) (2012) 144103. [22] S.R. Yeratapally, M.G. Glavicic, C. Argyrakis, M.D. Sangid, Bayesian uncertainty quantification and propagation for validation of a microstructure sensitive model for prediction of fatigue crack initiation, Reliab. Eng. Syst. Saf. 164 (2017) 110–123. [23] J. Goodman, J. Weare, Ensemble samplers with affine invariance, Commun. Appl. Math. Comput. Sci. 5 (1) (2010) 65–80. [24] J. Tlusty, F. Ismail, Basic non-linearity in machining chatter, CIRP Ann-Manuf. Technol. 30 (1) (1981) 299–304. [25] E. Budak, Y. Altintas, Analytical prediction of chatter stability in milling—part I: general formulation, J. Dyn. Syst. Meas. Control-Trans. ASME 120 (1) (1998) 22–30. [26] E. Budak, Y. Altintas, Analytical prediction of chatter stability in milling—part II: application of the general formulation to common milling systems, J. Dyn. Syst. Meas. Control-Trans. ASME 120 (1) (1998) 31–36. [27] D.J. Ewins, Modal testing: theory and practice, Vol. 15, Research studies press, Letchworth, 1984. [28] A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto, Sensitivity Analysis in Practice, John Wiley & Sons, 2004. [29] K. Worden, J.J. Hensman, Parameter estimation and model selection for a class of hysteretic systems using Bayesian inference, Mech. Syst. Signal Process. 32 (2012) 153–169. [30] G.A. Ortiz, D.A. Alvarez, D. Bedoya-Ruíz, Identification of Bouc-Wen type models using the transitional Markov chain Monte Carlo method, Comput. Struct. 146 (2015) 252–269. [31] P.E. Hadjidoukas, P. Angelikopoulos, D. Rossinelli, D. Alexeeva, C. Papadimitrioub, P. Koumoutsakos, Bayesian uncertainty quantification and propagation for discrete element simulations of granular materials, Comput. Meth. Appl. Mech. Eng. 282 (2014) 218–238. [32] J.M. Nichols, W.A. Link, K.D. Murphy, C.C. Olson, A Bayesian approach to identifying structural nonlinearity using free-decay response: application to damage detection in composites, J. Sound Vibr. 329 (15) (2010) 2995–3007. [33] H.F. Lam, J. Yang, S.K. Au, Bayesian model updating of a coupled-slab system using field test data utilizing an enhanced Markov chain Monte Carlo simulation algorithm, Eng. Struct. 102 (2015) 144–155. [34] R.C. Smith, Uncertainty quantification: theory, implementation, and applications. Vol. 12. Siam, 2013. [35] J. Christen, C. Fox, A general purpose scale-independent MCMC algorithm. Preprint (2007). [36] S. Sharma, Markov Chain Monte Carlo methods for Bayesian data analysis in astronomy, Annu. Rev. Astron. Astrophys. 55 (2017) 213–259. [37] H. Cao, B. Li, Z. He, Chatter stability of milling with speed-varying dynamics of spindles, Int. J. Mach. Tools Manuf. 52 (1) (2012) 50–58. [38] N. Grossi, L. Sallese, A. Scippa, G. Campatelli, Improved experimental-analytical approach to compute speed-varying tool-tip FRF, Precis. Eng. 48 (2017) 114–122. [39] M. Postel, O. Özsahin, Y. Altintas, High speed tooltip FRF predictions of arbitrary tool-holder combinations based on operational spindle identification, Int. J. Mach. Tools Manuf. 129 (2018) 48–60. [40] N. Grossi, A. Scippa, L. Sallese, R. Sato, G. Campatelli, Spindle speed ramp-up test: A novel experimental approach for chatter stability detection, Int. J. Mach. Tools Manuf. 89 (2015) 221–230. [41] S Mostaghim, J Teich, Strategies for finding good local guides in multi-objective particle swarm optimization (MOPSO). IEEE Service Center, Inidanapolis, Indiana, USA, 2003. April 26–33. [42] M. Wang, L. Gao, Y. Zheng, An examination of the fundamental mechanics of cutting force coefficients, Int. J. Mach. Tools Manuf. 78 (2014) 1–7. [43] A.E. Raftery, S.M. Lewis, [Practical Markov Chain Monte Carlo]: comment: one long run with diagnostics: implementation strategies for Markov Chain Monte Carlo, Stat. Sci. 7 (4) (1992) 493–497. [44] D. Foreman-Mackey, D.W. Hogg, D. Lang, J. Goodman, emcee: the MCMC hammer, Publ. Astron. Soc. Pac. 125 (925) (2013) 306.