Estimation of vibration stability in turning using operational modal analysis

Mechanical Systems and Signal Processing 130 (2019) 315–332

Contents lists available at ScienceDirect

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

Estimation of vibration stability in turning using operational modal analysis S. Kim, K. Ahmadi ⇑ Department of Mechanical Engineering, University of Victoria, Victoria, British Columbia V8W 2Y2, Canada

a r t i c l e

i n f o

Article history: Received 12 November 2018 Received in revised form 16 April 2019 Accepted 25 April 2019 Available online 15 May 2019 Keywords: Chatter Operational modal analysis

a b s t r a c t Unstable vibrations in machining operations, known as chatter, may cause damages to the tool, workpiece, or the machine tool. Chatter is usually avoided by selecting the appropriate spindle speed and depth of cut determined by the dynamic models of the machining process. Due to modeling uncertainties or variations of the process dynamics, chatter may happen even when stable cutting conditions are used in the process planning stage. Therefore, it is critical to detect chatter during the process and take corrective actions in the shortest possible time. The existing chatter detection methods process various signals such as sound, vibration, and force to extract chatter indicators during the process, however these methods detect chatter only after vibrations become unstable. In this paper, a new chatter detection method is presented based on the identification of the dynamics of turning processes using Operational Modal Analysis (OMA). Because this new method estimates the stability margin of the process, rather than the stability/instability of vibrations, it determines the level of stability of the cut before vibrations become unstable. The performance of the presented chatter detection method is studied using numerical simulations and machining experiments. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Unstable vibrations, known as regenerative chatter, may arise during machining operations due to the inherent self-excitation mechanism in the chip generation process [1]. Vibrations of the tool or the workpiece during one tool (or workpiece) revolution leave waves on the machined surface, which in turn modulate the cutting forces in the succeeding revolution. This feedback between vibrations in subsequent cuts result in self-excited vibrations during machining processes. If the spindle speed and depth of cut are not selected properly, such self-excited vibrations may become unstable and damage the tool and the workpiece. Stable spindle speeds and depth of cuts are usually selected according to the vibration models of the machining process. The dynamics of machining processes is described by infinite dimensional distributed parameter systems with a delay term in their equation of motion [2,3]. Various time and frequency domain methods have been presented to determine the stability of such systems when arbitrary spindle speed and depth of cuts are used. A review of the existing chatter analysis methods is available in [4,5]. Frequency domain methods such as [6–8] develop a closed loop model of the machining dynamics by incorporating the linear models of the cutting forces in the models of the machine tool dynamics, which are usually represented by the Frequency Response Function (FRF) at tool and workpiece contact points. The stability of ⇑ Corresponding author. E-mail address: [email protected] (K. Ahmadi). https://doi.org/10.1016/j.ymssp.2019.04.057 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

316

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

vibrations are then determined according to classical theories of the stability of closed loop systems [9]. Time domain methods such as [10–14] employ various discretization schemes to develop a finite dimensional approximation of the distributed parameter system representing machining dynamics. Subsequently, the stability of vibrations is determined with respect to the poles of the resulting lumped system. In both time and frequency domain approaches, the results of the stability analysis are presented in the form of Stability Lobe Diagrams (SLD) showing the border between the stable and unstable depth of cuts at each spindle speed. Although selecting stable spindle speed and depth of cut according to SLD has proved to be an effective method to avoid chatter without sacrificing productivity, SLD can be imprecise for various reasons such as un-modeled phenomena or the variation of the system parameters during the process (e.g. increased tool wear). Therefore, even when the spindle speed and depth of cut are selected according to SLD, it is critical to monitor the stability of the system during the process to avoid the damages that may occur because of chatter. The existing chatter detection methods monitor various signals such as vibration, force, sound, and acoustic emissions during the process to extract chatter indicators when vibrations become unstable [15]. Most of these methods include some form of signal processing that is performed either online during the process [16–18], or offline after completing the process [19,20]. Nonetheless, the vast majority of the existing methods detect chatter only after vibrations become unstable. Although the timely detection of chatter after its occurrence is important for preventing severe damages to the tool or the machine, damages to the workpiece will not be prevented by the existing methods, because chatter marks will leave an imprint on the finished surface. In this paper, a new approach for monitoring chatter is presented. In the presented approach, the stability margin of the system when it is still stable is determined from the vibrations measured during the process. A similar approach has recently been presented by Kiss et al. [21], where milling dynamics has been approximated by a lumped system, and then an experimental method is used to identify the dominant characteristic multiplier of the system from its response to an impulse applied during the milling process. Monitoring the modulus of the dominant characteristic multiplier of the stable system can be used to forecast the loss of stability of the system before it becomes unstable. In the case study presented in [21], the workpiece is the flexible component and the tool is assumed to be rigid. Therefore, the impulse is applied to the stationary workpiece and the response is also measured using accelerometers installed on the workpiece. Although the method in [21] can effectively determine the stability margin of the system while it is stable, applying impulse on the tool or the workpiece during the process may introduce technical challenges in practical applications. In the method presented in this paper, Operational Modal Analysis (OMA) is used to identify the dominant poles of the machining system when it is subjected to the ambient random excitation that inherently exists during the process[22–28]. Conventionally, poles of mechanical systems are identified by conducting Experimental Modal Analysis (EMA) [29–31] where controlled excitation (e.g. impulse hammer) to the structure and the resulting response (e.g. acceleration) are measured and used in various time [22,32–35] and frequency [36–40] domain methods to identify the poles (or modal parameters) of the system. Alternatively, OMA can be used when the measurement of the excitation force is not possible for practical reasons, and therefore only the response of the system to unmeasured operational loads is used to identify its poles. Because OMA only relies on the response of the system to random ambient excitation, which inherently exists during machining operations, the need to apply an impulse during the process is eliminated by using the presented OMA-based approach. The OMA of the measured vibration signals during the turning process leads to the identification of the frequency and damping of the dominant poles of the system. The damping of the dominant pole is then used to monitor the stability of the system while it is still stable. When the system moves towards the border of stability, the damping ratio of its dominant pole goes towards zero, signaling the loss of stability before the process becomes unstable. In the next section, the dynamics of turning processes is modeled as a distributed parameter system. In Section 3, Semi Discretization Method [41] is used to approximate the distributed parameter system with a finite dimensional lumped system. In Section 4 the two OMA methods that are used to identify the dynamics of the equivalent lumped system are explained. These methods are used in the numerical and experimental case studies described in Section 5 to study the performance of the presented method in quantifying the stability of the system. 2. Turning dynamics A schematic of generic turning operations is shown in Fig. 1 (a), and an image of the turning setup that is used in this paper is shown in part (b) of the figure. Also shown in Fig. 1 (a) is the Multi Degree Of Freedom (MDOF) system describing the dynamics of the tool in the feed (z) and cutting (y) directions. The tool in axial direction (x) and the workpiece in all of the three directions are assumed to be rigid. The tool deflection in the feed direction is represented by the deflections at M discrete points (z1 ; . . . ; zM ) and its deflection in the cutting direction (y) is represented by the deflections at N discrete points (y1 ; . . . ; yN ). In total, the model has D ¼ M þ N degrees of freedom. The equation governing the motion of the MDOF system is as follows:

€ þ Cy_ þ Ky ¼ f ðtÞ My

ð1Þ

where y contains the deflections at the D degrees of freedom

yðt Þ ¼ ½ z1 ðt Þ . . . zM ðt Þ y1 ðt Þ . . . yN ðtÞ T

ð2Þ

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

317

Fig. 1. (a) MDOF model of turning dynamics, and schematic of the test setup used in Section 5, (b) image of the turning setup, and (c) closed loop diagram describing self-excited vibrations in turning.

and M; K, and C are the corresponding D  D mass, stiffness, and damping matrices of the MDOF system, respectively. The force vector, f ðtÞ, consists of the cutting forces in feed (f z ðtÞ) and cutting (f y ðtÞ) directions that are applied at z1 and y1 :


f ðt Þ ¼ f z ðtÞ 01ðM1Þ

f y ðt Þ 01ðN1Þ

T

ð3Þ

The cutting forces are proportional to the thickness of the uncut chip, which is generated by the constant feed motion of the tool in every revolution of the workpiece (feedrate [mm/rev]). This constant chip thickness is also modulated by the oscillations of the tool around its equilibrium position. Therefore, the total cutting forces comprise a static component, generated by the feed motion, and the dynamic part, generated by tool vibrations. The static component of the forces do not affect the stability of the system, and its effect can easily be superimposed to the dynamic solution, therefore it is neglected in what follows. The dynamic component of the cutting forces is expressed in terms of the instantaneous, yðtÞ, and delayed, yðt  T Þ, vibrations as follows:

f ðtÞ ¼ Kc afyðt  T Þ  yðt Þg 8 > < K f i ¼ 1; j ¼ 1
 Kc ¼ kij ; kij ¼ K t i ¼ M þ 1; j ¼ 1 > : 0 otherwise

ð4Þ

where Kc is a square D  D matrix, and K f and K t are the constant cutting force coefficients in feed (z) and cutting (y) directions, respectively. These coefficients are usually determined empirically for each material and cutting tool [1]. The depth of cut, shown in Fig. 1(a), is denoted a. The delay term (yðt  T Þ  yðtÞ) in Eq. (4) represents the modulation of the chip thickness due to the vibrations in the feed direction (zðt Þ) and the waviness left on the surface by the vibrations at the previous pass (zðt  T Þ), where T is the spindle revolution period. Modulation of the constant chip thickness (feed per revolution) by the oscillations in the feed direction is depicted in Fig. 1(a). Substituting the cutting force from Eq. (4) in the governing equation of the MDOF system in Eq. (1) results in the following Delay Differential Equation describing the dynamics of self-excited (regenerative) vibrations in turning.

 ¼ Kc ayðt  T Þ; € þ Cy_ þ Ky My

K ¼ K þ Kc a

ð5Þ

Presence of the delay term on the right hand side of Eq. (5) converts the lumped MDOF system into a distributed parameter system. Part (c) of Fig. 1 shows the block diagram representation of the system in Eq. (5). Depending on the values of the depth of cut, a, and spindle revolution period, T, the closed loop system in Eq. (5) may become unstable and lead to regenerative chatter.

318

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

3. Discrete Time Domain model Various discretization schemes have been used in the literature to approximate the distributed parameter system in Eq. (5) by a finite dimensional lumped system [14,41]. In this paper, we use the Semi-Discretization Method (SDM)[41], which converts the DDE in Eq. (5) into a system of Ordinary Differential Equations. In this method, the duration of the delay period   (T) is discretized into equal time segments of length h such that T ¼ r þ 12 h, where r is an integer number. During each time interval ti 6 t < tiþ1 , the delay term in Eq. (5) is assumed to be constant yðtir Þ, converting the delay differential equation into an ordinary differential equation in each time interval:

 ¼ Kc ayðtir Þ; € þ Cy_ þ Ky My

t i 6 t < t iþ1

ð6Þ

It was shown in [41] that when r goes to infinity and h goes to zero, Eq. (6) converges to the original DDE in Eq. (5). The solution of Eq. (6) after applying the boundary values of yðt Þ at t i and t iþ1 results in the following equation that describes the finite dimensional equivalent of the distributed parameter system in Eq. (5):

uiþ1 ¼ Gui yi ¼ Pui ;

P ¼ ½ I 0    0 DðD:ðrþ2ÞÞ

ð7Þ

where I and 0 are D  D identity and zero matrices, and ui is the discrete time state vector consisting of the displacement and velocities at time steps within one delay period (i.e. spindle revolution period, T):


ui ¼ yTi

y_ Ti

yTi1

yTi2

   yTir

T

; yi ¼ yðt i Þ

ð8Þ

In this paper, the superscript T stands for the transpose of a matrix (or vector). In Eq. (7), G is the state transition matrix that advances the state vector, ui , by one time step to uiþ1 . The transition matrix depends on the mass, stiffness, and damping matrices of the system, depth of cut, cutting force coefficients, and the discretization time interval. The detailed composition of the transition matrix, G, is described in the Appendix A. Eq. (7) describes the free vibrations (or impulse response) of the regenerative system in discrete time domain. The impulse (or free) response of lumped systems can be expressed in modal form in terms of the poles of the system [9]:

yi ¼

2P X cp ap lip ;

lp ¼ ekp h ; ap ¼ Pvp

ð9Þ

p¼1

where lp are the eigenvalues of the transition matrix, h is sampling period, and the vectors v p are the corresponding eigenvectors. The eigenvalues, lp , represent the discrete time poles of the system, and ap represent the mode shapes. The correqffiffiffiffiffiffiffiffiffiffiffiffiffi sponding continues time poles are kp ¼ fp xp þ jxp 1  f2p , where xp and fp are the modal frequencies and damping ratios associated with each pole. The constants, cp , depend on the initial conditions of the free decay. The total number of poles considered in approximating the Impulse Response (IR) is 2P, consisting of P pairs of complex conjugate modes. The stability of the discrete time system in Eq. (7) is determined with respect to the modulus of the discrete time poles, lp . If all of the eigenvalues of the transition matrix (i.e. discrete time poles) are located inside of the unit circle on the complex plane, the system is stable, otherwise it is unstable. Equivalently, if the damping ratio of any of the corresponding continues time poles, fp , are negative, the system is unstable. Therefore, in this work we consider the smallest modal damping ratio of the system as the measure of its stability. The reduction of the smallest modal damping ratio indicates the declining trend of the stability of the machining process. Chatter prediction methods such as [10] employ the mass, stiffness, and damping matrices of the tool model along with the cutting force coefficients in Eq. (7) to determine the stability of vibrations based on the eigenvalues of the resulting transition matrix. In an inverse approach, in this work we will use System Identification methods to determine the frequency and damping ratio of the modes of the system from vibration signals measured directly during stable turning processes. However, the measurement of the IR of the system (yi ) during machining operations (e.g. using impulse hammer tests) is not practical and thus the IR in Eqs. (7) and (9) are approximated using the Correlation Function (CF) of the vibration signals measured during the process. When the lumped system in Eq. (7) is subjected to white noise excitation, the Correlation Function (CF) of the response is a scaled equivalent of the system IR (or free decay) [23]. Let us assume that the system in Eq. (7) is subjected to normally i ¼ y  ðt i Þ, is measured at i ¼ 1 . . . L discrete time instants. distributed white noise excitation, and the response of the system, y The CF of the measured response, Ry ðiÞ, in discrete time domain is defined as follows:

Ry ðiÞ ¼

Li 1 X k y Tkþi ; y ðL  iÞ k¼1

i ¼ 0L  1

ð10Þ

Each column of Ry ðiÞin Eq. (10) is a scaled equivalent of the response, yi , to an impulse applied at the corresponding DOF [23]. Therefore, similar to the IR in Eq. (9), the CF can also be decomposed in terms of the poles of the system in Eq. (7):

Ry ðiÞ ¼

2P X Bp lip ; p¼1

T

Bp ¼ bp bp

ð11Þ

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

319

where bp is the deflection shape vector corresponding to the mode shapes of the system, ap . In the next section, two OMA methods are used to identify the poles (lp ) of the lumped system in Eq. (11) from the CF of the vibration signals measured during the process. 4. System identification 4.1. Auto regressive model of correlation function The free decay (or IR) of lumped systems such as the system in Eq. (7) can be described by Auto Regressive (AR) models in discrete time domain[9]. Since the CF will be used to represent scaled IR in this work, similarly, an AR model is used to describe the CF in discrete time domain:

Ry ðiÞ ¼ A1 Ry ði  1Þ þ A2 Ry ði  2Þ þ . . . þ Ana Ry ði  naÞ

ð12Þ

where A1 . . . Ana are the coefficients of the AR model of order na. The AR model can also be expressed in companion form as follows:

Ud ði þ 1Þ ¼ UUd ðiÞ

ð13Þ

where Ud ðiÞ is composed of the CF matrices at na consecutive time steps:

h iT Ud ðiÞ ¼ RTy ði  na þ 1Þ    RTy ði  1Þ RTy ðiÞ

ð14Þ

and U is the companion matrix corresponding to the na order AR model in Eq. (12):

2

0

I

0

0

0

3

6 6 6 U¼6 6 6 4

0 .. . 0

0 .. . 0

I

0 .. . 0

0 .. . I

7 7 7 7 7 7 5

Ana

Ana1

   A2

A1

0

ð15Þ

Substituting one term from the modal expansion of the CF in Eq. (11) into Eq. (13) leads to the following eigenvalue equation, which shows that the eigenvalues of the companion matrix are also the discrete poles of the system in Eq. (7):

UZp ¼ lp Zp ;

h Zp ¼ BTp linaþ1 p

. . . BTp lpi1

BTp lip

iT

ð16Þ

In the next section, the coefficients of the AR model are determined from the CF of the vibrations measured during turning process, and then the poles of the system are determined as the eigenvalues of the corresponding companion matrix. 4.2. Time Domain Poly Reference (TDPR) identification In this method, the coefficients of the AR model of the CF are determined using the Least Squares approximation of the measured CF. The details of TDPR method are available in numerous sources such as [23,31], and a short description is also provided here. Because the machining forces, which include white noise excitation, are only applied on y1 and z1 , these DOFs are considered as the references. As a result, Ry ðiÞ in the AR model (Eq. (12)) are D  2 matrices consisting of the columns one and M þ 1 of the full CF matrix. Consequently, the coefficients in the AR model, Aj , are D  D matrices. By repeating Eq. (12) for i ¼ na þ 1; na þ 2; . . . ; np time steps, the following Hankel form matrix equation is obtained:

HT1 AT ¼ HT2

ð17Þ

where A contains the coefficients of the AR model, H1 is the Hankel block matrix consisting of the CF at discrete time steps, and H2 is the matrix of the CF at time steps na þ 1 to np:

A ¼ ½ Ana

Ana1

. . . A1 

ð18Þ

3 Ry ð2Þ ... Ry ðnp  naÞ R y ð 1Þ 6 R y ð 2Þ Ry ð3Þ Ry ðnp  ðna  1ÞÞ 7 7 6 7 H1 ¼ 6 . . .. 7 6 . . 5 4 . . . Ry ðnp  1Þ Ry ðnaÞ Ry ðna þ 1Þ

ð19Þ

H2 ¼ ½ Ry ðna þ 1Þ Ry ðna þ 2Þ    Ry ðnpÞ 

ð20Þ

2

320

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

The Least Squares (LS) approximation of A is obtained using the pseudo-inverse of the Hankel matrix, H1 :

b ¼ H2 Hþ A 1

ð21Þ

b is the estimate of A matrix and the superscript þ designates the pseudo-inverse of a matrix. Once the coefficients of where A the AR model are estimated from Eq. (21), the discrete poles of the system can be determined as the eigenvalues of the corresponding companion matrix. 4.3. Modified Least Squares Complex Exponential (LSCE) method In Section 2, the static component of the chip thickness was neglected because the cutting forces resulting from it do not influence the stability of the system. In practice, due to the presence of runout in the system, the static component of the chip thickness generates harmonic forces and consequently harmonic response at the spindle rotation frequency and its integer harmonics. The harmonic components usually dominate the response spectrum, and therefore they are identified as undamped poles when TDPR method is used. If the frequencies of the harmonic components are well separated from the modal frequencies of the system, one can simply neglect the poles with near zero damping at the harmonics of spindle rotation frequency, because they represent harmonic components of the response and not the system poles. However, if the frequency of the harmonic component is close to the modal frequencies of the system, the undamped nature of the harmonic response greatly affects the accuracy of the identification of the damped poles. To address this issue, Mohanty et al. [42] presented the modified Least Squares Complex Exponential method in which the known harmonic frequencies are included in the identification algorithm as undamped poles. The details of the identification method are available in [42], and a brief description of the method is also provided here. The discrete poles of the system, lp , satisfy the characteristic polynomial of the lumped system in Eq. (7):

b0 þ b1 l1p þ b2 l2p þ    þ b2P1 lp2P1 þ l2P p ¼ 0

ð22Þ

where b0 ; b1 ; . . . ; b2P , are the constant coefficients of the characteristic polynomial. It can be shown that the CF with respect to a single reference also satisfies the characteristic polynomial [31]:

b0 ry ðiÞ þ b1 ry ði þ 1Þ þ b2 ry ði þ 2Þ þ    þ b2P1 ry ði þ 2P  1Þ ¼ ry ði þ 2PÞ

ð23Þ

where ry ðiÞ is the column of the CF matrix Ry ðiÞ corresponding to the reference DOF. The harmonic components of the response can be treated as undamped poles, and therefore they also satisfy the characteristic polynomial (Eq. (22)). Repeating Eq. (23) for np consecutive time steps, and combining the resulting equations with the ones obtained by substituting the undamped poles (known harmonics) in the characteristic polynomial, leads to the following system of algebraic equations:

E1 b1 þ E2 b2 ¼ f 1

ð24Þ

and

E3 b1 þ E4 b2 ¼ f 2 where ½ b0



b2P1  ¼

ð25Þ h

T T b1 ; b2

i

, and the components of E1;2;3;4 and f 1;2 are discrete time CF values and discrete sinusoidal

functions of the harmonic terms. The components of these matrices are provided in Appendix A. While Eq. (24) enforces Eq. (23) at a set of time steps, Eq. (25) enforces the harmonic components to be undamped poles of the system. Therefore, in order to guarantee that the harmonic components are included as undamped poles, Eq. (25) must be exactly satisfied. The exact solution of Eq. (25) leads to the following:

b1 ¼ E1 3 ½f 2  E4 b2 

ð26Þ

Substituting Eq. (26) in Eq. (24) results in the following overdetermined system of equations:

  1 E2  E1 E1 3 E4 b2 ¼ f 1  E1 E3 f 2

ð27Þ

The least squares approximation of b2 can be obtained from Eq. (27), and along with the corresponding solution of b1 obtained from Eq. (26), they will form the complete solution of the coefficients of the characteristic polynomial, b0 ; . . . ; b2P1 . Poles of the system are then obtained as the roots of the characteristic polynomial. Note that the 2P roots of the polynomial include the considered harmonic terms and the damped poles of the system as pairs of complex conjugate poles. 5. Case study The experimental setup shown in Fig. 1 is used in this section to study the performance of the presented method in identifying the dynamics of turning systems. A turning tool holder with 38 mm overhang length and 19  19 mm cross section area is mounted on the table of a CNC milling machine. A cylindrical workpiece with 38 mm diameter and 51 mm overhang is mounted on the spindle. A triangular turning insert with zero degrees rake and approach angles is used to machine the

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

321

workpiece in the axial direction using 0.05 mm/rev feedrate and various spindle speeds and depth of cuts. Two PCB 352C22 accelerometers are installed on the tool in each of the feed (z) and normal (y) directions. In each direction, one of the accelerometers is installed close to the tool tip and the other one at 25 mm from the tool tip. A four channel NI9234 National Instrument Data Acquisition (DAQ) card is used to digitize the accelerometer signals at 16384 Hz sampling rate. The digitized signals are stored on a Personal Computer. Also, the sound pressure during machining is measured using a generic microphone. The direct and cross Frequency Response Functions (FRFs) between all four degrees of freedom were measured using impulse hammer tests and the resulting receptance (Hij ) plots are shown in Fig. 2. Modal Analysis was performed using Cutpro application [43] to identify the mass-normalized mode shapes, natural frequencies, and modal damping ratios of the two dominant modes at 1842 and 2445 Hz. The modal parameters of these two modes are shown in Table 1, and the FRFs that are constructed using the two identified modes are superimposed on the measured FRFs in Fig. 2. The two modes shown in Table 1 are used to obtain the stability diagrams of the system. The stability diagrams can be obtained by setting up a grid of spindle speed and depth of cut values and determining the stability of vibrations according to the eigenvalue of the transition matrix(Eq. (7)) associated with each point on the grid. The stability diagrams can also be obtained analytically by directly employing the FRFs in the frequency domain solution of the DDE as presented in [1]. The frequency domain solution was used in this paper and the resulting stability diagrams are shown in Fig. 3 5.1. Numerical simulations Before presenting the experimental results, a set of numerical simulations are conducted to verify that the dominant poles of the distributed system in Eq. (5) can be identified from its response to white noise excitation in discrete time domain, as explained in Section 4. The presented numerical simulations also demonstrate the evolution of the identified dominant poles of the system as it moves towards the border of stability. The block diagram shown in Fig. 4 outlines the algorithm used to numerically simulate the accelerations at the four measurement DOFs of the tool when arbitrary spindle speed and depth of cuts are used. The tool dynamics is represented by the two modes shown in Table 1. Assuming proportional damping, the mass-normalized mode shape vectors (W1 and W2 ) in Table 1 are used to convert the equation of motion of the tool from the physical space in Eq. (1) to modal space:

€ þ Cm q_ þ Km q ¼ f m ðt Þ; q y ¼ ½ W1 W2 q

Cm ¼ diagð2n1 x1 ; 2n2 x2 Þ   2 2 Km ¼ diag ð2pf 1 Þ ; ð2pf 2 Þ

f m ðtÞ ¼ ½ W1

ð28Þ

W2 T f ðtÞ

where q is the modal displacement vector, ni and xi ¼ 2pf i ; i ¼ 1::2, are the modal damping ratio and natural frequencies of the modes of the tool. Note that n denotes the modal damping ratio of the modes of the tool, and f denotes the damping of the poles of the closed loop system. The following state space form of Eq. (28) is used in the numerical simulations:

€ T ¼ A½ q q_ T þ Bf m ðt Þ ½ q_ q T q ¼ C½ q q_ 

Fig. 2. Direct and cross FRFs between the four DOFs (z1 ; z2 ; y1 ; y2 ) of the turning tool.

ð29Þ

322

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

Table 1 Modal parameters of the dominant modes of the turning tool. Mode

Frequency [Hz]

Damping ratio

Mode shape

r

fr

nr

1 2

1842 2445

0.021 0.009

Wr ¼ ½z1 ; z2 ; y1 ; y2 T [2.81,1.73,0.94,0.71]T [0.92,0.66,2.55,1.39]T

Fig. 3. stability diagrams of the turning system, and the damping ratio of the system identified from vibration signals obtained from numerical simulations (a), and from the vibration signals measured during turning operations (b); Circles and crosses indicate stable and unstable vibrations, respectively, and the shading of circles are proportional to the damping ratio of the most lightly damped pole; Frequency Spectrum of the sound signals recorded during turning at (c) 7500 rev/min spindle speed and 3.05 mm depth of cut, (d) 7500 rev/min and 2.54 mm depth of cut.

where A; B, and C are the system matrices in modal space:

A¼

0

I

Km Cm I ¼ I22 ; 0 ¼ 022

;B ¼



0 I

; C ¼ ½ I 0 ;

ð30Þ

Using orthogonal to oblique transformation method [1], the cutting force coefficients of the Al 6061 workpiece are obtained K t ¼ 635 and K f ¼ 159 MPa. The closed loop system shown in Fig. 4 is assumed to be initially at rest, and it is excited by normally distributed white noise at every 0.01 ms simulation time step. The resulting accelerations at the four DOFs are obtained by numerically integrating the modal state space equation (Eq. (29)). The system response is simulated for 10 s, and each simulation is repeated two times to obtain two data sets. The data set obtained from the first simulation is used for the identification of the system poles, and the data set from the second simulation is used for the cross-validation of the model identified using the first data set. For example, Fig. 5 shows the Singular Values of the Power Spectral Density (PSD) matrices of the responses used for the identification (data set 1) of the system poles when the spindle speed and depth of cut are (a) 7500 rev/min and 2.03 mm, (b) 7500 rev/min and 2.54 mm, (c) 7300 rev/min and 3.05 mm, and (d) 7300 rev/ min and 4.06 mm. According to the stability diagrams shown in Fig. 3, at 7500 rev/min the system has the lowest stable depth of cut, and 7300 rev/min coincides with one of the stability pockets of the system.Only the first two (out of four) singular values of the responses are shown in Fig. 5; the two smallest singular values are below 100 dB.At each frequency line,

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

323

Fig. 4. Block Diagram of the numerical simulations of self-excited vibrations during turning.

the number of dominant singular values of the PSD determines the number of dominant poles contributing to the response of the system at that frequency [44]. Consequently, the singular values of the PSD matrix plotted in the frequency range of interest exhibits peaks at the frequencies of the dominant poles of the system. The singular value plot of the PSD is therefore commonly used as the Complex Modal Indicator Factor (CMIF) in OMA algorithms[27]. Note that the PSD matrix and its singular values are not used in the identification methods in this work; they are only used as CMIF showing the number (and approximate frequencies) of the dominant poles of the system. The CF of the simulated acceleration signals is computed at np ¼ 1500 time steps using direct method[45]. The resulting CF is then used in the TDPR identification method discussed in Section 4.2 to identify the frequency and damping ratio of the dominant poles of the system. In order to determine the suitable AR model order, na, stabilization diagrams and cross-validation methods are used. The stabilization diagram of each simulated case is obtained by increasing the model order from na ¼ 2, at unit increments, to na ¼ 40. After each model order increment, the frequencies and damping ratios of the identified poles are compared with their corresponding values at the preceding model order, and the poles with a positive damping ratio that vary less than 0.2% in their frequency are identified as stable poles and shown with a circle on the PSD diagrams. Therefore, a streak of circles on the stabilization diagram indicates an identified system pole. The shading of the circles is proportional to their damping ratio. Darker circles are the poles with a higher damping ratio and the light circles have a lower damping ratio. The stabilization diagram of the simulated vibrations at 7500 rev/min spindle speed and 2.03 mm depth of cut is shown in part (a) of Fig. 5. At this point, two stable poles are identified, a highly damped pole at 1832 Hz and a lightly damped pole at 2465 Hz; the former is associated with mode 1 in Table 1 and the latter is associated with mode 2. When the depth of cut increases to 2.54 mm at the same spindle speed (7500 rev/min), as shown in part (b) of the figure, two stable poles are identified at 1835 Hz with a high damping (indicated by dark circles) and 2468 Hz with very low damping (nearly blank circles). By increasing the depth of cut to 3.05 mm at 7500 rev/min, vibration signals become unstable.The poles of the unstable system cannot be identified using TDPR method, because the post-chatter system is strongly nonlinear [46] and the resulting unstable vibration signals are non-stationary. Unstable points are indicated by crosses in Fig. 3 (a). At 7300 rev/min, the frequency of the second mode at 2445 Hz synchronizes with the spindle frequency (2445 Hz ¼ 20  7300=60 HzÞ, producing twenty complete waves in every spindle revolution[47,48]. Because of the synchronization of the spindle revolution frequency and the vibration frequency, a higher stability limit is obtained at this speed. As shown in parts (c) and (d) of Fig. 5, the PSD of the simulated signals at 7300 rev/min include two closely spaced peaks at around the second mode of vibrations (2445 Hz). The TDPR identification method identifies the frequencies and damping ratio of both of the peaks around the second mode, as well as the frequency and damping of the peaks around the first mode. In both cases in Figs. 5 (c) and (d), each tool mode shown in Table 1 splits into a lightly damped pole with a higher participation (converges at a lower model order) and a pole with a higher damping and a lower participation. At 3.05 mm depth of cut, the pole with the smallest damping ratio is at 2499 Hz; at 4.06 mm, the pole with the lowest damping ratio is at 2506 Hz. By comparing the shadings of the circles at these frequencies, it becomes clear that the lowest damping ratio decreases by increasing the depth of cut from 3.05 to 4.06 mm. At 4.32 mm depth of cut, the tool vibrations become unstable in agreement with the analytically obtained stability diagrams. The overall reduction of the damping ratios of the dominant poles when the depth of cut approaches the stability limit is clear in Fig. 5, and a similar trend is observed in the simulations performed at other spindle speeds. As also shown in the stabilization diagrams in Fig. 5, the frequencies of the stable poles do not change significantly (< 0:2%) by changing the order of the AR model. However, unlike the frequencies, the damping ratio of the identified poles vary by changing the model order. For instance, the variation of the damping ratio of the second pole in Fig. 5(a) when various model orders are used is shown in Fig. 6. Also shown in this figure is the resulting estimation errors when models with various orders are used. _

In order to obtain the estimation error, the identified poles (lp ) are employed in Eq. (11) to reconstruct the CF of the vibration signals in the second dataset:

Ry ðiÞ ¼

2P X _ Bp lip p¼1

ð31Þ

324

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

Fig. 5. Stabilization diagrams of the TDPR identification performed on the vibration signals obtained from numerical simulations of the turning process at (a)7500 rev/min and 2.03 mm depth of cut, (b) 7500 rev/min and 2.54 mm depth of cut, (c)7300 rev/min and 3.05 mm depth of cut, and (d) 7300 rev/min and 4.06 mm depth of cut.

Note that Ry in Eq. (31) are the CF of the vibration signals in the second dataset, and P is the number of identified stable poles. For example, in Fig. 5(a) to estimate the error when na ¼ 10, only two stable poles are included and therefore P ¼ 2. Unknown constant matrices, Bp , are composed of deflection shapes associated with the pole number p. These unknown matrices are determined by repeating Eq. (31) for the first i ¼ 0 . . . np time steps and then using the following Least Squares approximation:

h_ Bp

iT

_

2

   B2P _

l01

6_ 6 l1 6 1 Z¼6 6 .. 4 . _

_

l02 _ l12 .. .

_

T ¼ Zþ ½ Ry ð0Þ    Ry ðnpÞ  3 _    l02P 7 _    l12P 7 7 7 .. .. 7 . . 5

ð32Þ

_

lnp lnp    lnp 1 2 2P

_

_

where Bp stands for the least squares approximation of Bp . Reconstructed CF, Ry ðiÞ, are obtained by replacing Bp in Eq. (31) by

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

325

  Fig. 6. Estimation error and damping ratio fp of the dominant pole obtained from TDPR identification using the AR models with various orders; vibration signals are obtained from the simulation of the turning process at 7500 rev/min and 2.03 mm depth of cut.

_

its approximation, Bp : _

Ry ðiÞ ¼

2P _ X _ Bp lip

ð33Þ

p¼1

Estimation error at each time step is defined as the Euclidean norm of the difference between the reconstructed and measured (simulated) CF:

_



eðiÞ ¼ Ry ðiÞ  Ry ðiÞ

ð34Þ

The norm of the estimation error (e) at 7500 rev/min and 2.03 mm is shown in Fig. 6. According to this figure, estimation error is minimum when the AR model with order na ¼ 25 is used. At this model order, the damping ratio of the second pole at 2466 Hz is estimated at 0.16% (f2 ¼ 0:0016). Similar cross-validation method is used to determine the optimum model order and the corresponding damping ratio in all of the other simulated cases. The frequency and damping ratio of the poles identified from the simulated vibrations at 7500 rev/min and depth of cut 0 < a < 2:54mm are shown in Table 2. Similar analysis is performed at a set of spindle speed and depth of cut values and the results are demonstrated in Fig. 3(a). Simulations at the spindle speed and depth of cut values that result in stable vibrations are represented by a circle in Fig. 3; simulations at points that result in unstable vibrations are shown with a cross. At each simulation point, the shading of the corresponding circle is proportional to the damping ratio of the identified pole with the lowest damping. As seen in this figure, at all of the spindle speeds, as the depth of cut approaches the stability limit, the lowest damping ratio approaches zero. The variation of the smallest damping ratio and the corresponding frequency, as the depth of cut increases from zero to 2.54 mm at 7500 rev/min, is shown in Fig. 7. As shown in this figure, when the depth of cut is below 0.5 mm, the damping ratio is nearly equal the modal damping of the tool’s second mode (2445 Hz, and 0.98% damping). As the depth of cut approaches the stability limit, the damping ratio decreases to 0.06% and the frequency increases to 2468 Hz. 5.2. Experimental results The vibration signals measured during the conducted turning experiments are used in this section to identify the poles of the system and thereby estimate its stability margin. The singular values of the PSD of the measured acceleration signals at 7500 rev/min and 2.05 mm depth of cut are shown in Fig. 8. Only the two largest singular values are shown; the two smallest singular values are below 100 dB. The integer harmonics of the spindle revolution frequency (125 Hz) are marked by dashed vertical lines. As shown in this figure, sharp peaks are observed at the harmonics of the spindle frequency due to the strong presence of harmonic excitation. Also, a peak is observed next to the 20th harmonic of the spindle frequency, which indicates a damped pole at this frequency. No peaks are observed near the first mode at 1842 Hz, because this mode has a higher damping and therefore is suppressed by the harmonic terms, which are virtually undamped with a higher participation in the response. Before applying TDPR identification method on the measured response, detrending FFT filter [23] is applied on the measured signals. Also, as shown in Fig. 8, unlike the simulated data in Section (5.1), the measured signals include many structural poles and harmonic components, and therefore a very high order AR model is required to identify all of these poles simultaneously. In order to avoid the application of high order AR models (which may cause numerical issues),

326

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

Table 2 Poles obtained from TDPR identification of vibration signals simulated at 7500 rev/min spindle speed as the depth of cut increases from zero to the stability limit. Model order na ¼ 25 was used. Pole 1 ðp ¼ 1Þ a[mm] 0 0.51 1.0 1.52 2.03 2.54





Pole 2 ðp ¼ 2Þ

Freq. xp [Hz]

  Damping fp [%]





Freq. xp [Hz]

  Damping fp [%]

1845 1845 1843 1842 1841 1840

2.16 2.16 2.1 2.03 1.95 1.84

2455 2455 2463 2465 2466 2468

0.98 0.98 0.52 0.3 0.16 0.04

  Fig. 7. Frequency and damping ratio fp of the dominant pole obtained from the TDPR identification of the vibrations simulated at 7500 rpm and various depth of cuts.

a bandpass FFT filter is applied around the peak at 2476 Hz. The flat range of the bandpass filter is 25% of the spindle frequency, and the total width of its pass band is 75% of the spindle frequency. The largest singular value of the PSD of the measured acceleration signal after the application of bandpass filter is also shown in Fig. 8. The CF of the filtered signals are computed using direct method [45] at np ¼ 1500 time steps. The computed CF is then used in the TDPR method to identify the poles of the system that are located within the pass band of the bandpass filter. The resulting stabilization diagram of the identification at 7500 rev/min and 2.05 mm depth of cut is shown in Fig. 8(a). Similar to the method used in Section 5.1, the stabilization diagram is developed by increasing the model order at unit increments. The identified poles at each increment are considered stable and depicted by a circle in the stabilization diagram if their frequency is within the 0.2% of their corresponding frequency in the previous model order. The color shading of the circles are proportional to their corresponding damping ratio. Darker circles indicate higher damping. As shown in Fig. 8 (a), the harmonic response peak at 2500 Hz is identified as a pole with near zero damping, but the system pole next to the harmonic response peak does not converge to a stable frequency. Although TDPR method efficiently identified the poles of the system in numerical simulations, where the harmonic component of the response was neglected, it is not effective in identifying the modes from experimentally measured signals due to the presence of strong harmonic components. To compensate for the undamped effect of the harmonic components in the computed CF, the modified LSCE method is used. Three harmonic terms at 19, 20, and 21 times the spindle frequency are considered as undamped poles of the system in the CF. The CF with respect to excitations at z1 DOF are considered in identification, therefore only the first column of the CF matrix is used. The resulting stabilization diagram at 7500 rev/min and 2.05 mm depth of cut is shown in Fig. 8(b). Using the modified LSCE method, the three harmonic terms are identified as poles with zero damping and an additional damped system pole is identified at 2476 Hz. Similar to the numerical simulations in Section 5.1, cross validation method is used to determine the suitable model order (2P), except that in this section the CF with respect to y1 DOF (i.e. the third column of CF matrix) is used for cross-validation, instead of a second dataset. The resulting estimation error at each model order is shown in part (c) of Fig. 8. Estimation error is minimum at na ¼ 11 and therefore this number is used as the optimum model order. The corresponding damping ratio of the pole at 2476 Hz is estimated at 0.13%. Variation of the estimated damping ratio at various model orders is also shown in Fig. 8 (c). The PSD of the acceleration signals measured at (a) 7500 rev/min and 2.54 rev/min, (b) 7300 rev/min and 3.05 mm, (c) at 7150 rev/min and 2.29 mm, and (d) 7150 rev/min and 2.54 mm are shown in Fig. 9. Similar to Fig. 8, detrending and band pass filters around the observed non–harmonic peaks are applied before using modified LSCE identification method. In part

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

327

Fig. 8. (a) Stabilization diagram of TDPR identification performed on the vibration signals measured at 7500 rev/min and 2.03 mm depth of cut, (b) Stabilization diagram of modified LSCE identification performed on vibrations measured at 7500 rev/min and 2.03 mm depth of cut, and (c) the estimation error and the damping ratio of the dominant mode identified from the AR model (with various orders) of the turning system at 7500 rev/min and 2.03 mm depth of cut.

(a) of the figure, the stabilization diagram resulting from identification at 7500 rev/min and 2.54 mm depth of cut shows a stable pole at 2477 Hz with 0.066% damping (na ¼ 11 model order). By increasing the depth of cut to 3.05 mm at 7500 rev/ min, vibraions becomes unstable, which is also predicted by the stability diagrams in Fig. 3. At 7300 rev/min, where the 20th

328

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

harmonic of the spindle frequency coincides with the second vibration mode, the cut is stable at 3.05 mm but becomes unstable at 3.3 mm depth of cut. The higher stability at this spindle speed was also observed in numerical simulations, and it was predicted by stability diagrams. However, the stability limit (3.3 mm) is smaller than the simulation result (4.06 mm). This discrepancy is attributed to unmodeled process disturbances. As shown in part (b) of Fig. 9, no non–harmonic peaks are observed even at the highest stable depth of cut (3.05 mm) at 7300 rev/min, and therefore system identification cannot be performed at this speed. When the spindle speed is within one of the stability pockets of the stability diagram, the frequency of the dominant system pole coincides with one of the harmonic components of the signal and thus

Fig. 9. Stabilization diagram of modified LSCE identification performed on the vibration signals measured at (a) 7500 rev/min and 3.05 mm depth of cut, (b) 7300 rev/min and 3.05 mm depth of cut, (c) 7150 rev/min and 2.29 mm depth of cut, and (d) 7150 rev/min and 3.05 mm depth of cut.

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

329

it cannot be identified. To examine the consistency of the identification results, poles of the system at 7150 rev/min and 2.29 mm and 2.54 mm depth of cut are identified and the results are shown in parts (c) and (d) of Fig. 9, respectively. As shown in these figures, at 2.29 mm, a stable pole is identified at 2474 Hz with 0.12% damping ratio, and at 2.54 mm, the damping of this pole reduces to 0.076% and its frequency increases to 2477 Hz. At 3.05 mm depth of cut vibrations become unstable, similar to 7500 rev/min. Similar analysis is performed on the vibration signals recorded at a set of spindle speed and depth of cut values and the results are summarized in Fig. 3(b). Also, the stability diagrams are duplicated in this figure to enable a comparison between experimental results, analytically obtained stability limits, and numerical simulations presented in part (a) of the figure. The black circles in part (b) of the figure are the points at which the damping ratio of the system poles are not low enough to overcome the dominance of the harmonic components in their vicinity and therefore no non–harmonic peaks were observed in the PSD. Because the damping of the system poles cannot be identified at these points, the modal damping of the second tool mode in idle condition (0.9%) is assigned. The circles with light shading are the points at which the damping and frequency of the dominant system pole were successfully identified using modified LSCE method. Darker shading indicates a higher damping. Crosses indicate points at which chatter was detected both in the measured vibration signals and the recorded sound signals. According to Fig. 3, as the depth of cut approaches the stability limit, the damping of the dominant system pole becomes small enough to overcome the harmonic terms in the vibration signal. As a result, the modified LSCE method successfully identifies the damping and frequency of the dominant pole of the system. The parameters of the identified poles are in a close agreement with the ones obtained in numerical simulations shown in part (a) of the figure. At lower depth of cuts, where the cut is highly stable, no non–harmonic peaks are observed in the PSD plot and therefore the presented method only identifies the harmonic terms of the response. In summary, the process that was used to identify the poles of the machining system can be described as follows: 1. Record vibration (e.g. acceleration) signals at N points in the feed (Z) and M points in cutting (Y) directions. 2. Apply de-trending filter on the measured signals to remove their DC offset (usually caused by the thermal drift of electronic components) 3. Plot the singular values of the PSD matrix of the detrended signals, as shown in Figs. 5, 8, and 9. 4. Apply a bandpass filter around the dominant peak in the PSD diagram (obtained in the previous step) that is not a harmonic of the spindle rotation frequency 5. Compute the CF of the filtered signals, for example, using direct method. Either Z1 or Y1 can be used as the reference DOF. 6. Use the resulting CF and the two closest spindle harmonic frequencies to the dominant peak in the modified LSCE algorithm, described in Section 4.3, to identify the poles of the system. The suitable model order can be determined by performing the cross-validation process described in Section 5.2. 7. Neglect the numerical (or noise) poles with unrealistically high or negative damping ratio. Among the remaining poles, the one with the lowest positive damping ratio is the critical pole that is prone to becoming unstable. The frequency spectra of the recorded sound signals are commonly used to detect chatter in machining operations. Parts (c) and (d) of Fig. 3 show the frequency spectra of the sound signals recorded at 7500 rev/min and 2.54 mm (stable point) and at 7500 rev/min and 3.05 mm (unstable point). According to these figures, at the stable point, no peaks close to the frequency of the tool modes (1842 Hz and 2445 Hz) are detected in the sound signal. When the depth of cut exceeds the stability limit, a sharp peak appears at 2471 Hz, indicating chatter. While the spectrum of the sound signal detects chatter only after vibrations become unstable, the presented method can be used to monitor the declining trend of the damping of the dominant pole of the system to indicate the loss of stability while the cut is still stable.

6. Conclusions Regenerative vibrations in turning processes were modeled using a finite dimensional lumped system in discrete time domain, and the poles of the resulting system were identified using Operational Modal Analysis of the vibrations measured during the process. The presented method was used to identify the damping ratio of the dominant poles of the turning system, which decrease as the system moves towards the border of stability. The variation of the damping ratio of the dominant system pole can be used as a monitoring parameter to prevent the occurrence of chatter before the system becomes unstable. The performance of the presented method was studied using numerical simulations and an experimental case study. The method was shown to successfully identify the structural poles of the system when the cut is close to the border of stability. When the cut is highly stable or when the harmonics of spindle rotation frequency coincides with one of the modal frequencies of the system, the presented method cannot identify the structural poles. However, monitoring the process for chatter occurrence is not critical in those highly stable points. In order to use the presented approach as online chatter monitoring system, the OMA algorithm needs to be automated and its computational time needs to be reduced. Automated OMA is active research and substantial progress have been made in the application of automated OMA for structural health monitoring [49]. While modification of the presented approach for online monitoring of chatter remains a subject for future studies, the presented approach can be used in its current offline

330

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

form to quantify the stability of machining processes. This capability is particularly useful in machining chatter research where the accuracy of the developed stability lobe diagrams is verified according to the stable and unstable cutting conditions that are identified through machining experiments. In those machining experiments, usually, sound, vibration, and force signals or surface roughness analysis are used to determine if the cut is stable or unstable. The border between stable and unstable cutting parameters determined from the machining experiments are then compared against the stability lobe diagrams. Although such tests provide an approximation of the border between stable and unstable machining parameters, they don’t provide any information about the dynamics of the system when it is stable, nor do they provide any numerical metric of the level of stability of the system. The dominant closed-loop system poles that are identified using the presented approach can be used not only to determine the border between stable and unstable machining parameters, but also to determine the damping and frequency of the dominant poles that governs the dynamics of the system when it is stable. Also the damping of the dominant pole can be used as a numerical metric to evaluate the stability of the system. Acknowledgements This research was financially supported by the National Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grant program and Collaborative Research and Development (CRD) Grant No. CRDPJ/490628. Also, the authors thank Dr. Martin Juul for contributions to the initial stages of the project. Appendix A A.1. Matrices of modified LSCE method Consider m harmonic terms with x1 . . . xm rotational frequencies in the vibration signals measured at t s sampling periods. The components of E1;2;3;4 and f 1;2 in Eqs. (24) and (25) are expressed in terms of the CF of the measured signals, as follows:

3  ry ð2m  1Þ ry ð0Þ 7 .. .. 7 .  . 5 ry ðnp  1Þ    ry ðnp þ 2m  2Þ 2 ry ð2mÞ  ry ð2P  1Þ 6 . .. 6 E2 ¼ 4 ..  . 2

6 E1 ¼ 6 4

ry ðnp þ 2m  1Þ    ry ðnp þ 2P  2Þ 3 0    sinðx1 ð2m  1Þts Þ 6 1    cosðx ð2m  1Þt Þ 7 6 1 s 7 7 6 7 6 .. . . E3 ¼ 6 .    7 . 7 6 7 6 4 0    sinðxm ð2m  1Þt s Þ 5

3 7 7 5

2

1    cosðxm ð2m  1Þt s Þ 3 sinðx1 2mt s Þ    sinðx1 ð2P  1Þt s Þ 6 cosðx 2mt Þ    cosðx ð2P  1Þt Þ 7 6 1 s 1 s 7 7 6 7 6 . . . . E4 ¼ 6 7 .    . 7 6 7 6 ð x 2mt Þ    sin ð x ð 2P  1 Þt Þ sin 4 m s m s 5 cosðxm 2mts Þ    cosðxm ð2P  1Þts Þ 3 2 sinðx1 2Pts Þ 2 3 6 cosðx 2Pt Þ 7 ry ð2PÞ 6 1 s 7 7 6 6 7 7 6 . . 6 7 . . f1 ¼ 4 7 . . 5; f 2 ¼ 6 7 6 7 6 ry ðnp þ 2P  1Þ 4 sinðxm 2Pt s Þ 5 cosðxm 2Pts Þ 2

ð35Þ

A.2. Transition matrix Eq. (6) represents the dynamics of the delayed system during one discrete time interval t i 6 t < t iþ1 by assuming that the delayed term is constant Kc ayðt ir Þ within each interval. The complete solution of the ODE in Eq. (6), including the homogeneous and particular solutions, is as follows:

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

 1 Kc ay yðt Þ ¼ WKc1 þ WK c2 þ K ir K ¼ diagðek1 t ; . . . ; ekN t Þ; W ¼ ½ W1

k1 ¼ f1 x1 þ jx1

qffiffiffiffiffiffiffiffiffiffiffiffiffi 1  f21

331

ð36Þ

   WN 

where Wl ; l ¼ 1::N are the mass normalized mode shapes of the tool and xl and fl are the corresponding modal frequency and damping ratio. The total number of modes considered in the solution are N, and proportional damping is assumed. The constant vectors, c1 and c2 are determined by applying the boundary values of yðtÞ and its time derivative at t ¼ t i :

      1 1 W1 y_ i  WK  1 W1 y  K1 Kc ay 2  K c2 ¼ K1 K i ir i        1 2  K  1 1 W1 y_ i  WK  1 W1 y  K1 Kc ay yi  W K  K1 Kc ayir c1 ¼ K1 i ir i W

ð37Þ

Ki ¼ Kðt i Þ 1 ¼ K   ¼ diagðk1 ;    ; kN Þ. The displacement vector and its derivative at the next time step, tiþ1 , are obtained by subwhere K 2 stituting the constants, c1 and c2 , from Eq. (37) in Eq. (36):

yiþ1 ¼ P11 yi þ P12 y_ i þ R1 yir y_ iþ1 ¼ P21 yi þ P22 y_ i þ R2 yir

ð38Þ

where

    1 W1 ;   2  K K  P11 ¼ W Kh K h 1 K2  K1     1 1 W1 ; 2  K P12 ¼ W Kh  Kh K     1 1 W1 ; 2  K  2 Kh  K K P21 ¼ WK1 K h     1 1 W1 ; 2  K  2  Kh K 1 K P22 ¼ W Kh K       1 1 W1 K1 Kc a; 2  K  1  Kh K 2 K R1 ¼ W 1 þ Kh K       1K  1 1 W1 K1 Kc a 2  K  2 K  Kh K R2 ¼ W 1 þ K h

ð39Þ

According to Eq. (39), the components of the transition matrix,G, in Eq. (7) are obtained as follows:

2

P11

6P 6 21 6 6 I 6 G¼6 6 0 6 . 6 . 4 . 0

P12 P22 0 0 0

0 0 . . . 0 R1

3

0 0 . . . 0 R2 7 7 7 0 0 ... 0 0 7 7 I 0 ... 0 0 7 7 .. . 7 7 . . . 5 0 0 ...

I

ð40Þ

0

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

Y. Altintas, Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design, Cambridge University Press, 2012. J. Tlusty, M. Polacek, The stability of machine tools against self-excited vibrations in machining, 1963, Proceedings of the ASME International. S.A. Tobias, Machine-tool Vibration, J. Wiley, 1965. M. Siddhpura, R. Paurobally, A review of chatter vibration research in turning, Int. J. Mach. Tools Manuf. 61 (2012) 27–47. G. Quintana, J. Ciurana, Chatter in machining processes: a review, Int. J. Mach. Tools Manuf. 51 (5) (2011) 363–376. E. Budak, Y. Altintas, Analytical prediction of chatter stability in milling-part i: general formulation, J. Dyn. Syst., Measure., Control 120 (1) (1998) 22– 30. B. Yanushevsky, A new theoretical approach for the prediction of machine tool chatter in milling, J. Eng. Ind. 115 (1993) 1. D. Bachrathy, G. Stepan, Improved prediction of stability lobes with extended multi frequency solution, CIRP Ann.-Manuf. Technol. 62 (1) (2013) 411– 414. C.-T. Chen, Linear System Theory and Design, Oxford University Press Inc, 1998. T. Insperger, B.P. Mann, G. Stépán, P.V. Bayly, Stability of up-milling and down-milling, part 1: alternative analytical methods, Int. J. Mach. Tools Manuf. 43 (1) (2003) 25–34. P. Bayly, J. Halley, B.P. Mann, M. Davies, Stability of interrupted cutting by temporal finite element analysis, J. Manuf. Sci. Eng. 125 (2) (2003) 220–225. E.A. Butcher, H. Ma, E. Bueler, V. Averina, Z. Szabo, Stability of linear time-periodic delay-differential equations via chebyshev polynomials, Int. J. Numer. Meth. Eng. 59 (7) (2004) 895–922. T. Insperger, G. Stépán, Semi-discretization method for delayed systems, Int. J. Numer. Methods Eng. 55 (5) (2002) 503–518. Y. Ding, L. Zhu, X. Zhang, H. Ding, A full-discretization method for prediction of milling stability, Int. J. Mach. Tools Manuf. 50 (5) (2010) 502–509. R. Teti, K. Jemielniak, G. O’Donnell, D. Dornfeld, Advanced monitoring of machining operations, CIRP Ann.-Manuf. Technol. 59 (2) (2010) 717–739. E. Kuljanic, G. Totis, M. Sortino, Development of an intelligent multisensor chatter detection system in milling, Mech. Syst. Signal Process. 23 (5) (2009) 1704–1718. T.L. Schmitz, Chatter recognition by a statistical evaluation of the synchronously sampled audio signal, J. Sound Vib. 262 (3) (2003) 721–730. A. Honeycutt, T.L. Schmitz, A new metric for automated stability identification in time domain milling simulation, J. Manuf. Sci. Eng. 138 (7) (2016) 074501. T. Delio, J. Tlusty, S. Smith, Use of audio signals for chatter detection and control, J. Eng. Ind. 114 (2) (1992) 146–157.

332

S. Kim, K. Ahmadi / Mechanical Systems and Signal Processing 130 (2019) 315–332

[20] Y. Altintas, P.K. Chan, In-process detection and suppression of chatter in milling, Int. J. Mach. Tools Manuf. 32 (3) (1992) 329–347. [21] A.K. Kiss, D. Hajdu, D. Bachrathy, G. Stepan, Operational stability prediction in milling based on impact tests, Mech. Syst. Signal Process. 103 (2018) 327–339. [22] A. Brandt, A signal processing framework for operational modal analysis in time and frequency domain, Mech. Syst. Signal Process. 115 (2019) 380– 393. [23] R. Brincker, C. Ventura, Introduction to Operational Modal Analysis, Wiley, 2015. [24] I. Zaghbani, V. Songmene, Estimation of machine-tool dynamic parameters during machining operation through operational modal analysis, Int. J. Mach. Tools Manuf. 49 (12–13) (2009) 947–957. [25] K. Ahmadi, Output-only modal analysis of micromilling systems, J. Micro Nano-Manuf. 6 (1) (2018) 011006. [26] M. Wan, J. Feng, Y.-C. Ma, W.-H. Zhang, Identification of milling process damping using operational modal analysis, Int. J. Mach. Tools Manuf. 122 (2017) 120–131. [27] B. Peeters, G. De Roeck, Stochastic system identification for operational modal analysis: a review, J. Dyn. Syst. Meas. Contr. 123 (4) (2001) 659–667. [28] E. Reynders, System identification methods for (operational) modal analysis: review and comparison, Arch. Comput. Methods Eng. 19 (1) (2012) 51– 124. [29] D.L. Brown, R.J. Allemang, R. Zimmerman, M. Mergeay, Parameter estimation techniques for modal analysis, SAE Trans. (1979) 828–846. [30] D.J. Ewins, Modal Testing: Theory and Practice, vol. 15, Research Studies Press, Letchworth, 1984. [31] N. Maia, J. Silva, Theoretical and Experimental Modal Analysis, Engineering Dynamics Series, Research Studies Press, 1997. [32] Z. Dombovari, Dominant modal decomposition method, J. Sound Vib. 392 (2017) 56–69. [33] P. Van Overschee, B. De Moor, Choice of state-space basis in combined deterministic-stochastic subspace identification, Automatica 31 (12) (1995) 1877–1883. [34] H. Vold, J. Kundrat, G.T. Rocklin, R. Russell, A multi-input modal estimation algorithm for mini-computers, SAE Trans. (1982) 815–821. [35] R.J. Allemang, D.L. Brown, Experimental modal analysis and dynamic component synthesis. volume 3. modal parameter estimation, Tech. rep, Cincinnati Univ OH Dept of Mechanical and Industrial Engineering (1987). [36] R.R. Craig Jr, A.J. Kurdila, H.M. Kim, State-space formulation of multi-shaker modal analysis. [37] F. Lembregts, J. Leuridan, H. Van Brussel, Frequency domain direct parameter identification for modal analysis: state space formulation, Mech. Syst. Signal Process. 4 (1) (1990) 65–75. [38] H. Van der Auweraer, R. Snoeys, J. Leuridan, A global frequency domain modal parameter estimation technique for mini-computers, J. Vib. Acoust.Trans. ASME. [39] P. Verboven, P. Guillaume, B. Cauberghe, E. Parloo, S. Vanlanduit, Frequency-domain generalized total least-squares identification for modal analysis, J. Sound Vib. 278 (1–2) (2004) 21–38. [40] J. Busturia, J. Gimenez, Multiexcitation multiresponse non-linear least squares algorithm, Proc. 10-th ISMA. [41] T. Insperger, G. Stépán, Semi-Discretization for Time-Delay Systems: Stability and Engineering Applications, 178, Springer Science & Business Media, 2011. [42] P. Mohanty, D.J. Rixen, Operational modal analysis in the presence of harmonic excitation, J. Sound Vib. 270 (1) (2004) 93–109. [43] Cutpro, Manufacturing automation inc., www.malinc.com. [44] R. Brincker, L. Zhang, P. Andersen, Modal identification from ambient responses using frequency domain decomposition, in: Proc. of the 18thInternational Modal Analysis Conference (IMAC), San Antonio, Texas, 2000. [45] J.S. Bendat, A.G. Piersol, Engineering Applications of Correlation and Spectral Analysis, Wiley-Interscience, New York, 1980, p. 315. [46] J. Tlusty, F. Ismail, Basic non-linearity in machining chatter, CIRP Ann. 30 (1) (1981) 299–304. [47] A.K. Kiss, D. Bachrathy, G. Stépán, Laser scanned patterns of machined surfaces, Procedia CIRP 77 (2018) 355–358. [48] Y. Altintas, M. Weck, Chatter stability of metal cutting and grinding, CIRP Ann. 53 (2) (2004) 619–642. [49] E. Reynders, J. Houbrechts, G. De Roeck, Fully automated (operational) modal analysis, Mech. Syst. Signal Process. 29 (2012) 228–250.