Open access

Finite Element Analysis of Machining Thin-Wall Parts: Error Prediction and Stability Analysis

Written By

YongAn Huang, Xiaoming Zhang and Youlun Xiong

Submitted: December 13th, 2011 Published: October 10th, 2012

DOI: 10.5772/50374

Chapter metrics overview

3,479 Chapter Downloads

View Full Metrics

1. Introduction

Thin-wall workpieces are usually low rigidity and complex shapes, which results in great challenges in machining. The machining accuracy, both geometric accuracy and surface integrity, plays a significant role in achieving overall product’s functional performance. Machining is a key process since it is easy to result in deflection and chatter[1, 2]. In stable machining, surface dimensional error due to the workpiece deflection affects machining precision[1, 3]. In unstable machining, the chatter becomes a critical problem for high surface quality. Currently, the existing methods for predicting the cutting forces, deflections and stability can be divided into experimenting and modeling. The direct trial-and-error approach is often expensive and time consuming. Producing the right profile in such parts increasingly depends on specialized CAD/CAE/CAM packages for defining appropriate cutting strategies and tool paths[4-7]. However, most of the existing techniques and models are based on idealised geometries and do not take into account factors such as variable cutting force, part/tool deflection, machining stability[8].

In many cases, the parameters of milling system are uncertain derived from the measurement errors, system nonlinear behavior, use or not use of coolants and other environmental noise. In addition, the stability boundaries are highly sensitive to the milling system parameter uncertainties. Therefore it is questionable that the usefulness of stability Lobes for estimating the milling stability obtained by applying deterministic parameters. Nominal machining parameters from the deterministic machining parameters optimization cannot guarantee the stability of milling process and lead to an actual maximization of material removal rate(MRR) and minimization of surface location error (SLE

Surface location error is defined as the error in the placement of the milling cutter teeth when the surface is generated.

). In practice, the cutting forces that statically and dynamically excite the tool and part, reducing the validity of the CAD/CAE/CAM output and leading to additional machining errors. For example, it induces chatter and large deformation in thin-wall workpiece, and very high cutting temperature which causes excessive thermal expansion specifically under dry machining conditions. Finite element analysis (FEA) has been widely adopted in numerical simulation of the machining process. FEA-based simulation, considering physical factors, such as material properties, tool geometry etc., is required to accurately predict the deflection and stability. Normally, the aim of cutting parameters optimization is to improve the part quality, maximize MRR and the constraint conditions are to keep the spindle speed below a predifined one and require the cutting process stable.

The cutting force and the part deflection are usually solved by an iterative simulation algorithm. Wang et al[9] studied force-induced errors, and developed a quasi-static error compensation method. Law et al.[10] calculated the cutting force from the measured milling torque, and integrated both the force and deflection models to develop a compensation methodology. However, currently the force modeling research has been mainly focused on theoretical rigid force models or mechanistic force models [11]. Budak and Altintas[12] considered the peripheral milling of a very flexible cantilever plate that incorporate a mechanistic force model and finite element methods. Feng and Menq [13, 14] developed a cutting force model taking into account the engaged cut geometry, the undeformed chip thickness distribution and the effect of the cutter axis offset. Budak et al.[15] proposed an oblique cutting mechanics model, in which the oblique cutting force are obtained from orthogonal cutting force. Ratchev et al.[16] proposed a flexible force model to study the tool deflection based on an extended perfect plastic layer model. They considered the end-milling cutter as a cantilever with the force acting at the cutter tip centre position[11]. However, the dynamic effect and the generated heat during machining were not considered, and the tool and the workpiece were assumed to deform to their static equilibrium position at any milling instant[1, 17]. Spence et al.[18] indicated that most existing machining simulation techniques were geometric and ignored the physical aspects of the process. Therefore, for making an appropriate choice of the cutting parameters and the operation sequence, deep understanding of the induced cutting deformation and the heat transfer is necessary. In the recent past, various techniques have been developed for studying the force- and temperature-induced aspects of machining, but separately. There is still a lack of a comprehensive milling simulation model which, taking into account the effects of the tool-path, and the cutting variables, simulates the thermomechanical aspects of machining[19]. Their applicability to model force in machining of thin-wall workpiece is limited due to the variability of material properties, cutting force, non-linear dependency on tool immersion angle and chip thickness[11]. Increased attention is being focused on the development of a computationally efficient milling process model, well capable to perform thermomechanical analysis of the metal cutting process.

Although there are many mechanisms of instability or chatter as mentioned in Wiercigroch and Budak[20], instability due to regeneration of surface waviness is by far the primary cause of instability. The stability analysis of the milling system can be performed only by applying approximated numerical methods since there is no closed solution to time-delay differential equation of milling dynamics. Alternatively it can be carried out by means of time-domain simulations[21], in the frequency domain[22] or by applying methods based on delay differential equation theory, such as the semi-discretization method[23] and the time finite element approach (TFEA)[24]. Budak and Tekeli[25] proposed a method to determine the optimal combination of depths of cut, so that chatter free material removal rate is maximized, with application on a pocketing example and significant reduction in the machining time is obtained. Altintas and Merdol[26] developed a generalized process simulation and optimization strategy for increasing MRR while avoiding machining errors and considering the chatter, and spindle's torque or power limits. However, the works of milling stability and cutting parameters optimization were addressed little to take into account the parameters uncertainties. Duncan et al.[27] used the Monte Carlo method to determine the associated uncertainties in the stability limit at each spindle speed. However the estimated stability intervals are too large to supply a useful advice to the parameters selection in an actual milling process. Kuidi[28] studied the robust optimization of SLE and MRR in milling process with uncertainties. The deterministic optimization formulation was modified to account for the axial depth uncertainty. But the uncertainties of Lobe diagram and SLE, which determine the milling stability and part quality were not taken into account. And the detailed optimization procedures are absent. Totis[29] used a new probabilistic algorithm for a robust analysis of stability in milling process, which performs the stability analysis on an uncertain model. The main objection to the general use of probabilistic analysis techniques is that non-deterministic properties cannot always be exactly represented using the probabilistic concept. Indeed, probabilistic methods cannot deliver reliable results unless sufficient experimental data are available to validate the assumptions.

In this chapter, we describe the force and temperature-induced error prediction and stabilty analysis of milling process by using the finite element method. Firstly a thermo-mechanical analysis is established to predict the force- and temperature-induced deflection. Then, TFEA is adopted as the deterministic model to obtain SLE and milling stability Lobe diagram. The uncertainties of modal shape parameters of spindle-tool system are investigated and sensitivity analysis is used to evaluate the upper and lower bounds of SLE and stability Lobe diagram. Finally, the formulation of robust spindle speed optimization is given to minimize the maximal SLE in a milling process and constraint condition is to maximize spindle speed and keep the machining parameters below the lower bound of Lobe diagram. With two optimization results, derived from robust and deterministic optimization formulations, experimental verifications are given.


2. Flexible thermo-mechanical model for error prediction

The surface error is induced by the deflection of tool and workpiece. Most of the existing techniques are based on idealized geometric profile and do not take into account tool/workpiece deflection, which results in a significant deviation between the planned and machined workpiece profiles[11]. This chapter focuses on deflection of thin-wall workpieces induced by cutting force and temperature. Reliable quantitative predictions of the cutting force and temperature components are essential for determining geometric errors of machined workpieces. Force predictions are required for arriving at constrained optimization strategies in computer-aided process planning[15].

2.1. Flexible force model for thin-wall workpiece

The research in the area of cutting covers a very wide range, since there are many independent influencing factors, including the cutting parameters, material properties, properties of the machine–tool–workpiece system, and tool geometry. A general model for the determination of cutting forces in ball-end milling operations was presented in reference[30]. The mathematical model represents the relations between the cutter and the workpiece, the change of chip thickness and the milling cutter rotation angle. The cutting forces are divided into several parts which depends on the number of cutting edges, cutting edge length and milling cutter rotation angle. Figure 1 shows the schematic diagram of a cutter of a ball-end milling and its configuration parameters.

Figure 1.

Constant lead spiral cutting edge for ball-end milling

The following is the main idea how to calculate it [31, 32]. The lag between the tip of the flute at z=0 and at axial location z is[30]


where, R0is the radius of the hemispherical part, i0is the helix angle of the cutting edge. A point on the flute j at height z is referenced by its angular position in the global coordinate system,


where Nf is the number of flutes, θis the tool rotation angle.

The angle position in the direction of z-axis from the center of the hemispherical part to the point on the cutting edge can be written as


The tangential, radial and binormal components are calculated as[15, 32, 33],


where, tn(Ψ,θ,κ)=fzbsinΨsinκis the uncut chip thickness normal to the cutting edge, and varies with the position of the cutting point, Ktc, Krc, Kac(N/mm2) are the shear specific coefficients, Kte, Kre, Kae(N/mm) are the edge specific coefficients. dS=(R'(φ))2+R2(φ)+R02cot2i0dφ(mm) is the length of each discrete elements of the cutting edge, fzbis the feeding per tooth, κ=arcsin(R(Ψ)/Rb), and db (mm) is the differential length of cutting edge. In many mechanistic models for the milling process, the milling force coefficientsKtc, Krc, Kac, Kte, Kre, Kaeare established from specially devised milling tests (e.g. average cutting force coefficient model) followed with linear regression analysis. Usually there are two methods to predict the parameters which are mechanistic evaluation and prediction from an oblique cutting model[15].

The cutting forces in Eq.(4) are modeled in terms of two fundamental phenomena, an edge force component due to rubbing or ploughing at the cutting edge, and a cutting component due to shearing at the shear zone and friction at the rake face[15]. The cutting force model including explicitly the ploughing component can obtain more precise prediction accuracy. Once the tangentialFt(θ,z), radialFr(θ,z), and axial Fa(θ,z) cutting force at the tooth-workpiece contact point are determined, the resultant forces in Cartesian coordinates are obtained by introducing the transformation matrix[32],


where, {dFx,y,z}=[dFxdFydFz], [T]=[sinκsinΨcosΨcosκsinΨsinκcosΨsinΨcosκcosΨcosκ0sinκ],{dFr,t,a}=[dFrdFtdFa].

Then one can get [32]


Since the cutting force coefficients (Ktc, Krc,Kac) may be dependent on the local chip thickness, the integrations given above should be calculated digitally by evaluating the contribution of each discrete cutting edge element at dz intervals. There are many models to calculate the coefficients, such as bi-linear force model, exponential chip thickness model, high-order force model, and semi-mechanistic model. Budak et al[15] presented the cutting transformation model.


where, τis shear stress at the shear plane, ϕnand βn are normal shear and friction angles in oblique cutting respectively, αnis normal rake angles, ηcis chip flow angle in the rake face, βis friction angle at the rake face. For example, when Cutter radius is 3mm, Feed is 0.02mm/rev, Depth of cutting is 0.02m, the cutting force can be calculated as Figure 2.

Figure 2.

Predicted cutting forces for slot cutting tests

Most of reported papers in the area of cutting force- induced error belongs to those caused by large deformation of thin-wall workpieces under load[34, 35]. Peripheral milling of flexible components is rather complicated due to periodically varying cutting-forces. These forces statically and dynamically excite the tool and workpieces, which leads to significant and often unpredictable deflections. Static deflections produce dimensional profile errors, and dynamic displacements lead to poor surface finish quality and machining stability problems in addition to dimensional profile errors[16, 36]. In the known condition of cutting force, the deflection of the thin-wall workpiece can be calculated through finite element analysis in real time.

There are two kinds of force model, theoretical rigid force model and adaptive theoretical flexible force model. High complexity is associated with modeling of cutting forces in machining due to the variable tool/workpiece deflection and changing tool immersion angle. To address this complex dependency an interactive approach integrating an extended perfect plastic layer force model is adopted to link force prediction with workpiece deflection modeling[8]. The predicted profile of the workpiece is utilized to identify the “real” material volume that is removed during machining, instead of the “ideal” one defined by the “static” NC simulation packages[11]. In milling a thin-wall workpiece, the differential cutting force on the engaged infinitesimal tool cutting edge varies with the cutting depth that is effected by workpiece deflection[17]. The force is calculated by taking into account the changes of the immersion angle, φ, of the engaged teeth. As soon as the deflection, uyand the coordinate (x, y, z) of point, a, are known (Fig. 4), the instant immersion angle ϕ in milling after deflection can be calculated using[1]


where R is the cutter radius, and hr is the designed milling depth in the workpiece thickness direction, while uy is the deflection in the corresponding point predicted through FEA. Ratchev et al[3, 11] proposed a flexible force model for machining dimensional form error prediction of thin-wall components. Here, the thermal deflection is introduced into the the flexible force model to consider the temperature effects.

FEA softwares are used to calculate the deflection caused by the cutting force at each sampling point through the following equation,


where [K] is the stiffness matrix of the workpiece, {U}=[Ux,Uy,Uz]Tand {F}=[Fx,Fy,Fz]T are nodal displacements of workpiece and the external cutting force acting on the tool-workpiece transient contact surface, respectively. As the boundary conditions are specified, the nodal displacements can be obtained through solving the above equations.

2.2. Thermo-mechanical analysis for thin-wall workpiece

The errors are usually caused by excessive deformation at the interface due to the cutting force and temperature, so they have to be considered simultaneously. With the recent developments in machining automation, various cutting temperature measurement techniques, including tool-work thermocouple, embedded thermocouple, and infrared pyrometer, emerged[37-41]. The aforesaid experimental techniques have been widely applied in machining due to its simplicity. Lin[42] and Kwon[43] studied the transient interfacial temperature and heat dissipation in the workpiece during a slot milling process. Fang and Zeng[44] utilized FEM to develop a 3D model of the oblique cutting process for the prediction of temperature distribution in the workpiece, tool and chip. Temperature distribution in the workpiece was estimated for a simple single pass slot milling operations in only a few reported works[45, 46]. However, these studies ignore the structural analysis for predicting part deformations under coupled thermomechanical loading conditions.

The milling techniques still face to a severe problems of inducing very high cutting temperatures causing excessive thermal expansion of the workpiece, especially in dry machining. The cutting parameters, namely, cutting speed and feed rate, have the greatest influence on the cutting temperature. For analyzing the phenomenon of heat dissipation into the workpiece and its influence on part deformation, a 3D model of the transient milling process was developed based on commercial FEM program, such as Abaqus, Ansys and Comsol. These systems allow: (1) creating 3D FEM models of the fixture–workpiece configurations, (2) applying appropriate materials for the workpiece, (3) applying appropriate machining boundary conditions, and (4) performing transient thermal and structural analysis where the transient temperature distribution profiles are applied along with the cutting forces to predict part deflections.

If the thermal conductivity (K) and heat capacity (ρCp) of the work material is higher, then the generated heat is more readily conducted into workpieces and causes thermal expansions which produce severe irregularities. The cutting force more easily induces large deformation in the workpiece. In return, the deformation resulted from the cutting force and the generated heat will change the cutting parameters. However, most models are based on idealized geometries and do not take into account the factors, such as variable cutting force, thermal load, part/tool deflection, etc[1]. Before calculating the temperature distribution, cutting temperature at the interface need to be predicted. For speed-, and feed-dependent boundary conditions in machining, the thermal source problem in machining thin-wall workpiece is very difficult to be solved analytically. Therefore, the temperature prediction at the interface is usually achieved by non-linear empirical modeling approaches. The average interface temperature, as measured by the tool-work thermocouple, is[47]


where V is cutting speed (m/s), dis depth of cut (mm), and f is feed (mm/rev). The above empirical model of the interface temperature is developed for turning of 4140 Steel alloy with tungsten carbide tools. The relation of cutting temperature and feed is plotted in Fig.3 where the temperature at the interface is very high.

Figure 3.

Cutting temperature versus cutting feed, depth of cut 0.763mm, cutting speed 3m/s

The temperature variations are related to the heat source movement, heat source intensity, and thermal resistance coefficient. The fundamental generalized problem to be solved analytically is the heat conduction in a thin infinite plate with a convective and radial boundary condition on the face. The time-dependent heat transfer process is governed by the following differential equation[48],


where, η=h/kw, his the convection coefficient of heat transfer, kis the thermal conductivity of the material, wis the plate thickness, αis the thermal diffusivity of the material, and g(r,t) is the internal heat generation rate per unit volume, the variable r is the radial distance from the heat source.

Figure 4.

Temperature-effected machining: (a) Top view;(b) Side view

The dynamic non-uniform temperature distribution roots in time/position-dependent thermal deformations. Eq.(13) is to be solved by finite element method here. Ref. [48] proposed another computational methods. There are two basic thermal error modes, namely, thermal expansion and thermal bending. Figure 4 is the schematic diagram of machining results due to temperature effects. Figure 4a-4b show the deflections from the top and side views respectively. It can be observed from Figure 4 that the machining topography will be very complex when the temperature effects are considered. In return, the thermal expansion and thermal bending affect the cutting force. It is a complex, back and forth cycle, which has to be simulated by interactive algorithm.

2.3. Interactive algorithm for thermo-mechanical analysis

The cutting forces depend on the chip thickness which is a function of the tool immersion angle. The tool immersion angle is a function of the part deflection which itself depends on the cutting forces[47]. The deflection is determined by the cutting force, as well as temperature distribution. The cutting forces depend on the chip thickness which is a function of the tool immersion angle, and the machining temperature is a function of the cutting speed, depth of cut and the feed. An iterative procedure is used to determine the milling error. The predicted deflection from CAE is used to identify the practical material volume that is removed during machining instead of the ideal one defined by CAM. The thermo-mechanical analysis is outlined in Figure 5. While milling a thin-wall workpiece, as soon as the cutter is engaged, the workpiece deflects to a new position, and the cutting temperature also changes at the same time. In more detail, the tool-workpiece intersection line can be used as an example to explain the impact on the cutting force and deflection[16].

Figure 5.

Flowchart of complex surface machining based on thermo-mechanical analysis

Although commercial FEA softwares are also used to simulate manufacturing processes, they cannot be solely used to simulate multi-step cutting processes of thin-wall parts. The main difficulty is that material removal and remeshing of part model are very complex for the multi-step processes, and all these FEA softwares do not integrate an appropriate theoretical force model for workpiece/tool deflection prediction, so varieties of models and software are involved[1]. There is a need to link the mainstream commercial FEA software with force prediction models, thermal prediction models and material removal models in order that data exchange among them can be achieved.

2.4. Error definition and thermo-mechanical deflection

The surface dimensional error is the normal deviation of the actual machined surface from the desired machined surface. For example, at point P in Figure 6, the corresponding surface dimensional error iseP=δt,P+δf,P, where δt,P and δf,P are the normal projections of the temperature- and force-induced deflection corresponding to PointP, respectively. For the convenience of study, the distance between the initial surface to be machined and the desired machined surface is named as the nominal radial depth of cut symbolized byRN. In actual machining, to ensure that the surface dimensional error does not violate the tolerance, RAis often specified to be different from RN(RARN)[3]. In this case, eP=δt,P+δf,Phas to be adjusted to the calculation of surface dimensional error


Note that RN and RA are the nominal and specified radial depth of cut, respectively. For a certain surface generation line, the steps adopted to calculate the error distributions can be found in Ref.[3].

Figure 6.

Definition of the surface dimensional error

2.5. Examples of thermomechanical models

For simplification the part is assumed to be a thin-wall rectangular workpiece. The required machined profile is a flat surface parallel to the planeOXY. During milling, the workpiece deflection in its thickness direction has a significant impact on forming the surface profile error. The contributions of the workpiece deflection in the feed direction and the tool axial direction can be ignored. Therefore, the investigations focus only on error prediction in y-axis direction. The simulations were based on clamped-free-free-free cantilever plates with dimension 150×120×5mm3and Aluminium alloy 6063 T83. The quasi-static cutting force is treated as a moving-distributed load acting on the workpiece-tool contact zone in the milling process. To compute the workpiece/tool response to the cutting force, the continuous machining process was simulated by multi-step cutting processes[18]. Comsol is employed to estimate the force- and temperature-induced deflection.

The input to finite element model is the chip-related cutting force and temperature, and a set of parameters describing material properties, boundary conditions and other constraints. In order to simplify the complex simulation, we employ an assumption that the instantaneous stiff variation due to material removal will be not taken into account. To compute the workpiece response, the continuous machining process was simulated by multi-step cutting processes. During machining, the tool moves along the machining surface. Here, we took a different approach that uses a moving coordinate system fixed at the tool axis. After making the coordinate transformation, the heat transfer problem becomes a stationary convection-conduction problem in COMSOL. With the workpiece undergoing deformation, the cutting force is also changing considerably. These changes is taken into account by computing the cutter on a moving mesh attached to the workpiece.

2.5.1. Machining based on a designed tool path

The results based on the rigid force model are discussed. Supposed that the force is0.23N/mm2, the temperature is 1200℃, and the two loads act on the interface workpiece and cutter at neighborhood of x=75mm of the workpiece. It can be noted that the maximal deflections are 0.348mm and 0.384mm in Figure 7-10 in the positive direction and 1.5mm in the passive direction, respectively. Figure 9 is the force-induced and temperature-induced deflection of the top edge of the workpiece. It can be observed that the deflections at x=75mm in the two figures have different directions. The force-induced deflection is along the positive direction, however the temperature-induced deflection along the passive direction. In reported studies, the force-induced and the temperature-induced deflection are investigated dividually, which will leads large error in machining.

Figure 7.

Deflection results of the rigid force model

Figure 8.

Thermal deflection results

Figure 9.

Deflection of the top edge of the workpiece due to force and temperature

2.5.2. Machining based on flexible force model

Figure 10 is the results based on the flexible force model. The maximal deflection is 0.278mm. Comparing with Figure 7, one can observe that the max stress and the max deflection in Figure 10 are all smaller. The simulation shows the deflection of workpiece affect the cutting force, and the cutting force also affect the deflection in return. We can get the following conclusion: (1) if the workpiece is machined based on CAM, there exists large errors, and (2) the prediction from the flexible force model is smaller than that from the rigid force model. Obviously, the flexible force model is closer to the real machining that the rigid force model. The temperature is not considered in both models.

Figure 10.

Stress and deflection results of the flexible force model

2.5.3. Machining based on thermo-mechanical analysis

The results plotted in Figure 11 and Figure 12 are based on the thermo-mechanical analysis. The maximal deflection in Figure 11 is 0.527mm. It is larger than the deflection in the rigid force model and the flexible force model, but smaller than the sum of 0.278mm of the flexible force model and 0.384mm of the temperature-induced deflection. So the thermo-mechanical model is not a simple combination of the flexible force model and the thermal deflection. The temperature has heavy effect on the deflection that the curve in Figure 12 is very similar with the curve in Figure 9. In practical error compensation, there will overcut or undercut if the force and the temperature studied separately. In order to reduce the machining error, it is important to reduce the cutting temperature during machining.

Figure 11.

Deflection results of the thermo-mechanical analysis

Figure 12.

Deflection of the top edge of the workpiece due to thermo-mechanical action


3. Stability analysis of machining of thin-wall parts [49]

3.1. Milling dynamics

The standard two degree of freedom milling process is shown in Fig.1. The tool is assumed to be compliant relative to the rigid workpiece. The vibration is excited by the summation of cutting force. The governing equation of motion has the following form


where, M=[mx00my]C=[cx00cy]K=[kx00ky]and F are the modal mass, damping, stiffness, and cutting force matrices, respectively. The termsmx,y, cx,y, kx,yand Fx,y are the corresponding components in the flexible directions of the system. bis the axial depth of cut. τ=60/NΩis the tooth passing period in seconds, in which N is the number of teeth on the cutting tool and Ω the spindle speed in rpm.x(t)=[x(t)y(t)]T is the dynamic response vector and x(t)x(tτ) the dynamic chip thickness, as shown in Figure 13.

Figure 13.

Dynamic chip thickness

Kcis a matrix given as follows, which represents the component of cutting forces that depend on the position vector


where Kt and Kn are the tangential and normal cutting force coefficients components, respectively.s=sinφj(t), c=cosφj(t)and the function g[φj(t)] acts as a switching function, which is equal to 1 if the jth tooth is active and 0 if it is not cutting.

g[φj(t)]={ 1     φe<φj(t)<φa 0     otherwiseE17

where φe and φa are the angles where the jth tooth enters and exits the cut, respectively. For down-milling operation, φa=π, for up-milling,φe=0. Note that the entry and exit angles may vary due to heavy vibrations of the tool. This effect is neglected here, and the angles φe and φa are approximated by constant values as it is usually done in the literature. fo(t)is the stationary cutting force vector (fzis the feed per tooth):


3.2. Deterministic model for predicting milling stability and calculating SLE

Here the TFEA is introduced for brevity as a basis for the further uncertain analysis. The initial work of applying TFEA to the delay equations can be found in Ref. [5]. The main idea of TFEA is that the dynamic behavior of the milling process is governed by TFEA as a discrete linear map which relates the vibration response while the tool tooth is engaged in the cut, which depends on previous tooth passages and therefore includes the time delayτ, to free vibration while the tooth is not engaged in the cut. The dynamic map is expressed as


where Ais the state transition matrix, the size of which depends on the number of time finite elements and polynomial order representing one time period. Bis a vector that depends on the process parameters. qand q˙ are the sets of x and ypositions and velocities for all nodal times in one tooth passage, respectively. Stability of the milling process is determined from the eigenvalues ofA, i.e.λ(A). The maximum magnitude of the map eigenvalues is described by


where λk denotes the kth eigenvalue of the dynamic map matrixA, which is a function of the cutting conditions. Unstable conditions exist if


The stability boundary is defined by the boundary curve of axial depth b and spindle speedΩ. A combination of band Ω values below the stability boundary, blim, gives stable cutting conditions, whereas a combination above the stability boundary leads to an unstable cut. The stability boundary corresponds to the cutting conditions at which


When the milling process is stable, SLE can be obtained from fixed points of the dynamic map as given in Eqs. (21):


Substituting Eqs.(23) into Eqs.(19) gives the solution of fixed point map:


The solution of fixed point displacement can be obtained and used to specify SLE as a function of cutting parameters.

Figure 14.

Nyquist plot of FRF of spindle-tool system with uncertainty regions[51]

3.3. Upper and lower bounds of Lobe diagram and SLE

Since the uncertainties occur in the experimental modal analysis, the measured values for the structural parameters given in Eq.(15) are different from the actual ones. Much experiments had been done so that the errors can be observed, see Ref.[50, 51]. A detailed study was carried out in the work[51] to discuss the modal test. As shown in Figure 14, Nyquist plot of FRF of the spindle-tool system with uncertainty regions was given.

In this section, we address the uncertainty problems in model parameters, which are modeled as the interval parameters as follows:


where parametersM, Cand Kare assumed to be bounded in the intervals and can be stated as:


The problem can be formulated as given the uncertainties as shown in Eqs.(26) how to estimate the stability of the milling system and SLE (Eqs.(25)). Noted the linear transform matrix A given in Eq.(19), the right eigenvalue problem can be stated as


where λ is the eigenvalue and u the corresponding eigenvector. The accompanying left eignenvalue problem is


Combined Eqs.(27) and (28), we can obtain


For a parameterz, i.e., a spindle-tool modal parameter or a part geometrical one, the sensitivity of system eigenvalue with respect toz, λ/zcan be expressed as follows


whereA/z, the derivative of the matrix A is obtained by the central difference method. So the sensitivity expression in Eq.(30) can be treated as a semi-analytical one[52]. It should be noted that


where λmax'is the complex conjugate ofλmax, then


Substituting Eq.(30) into Eq.(32), we can get the explicit expression for|λmax|/z. More details of sensitivity of stability boundary, especially the error analysis can be referred with the reference[52]. When there is a perturbation for the milling system parameterz, i.e., zz+δz, the increment of λmax is


Denote λmaxU and λmaxL the upper and lower bounds of the eigenvalues of milling systems, respectively. They read

λmaxU=max{|λmax+δλ|,  |λmaxδλ|}λminL=min{|λmax+δλ|,  |λmaxδλ|}E34

Then the upper and lower bounds of the milling stability Lobes, cUand cL can be obtained as follows

{cU  when  λmaxL(blim,Ω)=1cL   when  λmaxU(blim,Ω)=1E35

The central difference method used for the sensitivity analysis of system eigenvalues can be also applied to the sensitivity analysis of SLE. Thus we get the upper and lower bound of SLE, i.e. fSLEUandfSLEL. The detailed discussions are omitted here.

3.4. Robust machining parameters optimization

In general, a deterministic cutting parameters optimization problem can be stated as follows:

min {fSLE(b,Ω), Ω}s.t{λmax(A(b,Ω))1ΩΩ0 E36

where the optimization objects are to minimize SLE and maximize the spindle speed, and constraint conditions are to keep the milling process stable and the spindle speed is less than the predefined one, which is set by the spindle system and tool suppliers. Since there are uncertainties between the real structure and the model, the predicted performance including the milling stability and SLE will not be guaranteed. When the uncertainties exist, as given in Section 5, the optimized machining parameters, i.e. axis depth band spindle speed Ω obtained from the deterministic optimization formulation (36) cannot guarantee the stability of the milling process, and minimization of SLE and maximization of spindle speed cannot be achieved. As a comparison, we give a formulation of robust machining parameters optimization as follows:

min {max fSLE(b,Ω),Ω}s.t{λmaxU(A(b,Ω))1ΩΩ0 E37

where the optimization object is to minimize the maximal SLE and maximize the spindle speed, and constraint conditions are to keep the machining parameters band Ω below the lower bound of Lobe diagram given by λmaxU(A(b,Ω))1 and the spindle speed is less than the predefined one. The formulation of robust machining parameters optimization ensures the machining process stable and can lead to minimization of SLE under uncertainties.

3.5. Implementation of robust optimization formulation

The robust optimization formulation presented in Eqs.(37) can be transformed into the sequential optimization ones, that is:

min max fSLE(b,Ω)s.t{λmaxU(A(b,Ω))1ΩΩi  i=1,2,k for series of ΩiΩ0E38

The optimization problem given in Eq.(38) is solved by the augmented Lagrangian function method. The basic idea of this method is that it transforms the nonlinear optimization problem into an unconstrained optimization one by introducing a penalty function, named augmented Lagrangian function[53]. The augmented Lagrangian function LA(x,λ;μ) achieves these goals by including an explicit estimate of the Lagrange multipliers λ in the objective. The augmented Lagrangian function can be defined as


wherex=[b,Ω]T,c1=l(x),c2=ΩΩi, l(x)is a cubic polynomial curve to interpolate the obtained Lobe diagram. It can be easily to obtain the differentiation of the augmented Lagrangian function with respect to the machining parameters,xLA(x,μ;λ). Then we can obtain the following algorithmic framework.

Algorithm: the augmented Lagrangian function method.

Step 0: Given initial points x0andλ0, μ0>0, toleranceε0;

Step 1: Find an approximate minimizer xk ofLA(x,μk;λk), i.e.xk=arg minLA(x,μk;λk);

Step 2: IfxLA(x,μk;λk)>ε, set k=k+1and update Lagrange multipliers λk+1=λkci(xk)/μk i=1,2 to obtain λk+1 and go to Step 1; else exit and report the minimizerx*.

3.6. Experimental verification

In this section, experiments of cutting parameters optimization are carried out to verify the numerical optimization modeling results. First the modal shape parameters of the spindle-tool system are identified, with up and lower bounds, and then the optimization modeling results are presented by using the system parameters as input parameters. As a comparison, the deterministic optimization results are also presented. With the optimized machining parameters obtained from the robust optimization formulation and deterministic one, cutting experiments are implemented to check if the optimized cutting parameters are really “optimized”, which means there is no chatter arising with the optimized spindle speed.

3.6.1. System parameters identification

The system parameters as input variables to the optimization formulations are the structural dynamics parameters, i.e. modal shape parameters of the spindle-tool system. An impact hammer test is implemented to obtain the structural dynamics parameters of spindle-tool system. To obtain a reasonable interval, normally the repeated modal hamper impact experiments (see Fig.3) are necessary. Based on these experimental results, the modal parameters are identified by a Rational Fraction Polynomial (RFP) method[54], as given in Figure 16 and Figure 17. The mean values with upper and lower bounds for the modal parameters are given as follows (Table 1), which are derived from the repeated modal hamper impact experiments. MxandMy (kg) are the modal mass in x and y directions, respectively. ξxandξy (kg/s) are the modal damping coefficients in x and y directions, respectively. KxandKy (MN/m) are the modal stiffness, in x and y directions, respectively.

Figure 15.

Experiment setup for modal shapes of spindle-tool system

Figure 16.

Fitting of mode parameters of tool system in X-axis

Figure 17.

Fitting of mode parameters of tool system in Y-axis


Table 1.

Modal parameters of spindle-tool system

The tool for the simulation and experiments is a 10mm cylindrical cutter with 4-tooth and the cutter parameters are shown in Table 2.

Number of toothHelix angle
Edge length
Total length
104308.533Hard alloy

Table 2.

Cutter parameters

3.6.2. Modeling results

In the experiment of machining of an impeller blade(single blade), the upper limit of the spindle speed is 20000rpm. The cutting depth is 0.560mm, which is derived from the geometry of the semi-finished and finished workpieces and the given tool path[55].

In table 1, we have presented the upper and lower bounds of the modal shapes parameters of spindle-tool system with mean values. The deterministic Lobe diagram is calculated by the TFEA[24]. Using the proposed method given in Section 4, we can obtain the upper and lower bounds of the boundary curve of Lobe diagram and SLE, as given in the upper part of fig.8. After the lower bound of the Lobe diagram and upper bound of SLE have been obtained, they are selected as the constraint conditions in the robust optimization formulations given as:

min max fSLE(b,Ω)s.t{λmaxU(A(b,Ω))1ΩΩi  i=1,2,k for series of Ωi20000E40

As a comparison, the mean values of the modal shapes parameters and SLE are adopted in the deterministic ones, stated as:

min fSLE(b,Ω)s.t{λmax(A(b,Ω))1ΩΩi  i=1,2,k for series of Ωi20000E41

We solve the sequential optimization formulations as shown in Eqs. (40)-(41) by setting the spindle speedΩi=20000100*(i1),  i=1,2,51. That means we are interested in the spindle speed interval from 15000rpm to 20000rpm, and in this interval the optimization formulations are solved for 51 times. Augmented Lagrangian function method is adopted to solve the above two constrained optimization problems and the detailed procedures are given in Section 6. The modeling results are shown in Table 3.

SLESpindle speed
Deterministic optimization[1.88E-5m, 2.08E-5m]17900 rpm
Robust optimization1.99E-5 m19100 rpm

Table 3.

Comparison of the optimization results

3.6.3. Experiment of machining stability

The cutting experiments of impeller blade milling are implemented in Mikron HSM 600U, a five-axis NC machining center. Experiment setup including the machining setup and Labview signal acquisition interface are shown in Figure 18. First, we prepare two semi-finished blades, using the same workpiece geometry, tool path as well as machining parameters to ensure that the two blades are almost the same. Then the two semi-finished blades are used for the following finish milling with different spindle speeds.

The optimized spindle speeds derived from the robust and deterministic optimization formulations are used respectively to check the chatter occurrence. When the first blade is machined with spindle speed 19100 rpm, we observe that the milling process is stable and the quality of resulted workpiece is quite good, as shown in Figure 19(R). Then the second blade is machined and the spindle speed is set at 17900rpm, and the noise level dramatically becomes very high, which indicates the energy of the vibration at the frequencies related to chatter increases, see lower part of Figure 20. And we can see that chatter deteriorates the quality of the workpiece, as shown in Figure 19(L).

Spindle speed (rpm)Resonance (Hz)Resonance source
179001480natural frequency of structures
191001271spindle revolution frequency

Table 4.

Resonances of sound pressure signal from different spindle speeds

Further frequency response function (FRF) analysis of sound pressure signal indicates that when the spindle speed is at 17900 rpm, the resonance is about 1480Hz, which is approaching the first modal shape frequency of the spindle-tool system given byω01ξ2/2π. This indicates the occurrence of chatter. And when the speed is 19100rpm, the resonance is about 1271Hz, which is the spindle revolution frequency and its higher harmonicsknΩ/60, with k+ and n is the number of the tooth. We list the results in Table 4. From Table 4, we can see that the FRF analyses also suggest the chatter occurrence when the spindle speed is at 17900rpm and no chatter at 19100rpm.

Figure 18.

Experiment setup, (L)Machining setup (R)Labview signal acquisition interface

Figure 19.

Comparison of finished blades using different spindle speeds, (L) from deterministic optimization (R) from robust optimization

Figure 20.

Signal of sound pressure with different spindle speeds. (U) Upper and lower bounds of Lobe diagram (L) Signal of sound pressure


4. Conclusions

In machining thin-wall workpieces, the temperature and force-induced deflection contribute significantly to the surface error. The proposed methodology is based on coupling effects between cutting forces and temperature and their induced deflection during machining. There is still a knowledge gap in identifying the impact of deflection on the process of metal removal, and hence there are not systematic approaches to modeling, prediction of the component errors due to thermo-mechanical deflection in thin-wall structures. And we develop a robust spindle speed optimization formulation. The quantitative analysis on how the uncertainties in milling process affect the milling stability and SLE are presented. Comparing with the traditional deterministic spindle speed optimization formulation, our model can take into account the uncertainties, i.e. modal shape parameters of spindle-tool system and the resulted optimized spindle speed ensures the milling stability.


  1. 1. RatchevS.LiuS. L.HuangW.BeckerA. A.2007Machining simulation and system integration combining FE analysis and cutting mechanics modelling. Int J Adv Manuf Tech. 355565
  2. 2. HerranzS.CampaF. J.LacalleL. N. L. d.RiveroA.LamikizA.UkarE.SánchezJ. A.BravoU.2005The milling of airframe components with low rigidity: A general approach to avoid static and dynamic problems. Proceedings of the Institution of Mechanical Engineers. Part B. Journal of engineering manufacture. 219789801
  3. 3. WanM.ZhangW. H.QinG. H.WangZ. P.2008Strategies for error prediction and error control in peripheral milling of thin-walled workpiece. Int J Mach Tool Manu. 4813661374
  4. 4. RatchevS.NikovS.MoualekI.2004Material removal simulation of peripheral milling of thin wall low-rigidity structures using FEA. Advances in Engineering Software. 35481491
  5. 5. RatchevS.LiuS.HuangW.BeckerA. A.2004Milling error prediction and compensation in machining of low-rigidity parts. International Journal of Machine Tools and Manufacture. 4416291641
  6. 6. ZhangX. M.ZhuL. M.ZhengG.DingH.2010Tool path optimisation for flank milling ruled surface based on the distance function. International Journal of Production Research. 4842334251
  7. 7. DingH.BiQ. Z.ZhuL. M.XiongY. L.2010Tool path generation and simulation of dynamic cutting process for five-axis NC machining. Chinese Science Bulletin. 5534083418
  8. 8. RatchevS.HuangW.LiuS.BeckerA. A.2004Modelling and simulation environment for machining of low-rigidity components. J Mater Process Tech. 153-154: 67-73.
  9. 9. WangS. M.LiuY. L.KangY. A.2002An efficient error compensation system for CNC multi-axis machines. Int J Mach Tool Manu. 4212351245
  10. 10. LawK. M. Y.GeddamA.2003Error compensation in the end milling of pockets: a methodology. J Mater Process Tech. 1392127
  11. 11. RatchevS.LiuS.BeckerA. A.2005Error compensation strategy in milling flexible thin-wall parts. J Mater Process Tech. 162673681
  12. 12. BudakE.AltintasY.1995Modeling and Avoidance of Static Form Errors in Peripheral Milling of Plates. Int J Mach Tool Manu. 35459476
  13. 13. FengH. Y.MenqC. H.1994The Prediction of Cutting Forces in the Ball-End Milling Process.1. Model Formulation and Model-Building Procedure. Int J Mach Tool Manu. 34697710
  14. 14. FengH. Y.MenqC. H.1994The Prediction of Cutting Forces in the Ball-End Milling Process.2. Cut Geometry Analysis and Model Verification. Int J Mach Tool Manu. 34711719
  15. 15. BudakE.AltintasY.ArmaregoE. J. A.1996Prediction of milling force coefficients from orthogonal cutting data. J Manuf Sci E-T Asme. 118216224
  16. 16. RatchevS.LiuS.HuangW.BeckerA. A.2004A flexible force model for end milling of low-rigidity parts. J Mater Process Tech. 153134138
  17. 17. HuangY.LiuH.YinZ.XiongY. L.2009Complex Surface Machining: Thermo-mechanical Analysis for Error Prediction of Low-Rigidity Workpiece. Intelligent Robotics and Applications: 666677
  18. 18. SpenceA. D.AbrariF.ElbestawiM. A.2000Integrated solid modeller based solutions for machining. Comput Aided Design. 32553568
  19. 19. RaiJ. K.XirouchakisP.2009FEM-based prediction of workpiece transient temperature distribution and deformations during milling. Int J Adv Manuf Tech. 42429449
  20. 20. WiercigrochM.BudakE.2001Sources of nonlinearities, chatter generation and suppression in metal cutting. Philos T Roy Soc A. 359663693
  21. 21. CampomanesM. L.AltintasY.2003An improved time domain simulation for dynamic milling at small radial immersions. J Manuf Sci E-T Asme. 125416422
  22. 22. AltintasY.BudakE.1995Analytical prediction of stability lobes in milling. Cirp Ann-Manuf Techn. 44357362
  23. 23. InspergerT.StepanG.2004Updated semi-discretization method for periodic delay-differential equations with discrete delay. Int J Numer Meth Eng. 61117141
  24. 24. MannB. P.PhDWashingtonUniversity, 2003
  25. 25. BudakE.TekeliA.2005Maximizing chatter free material removal rate in milling through optimal selection of axial and radial depth of cut pairs. Cirp Ann-Manuf Techn. 54353356
  26. 26. AltintasY.MerdolS. D.2007Virtual high performance milling. Cirp Ann-Manuf Techn. 568184
  27. 27. DuncanG. S.KurdiM. H.SchmitzT. L.SnyderJ.2006Uncertainty propagation for selected analytical milling stability limit analyses. Transactions of NAMRI/SME. 341724
  28. 28. KurdiM. H.PhD.Universityof.Florida, 2005
  29. 29. TotisG.2009RCPM-A new method for robust chatter prediction in milling. Int J Mach Tool Manu. 49273284
  30. 30. MilfelnerA.KopacJ.CusF.ZuperlU.2005Genetic equation for the cutting force in ball-end milling. J Mater Process Tech. 16415541560
  31. 31. MilfelnerM.CusF.2003Simulation of cutting forces in ball-end milling. Robot Cim-Int Manuf. 1999106
  32. 32. LeeP.AltintasY.1996Prediction of ball-end milling forces from orthogonal cutting data. Int J Mach Tool Manu. 3610591072
  33. 33. LacalleL. N. L.SanchezJ. A.SalgadoM. A.2004Cutting force estimation in sculptured surface milling. Int J Mach Tool Manu. 4415111526
  34. 34. RameshR.MannanM. A.PooA. N.2000Error compensation in machine tools- a review Part I: geometric, cutting-force induced and fixture-dependent errors. Int J Mach Tool Manu. 4012351256
  35. 35. RameshR.MannanM. A.PooA. N.2000Error compensation in machine tools- a review Part II: thermal errors. Int J Mach Tool Manu. 4012571284
  36. 36. RatchevS.GovenderE.NikovS.PhuahK.TsiklosG.2003Force and deflection modelling in milling of low-rigidity complex parts. J Mater Process Tech. 143796801
  37. 37. M. B.da Silva and J. Wallbank (1999Cutting temperature: prediction and measurement methods- a review. J Mater Process Tech. 88195202
  38. 38. MillerM. R.MulhollandG.AndersonC.2003Experimental cutting tool temperature distributions. J Manuf Sci E-T Asme. 125667673
  39. 39. PotdarY. K.ZehnderA. T.2003Measurements and Simulations of Temperature and Deformation Fields in Transient Metal Cutting. Journal of Manufacturing Science and Engineering. 125645655
  40. 40. AbukhshimN. A.MativengaP. T.SheikhM. A.2006Heat generation and temperature prediction in metal cutting: A review and implications for high speed machining. Int J Mach Tool Manu. 46782800
  41. 41. DaviesM. A.UedaT.M’SaoubiR.MullanyB.CookeA. L.2007On the measurement of temperature in material removal processes. Cirp Ann-Manuf Techn. 56581604
  42. 42. LinJ.1995Inverse estimation of the tool-work interface temperature in end milling. International Journal of Machine Tools and Manufacture. 35751760
  43. 43. KwonP.SchiemannT.KountanyaR.2001An inverse estimation scheme to measure steady-state tool-chip interface temperatures using an infrared camera. International Journal of Machine Tools and Manufacture. 4110151030
  44. 44. FangG.ZengP.2005Three-dimensional thermo-elastic-plastic coupled FEM simulations for metal oblique cutting processes. J Mater Process Tech. 1684248
  45. 45. LinJ. M.1995Inverse Estimation of the Tool-Work Interface Temperature in End Milling. Int J Mach Tool Manu. 35751760
  46. 46. ChenM.SunF. H.WangH. L.YuanR. W.QuZ. H.ZhangS. Q.2003Experimental research on the dynamic characteristics of the cutting temperature in the process of high-speed milling. J Mater Process Tech. 138468471
  47. 47. LeshockC. E.ShinY. C.1997Investigation on Cutting Temperature in Turning by a Tool-Work Thermocouple Technique. Journal of Manufacturing Science and Engineering. 119502508
  48. 48. FraserS.AttiaM. H.OsmanM. O. M.1998Modelling, identification and control of thermal deformation of machine tool structures, part 1: Concept of generalized modelling. J Manuf Sci E-T Asme. 120623631
  49. 49. ZhangX.ZhuL.ZhangD.DingH.XiongY.Numericalrobust.optimizationof.spindlespeed.formilling.processwith.uncertaintiesInternational.Journalof.MachineTools.Manufacture2012
  50. 50. HasselmanT. K.ChrostowskiJ. D.1997Effects of product and experimental variability on model verification of automobile structures. P Soc Photo-Opt Ins. 3089612620
  51. 51. BalmesE.1998Predicted variability and differences between tests of a single structure. Imac- Proceedings of the 16th International Modal Analysis Conference, Vols 1 and 2. 3243558564
  52. 52. KurdiM. H.HaftkaR. T.SchmitzT. L.MannB. P.2008A Robust Semi-Analytical Method for Calculating the Response Sensitivity of a Time Delay System. J Vib Acoust. 130.
  53. 53. NocedalJ.WrightS. J.1999Numerical optimization: Springer verlag.
  54. 54. M. H. Richardson and D. L. Formenti, Parameter estimation from frequency response measurements using rational fraction polynomials, 1982.
  55. 55. AltintasY.ShamotoE.LeeP.BudakE.1999Analytical prediction of stability lobes in ball end milling. J Manuf Sci E-T Asme. 121586592


  • Surface location error is defined as the error in the placement of the milling cutter teeth when the surface is generated.

Written By

YongAn Huang, Xiaoming Zhang and Youlun Xiong

Submitted: December 13th, 2011 Published: October 10th, 2012