Nonlinear Dynamic Analysis of Turning Machine with and without Nonlinear Suspension

We performed a systematic analysis of the dynamic behavior of the cutting process of a machine tool with an unbalanced force induced by the mass eccentricity of the workpiece, a nonlinear cutting force, and a nonlinear suspension effect. Phase diagrams, power spectra, bifurcation diagrams, and a Poincaré section were applied to identify the dynamic motion. The simulation results show that nonperiodic dynamic responses are very prevalent in the cutting process of the machine tool. A stability analysis and modeling with multiple-degree-of-freedom coupling systems were also performed in this study. The results provide an understanding of the operating conditions under which undesirable dynamic motion takes place in this type of system and therefore serve as a useful reference for engineers designing and controlling such systems for precision manufacturing. Moreover, the results of this research will help clarify the dynamic characteristics of cutting machine tools and provide relevant and useful parameters to build or design better sensors.


Introduction
The complicated interaction mechanisms among the tool, workpiece, and chips in a cutting process make them difficult to analyze and to provide a comprehensive solution. Many significant and groundbreaking investigations on machine cutting have been reported. The dynamics of cutting has been analyzed on the basis of the model of Hastings et al. (1) in many studies. Grabec (2)(3)(4) presented a series of papers discussing the chaotic dynamic responses occurring in cutting machines and also found some mechanisms involving chaos in the cutting process. Hanna and Tobias (5,6) analyzed a machine tool structure in terms of the nonlinear stiffness effect, hysteretic damping effect, and cutting force (the variation of the chip thickness was modeled as a cubic polynomial) and found many significant types of nonperiodic motion. Wu and Liu (7) modeled a cutting system using the friction and empirical ploughing forces and performed experiments to verify their simulation results. Johnson and Moon (8) presented a onedegree-of-freedom model to simulate the cutting behaviors in the machining process and analyzed both the prechatter and chatter vibrations. Nosyreva and Molinari (9) also simulated a cutting model considering the velocity-dependent friction and ploughing force, and the experimental results were in accordance with the simulation results. Jemielniak and Widota (10) carried out a numerical simulation of nonlinear chatter vibration in turning machines. Other noteworthy studies (11)(12)(13) also found widespread nonperiodic motion or the so-called chaotic motion when simulating cutting machines. Therefore, it is known that the cutting process is a highly nonlinear phenomenon and that linearization or simplification when analyzing the process may result in large simulation errors.
In any automated process, sensors and their signal interpretation play an important role. A variety of sensors are used in systems for monitoring cutting tool conditions, such as power, vibration, and torque sensors. (14)(15)(16) Signal processing and analysis are important for enhancing the production capacity and processing quality. Most systems in nature are nonlinear systems, but humans linearize them to simplify their analysis. Linearization has less impact when the accuracy requirement is low, but it will produce considerable errors and misjudgments in precision production or in machines with extremely high accuracy requirements. Therefore, we need to clarify the nonlinear dynamic characteristics of cutting machine tools and provide relevant and useful parameters to build or design better sensors.
In the previous studies, some assumptions were made or linearization was performed to simplify the simulation model and reduce the simulation time. This may lead to large errors compared with the real situation. In this study, we consider the nonlinear dynamic responses in the cutting process of a machine tool with a nonlinear suspension effect and also take the nonlinear cutting force into consideration. The nonlinear dynamic equations are solved using the fourth-order Runge-Kutta and differential transform methods. The dynamic trajectories, power spectrum, Poincaré maps, and bifurcation diagrams are applied to analyze the dynamic motion. Figure 1 shows the model of metal cutting under a nonlinear suspension. K 1x and K 2x are the first and second equivalent stiffness coefficients in the vertical direction; K 1y and K 2y are the first and second equivalent stiffness coefficients in the horizontal direction; and C x and C y are the damping coefficients of the supported structure in the vertical and horizontal directions, respectively. Moreover, F x and F y are the components of the external excited cutting forces; F y is the cutting force dependence on the cutting speed and chip thickness; and F x is the thrust force, which is related to the main cutting force through a related frictional coefficient μ (F x = μF y ). The nonlinear parts of the dynamic equations include a nonlinear suspension term (hard spring case) and a nonlinear cutting and thrust force term. 3 2

Mathematical Modeling
where 1 2 is approximated as (1/2)[1 + tanh(h/ε)], and sgn(V f ) is approximated as tanh(v f /ε). The nonlinear dynamic equations presented in Eqs. (3) and (4) for the cutting system with nonlinear suspension effects and a strongly nonlinear cutting force were solved using the fourth-order Runge-Kutta method.
The differential transform method is an effective tool for solving differential equations owing to its rapid convergence and minimal computational error. Furthermore, compared with the integral transformation approach, the differential transform method has an additional advantage in that it can be used to solve nonlinear differential equations. This method has been widely used for solving nonlinear dynamic problems. Thus, we used this method to carry out the numerical analysis in this study (see Appendix A). The numerical data were then used to generate the dynamic trajectories, power spectrum, Poincaré maps, and bifurcation diagrams.
Equations (1) and (2) are of the second order and contain a cubic nonlinearity. Taking the differential transform of Eqs. (1) and (2) with respect to τ, they become These equations can be rewritten as

Stability Analysis of the Cutting System
The stability of the cutting system is analyzed using the characteristic equation of the system. If the chip thickness is evenly distributed, there will be periodic vibration and the cutting force F will remain unchanged. In addition, the deflection will remain fixed during the cutting process. When an impact occurs or the cutter passes through a notch, there will be a slight disturbance, and the output will become smaller and approach zero. This balance is called "stable" or "asymptotically stable". The stability is analyzed through the characteristic equation of the system. Let us take Eq. (1) into consideration, with F being the amplitude of the cutting force. Equation (1) is then replaced with We assume that Eq. (9) is under homogeneous initial conditions with (0) 0 x = and (0) 0 x =  . We perform the stability analysis by solving this differential equation with a fixed ω value and a small F value.

Numerical Results and Discussion
In the present study, the nonlinear dynamics of the cutting system shown in Fig. 1 are analyzed using the Poincaré maps, bifurcation diagrams, Lyapunov exponent, and fractal dimension. The dynamic trajectories of the cutting system provide a basic indication as to whether the behavior of the system is periodic or nonperiodic. However, they are unable to identify the onset of chaotic motion. Accordingly, some other analytical method is required.
The dynamics of the cutting system is analyzed using the Poincaré maps derived from the Poincaré section of the system. A Poincaré section is a hypersurface in the state space that is transverse to the flow of the system of interest. In non-autonomous systems, points on the Poincaré section represent the return points of a time series corresponding to a constant interval T, where T is the driving period of the excitation force. The projection of the Poincaré section on the y(nT) plane is referred to as the Poincaré map of the dynamic system. When the system performs quasi-periodic motion, the return points in the Poincaré map form a closed curve. For chaotic motion, the return points form a fractal structure comprising many irregularly distributed points. Finally, for nT-periodic motion, the return points have the form of n discrete points. The spectrum components of the motion performed by the cutting system are analyzed by using the fast Fourier transform to derive the power spectrum of the displacement of the dimensionless dynamic transmission error. In the analysis, the frequency axis of the power spectrum plot is normalized using the rotational speed ω.
In the present analysis, bifurcation diagrams are generated using two different control parameters: the dimensionless damping coefficient c and the dimensionless rotational speed ratio s. In each case, the bifurcation control parameter is varied with a constant step and the state variables at the end of one integration step are taken as the initial values for the next step. The corresponding variations of the y(nT) coordinates of the return points in the Poincaré map are then plotted to form the bifurcation diagram. The Lyapunov exponent of a dynamic system characterizes the rate of separation of infinitesimally close trajectories and provides a useful test for the presence of chaos. In a chaotic system, the points of nearby trajectories starting initially within a sphere of radius ε 0 form an approximately ellipsoidal distribution with semi-axes of length ε j (t) after time t. The Lyapunov exponents of a dynamic system are defined by The attractors in the embedding space were constructed using a total of 60000 data points taken from the time series corresponding to the displacement of the system. By trial and error, it was found that the optimum delay time when constructing the time series corresponds to onethird of a revolution of the system. The reconstructed attractors were placed in embedding spaces with dimensions of n = 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20, yielding 10 different log c(r) vs log r plots for each attractor. The number of data points chosen for embedding (i.e., 60000) reflects a compromise between the computation time and the accuracy of the results. The number of points used to estimate the intrinsic dimension of the attracting set in the current analysis is less than 42 M , where M is the greatest integer value less than the fractal dimension of the attracting set.
In this study, we compare the numerical results of applying the fourth-order Runge-Kutta and differential transform methods, which are found to be almost the same. The Runge-Kutta method is a well-developed numerical method and the differential transform method is demonstrated to be as effective as the Runge-Kutta method in this case. The time step in the iterative solution procedure was assigned a value of π/300 and the termination criterion was specified as an error tolerance of less than 0.0001.
In practical cutting systems, the dimensionless damping coefficient c is commonly used as a control parameter. Accordingly, the dynamic behavior of the current cutting system was examined using c as a bifurcation control parameter. Figure 2 presents bifurcation diagrams for the cutting system displacement plotted against c. The bifurcation diagrams show that the geometric center of the cutting system performs nonsynchronous motion in the horizontal and vertical directions. Strongly nonperiodic or even chaotic motion occurs at lower c and the nonperiodic motion becomes periodic at higher c. The vibration amplitude in the horizontal direction also decreased at higher c values (c > 0.0575). The above simulation result seems to be  consistent with real phenomena. The dynamic responses of the cutting system in the vertical direction are markedly different from those in the horizontal direction. As c increases, the dynamic motion still shows a nonperiodic response, even for c > 0.0575. Thus, we found very interesting nonsynchronous motion in the vertical and horizontal directions, especially at higher damping coefficients. Although the cutting system or another vibrating machine system is known to become steady at a higher damping coefficient, actually the suspension of this system is highly nonlinear. Thus, we have found synchronous behaviors in the vertical and horizontal directions of the cutting system, which may provide interesting or valuable information for analyzing or controlling this type of system. Figures 3-6 show the phase diagrams, power spectra, Poincaré map, Lyapunov exponent, and fractal dimension of the cutting system for chaotic motion at c = 0.015 (Figs. 3 and 4) and 0.025 (Figs. 5 and 6) in the vertical and horizontal directions. These simulation results show that the dynamic responses are synchronous in the vertical and horizontal directions. Also, the phase diagrams show disordered dynamic behaviors, the power spectra reveal numerous excitation frequencies, the return points in the Poincaré maps form geometric fractal structures, the maximum Lyapunov exponent is positive, and the fractal dimensions are found to be noninteger. Thus, we may conclude that the dynamic motion is chaotic for the above control parameters, with the simulation results obtained by the different methods corresponding to one another. The dimensionless rotation speed s is also an important control parameter for analyzing the dynamic responses of rotating machines.  Figure 7 shows bifurcation diagrams for the geometric center of the cutting system without a nonlinear suspension using the dimensionless damping ratio c as a bifurcation parameter. Compared with the cases with a nonlinear suspension, the simulation results show that the dynamic motion behaves differently in the horizontal direction, particularly at low c values, although similar behavior is observed in the vertical direction. This may be due to the increased absorbance of energy of the system and the suppression of vibration at higher damping ratios, even for highly nonlinear systems. Figure 8 presents bifurcation diagrams for the dimensionless displacement of the cutting system in the vertical and horizontal directions using the dimensionless rotating speed s as a bifurcation parameter. It can be observed that the cutting system exhibits periodic motion at low rotation speeds and nonperiodic or even chaotic motion at high rotation speeds. We also found that the motion is synchronous in the vertical and horizontal directions.
The results of a stability analysis of the cutting system in the vertical direction are shown in Fig. 9. We found that the system behaves stably in the red region and unstably in the blue region. A sudden and unexpected jump was also found at ω = 2.71. For ω > 2.71, the corresponding cutting force F on the boundary curve becomes large. This jump frequency may be the high resonance frequency of this system. We conclude that the cutting force should be greater with increasing frequency of the system.

Conclusions
We showed that chaotic behavior exists in a cutting system with a nonlinear suspension and nonlinear cutting force. Some interesting and useful simulation results were also found. In particular, we found that the dynamic responses behave nonsynchronously in the vertical and horizontal directions with increasing dimensionless damping coefficient. It is well known that if the behavior of a nonlinear dynamic system is chaotic, the resulting broad band vibration with a comparatively large vibrational amplitude will enhance the probability of fatigue failure. To increase the working life of a cutting system or enhance its performance, it is important not to operate the whole system under chaotic motion. Therefore, this study may help provide a theoretical understanding of nonlinear systems of cutting machine tools and avoid undesired dynamic responses in machining. Furthermore, the results of this research will help clarify the dynamic characteristics of cutting machine tools and provide relevant and useful parameters for building or designing better sensors.
This equation reveals that the kth term of the differential transformation of dx(t)/dt can be represented by the (k+1)th term of the differential transformation of the original function x(t).
The differential transformation of the product of two functions x(t) and y(t) is This equation reveals that the kth term of the differential transformation of the product of x(t) and y(t) can be represented by the first k terms of the differential transformations of the original functions x(t) and y(t). Using the differential transform method, a differential equation can be transformed into an algebraic equation in domain K. The differential transform X(k) of the unknown function x(t) can be achieved by solving the iterative equation. After that, a solution of the unknown function x(t) in the form of a power series can be obtained by the inverse differential transform of X(k), as shown in Eqs. (A4) and (A5). To speed up the convergence and improve the accuracy of the calculation, the entire domain D is usually split into sub-intervals. First, the differential transform method is applied to solve the original equation in the first sub-interval. After that, the final values of the first sub-interval are adopted as the initial values of the second interval and the original equation is solved under these new initial values. The same procedure is repeated in all the later sub-intervals until the solution of the entire domain D is obtained.