Cutting dynamics and surface quality
Dániel Bachrathy
Supervisor: Prof. Gábor Stépán
Department of Applied Mechanics
Budapest University of Technology and Economics
A thesis submitted for the degree of
Doctor of Philosophy
2013
To ∞.
A dolgozat magyar nyelvű címe:
Forgácsolás dinamikája és felületi minőség
Nyilatkozat
Alulírott Bachrathy Dániel Sándor kijelentem, hogy ezt a doktori értekezést magam
készítettem és abban csak a megadott forrásokat használtam fel. Minden olyan részt,
amelyet szó szerint, vagy azonos tartalomban, de átfogalmazva más forrásból átvettem, egyértelműen, a forrás megadásával megjelöltem.
Budapest, 2013. agusztus 29.
Bachrathy Dániel Sándor
ii
Acknowledgement
First of all, I would like to thank Prof. Gábor Stépán for his leadership and that
he endeared science to me.
I also thank my colleges Zoltán Dombóvári, Tamás Insperger and the whole staff
of the Department of Applied Mechanics.
I thank Imre Mészáros at Direct-Line Ltd., Jokin Muñoa at the Ideko, Danobat
Technological Centre, János Turi at the University of Dallas at Texas and for Prof.
Berend Denkena and Max Krüger at the Leibnitz University of Hannover who hosted
me for research.
This research was realized in the frames of TÁMOP 4.2.4. A/111-1-2012-0001
“National Excellence Program – Elaborating and operating an inland student and
researcher personal support system”. The project was subsidized by the European
Union and co-financed by the European Social Fund. The research was also supported
by the Hungarian National Science Foundation under grant no. OTKA K101714.
These supports are gratefully acknowledged.
iii
“If I have seen further it is by standing on the shoulders of Giants.”
(1675 - Sir Isaac Newton)
iv
Contents
1 Hard Turning Process
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Dimensionless Dynamical Model . . . . . . . . . . . . . . . .
1.3 Rectangular Step Cutting Force Function . . . . . . . . . . .
1.3.1 Experimental Validation . . . . . . . . . . . . . . . .
1.3.2 Machining Process Parameters . . . . . . . . . . . . .
1.3.3 Surface Topography . . . . . . . . . . . . . . . . . . .
1.4 Surface Pattern Prediction for General Workpiece Geometry
1.5 Continuous Force Function . . . . . . . . . . . . . . . . . . .
1.5.1 Dimensionless Smoothing Parameters . . . . . . . . .
1.5.2 Analysis of Smoothing Parameters . . . . . . . . . .
1.5.3 Experimental Validation . . . . . . . . . . . . . . . .
2 Dynamics of Milling
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Surface Quality of Milled Surfaces . . . . . . . . . . . . . . .
2.2.1 Mechanical Model . . . . . . . . . . . . . . . . . . . .
2.2.1.1 Straight Fluted Tool . . . . . . . . . . . . .
2.2.1.2 Helical Fluted Tool . . . . . . . . . . . . . .
2.2.2 Stationary and Transient Motion of the Milling Tool
2.2.3 Forced Vibration and Machined Surface . . . . . . .
2.2.4 Surface Quality Parameters . . . . . . . . . . . . . .
2.2.5 Effect of Machining Process Parameters on Profile . .
2.2.5.1 Real-case Mechanical Parameters . . . . . .
2.2.5.2 Effect of Machining Process Parameters . .
2.2.5.3 Fourier transformation of the cutting force .
2.2.6 Experimental setup . . . . . . . . . . . . . . . . . . .
2.2.7 Results and Experimental Validation . . . . . . . . .
2.2.8 Multiple Degree of Freedom System . . . . . . . . . .
2.3 Surface Reconstruction for Combined FRF and FEM model
2.3.1 Development of a Dynamic Tool Model . . . . . . . .
2.3.1.1 Dynamic Model of the Spindle Part . . . . .
2.3.1.2 Dynamic Model of the Tool Part . . . . . .
2.3.1.3 Global System Matrices . . . . . . . . . . .
2.3.1.4 Forced Motion . . . . . . . . . . . . . . . .
v
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
2
3
4
6
7
7
11
13
14
15
16
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
19
19
22
22
23
24
26
27
29
30
30
31
33
35
35
37
43
43
44
45
47
48
2.3.2
2.3.3
2.4
2.5
Generation of Surface Topography in Peripheral Milling . . . .
Experimental Verification . . . . . . . . . . . . . . . . . . . .
2.3.3.1 Description of the Experimental Setup . . . . . . . .
2.3.3.2 Reconstruction of Surface Topographies . . . . . . .
Effects of Large Amplitude Vibrations on Milling Stability . . . . . .
2.4.1 Velocity Dependent Dynamic Chip Thickness in Milling Processes
2.4.2 State-Dependent Delay Model . . . . . . . . . . . . . . . . . .
2.4.2.1 Improved Chip Thickness Modelling . . . . . . . . .
2.4.2.2 Improved Screen Function Modelling . . . . . . . . .
2.4.2.3 Periodic Solution . . . . . . . . . . . . . . . . . . . .
2.4.2.4 Linearization Around the Periodic Solution . . . . .
2.4.2.5 Stability by the Semi-Discretization . . . . . . . . . .
2.4.2.6 Real-Case Parameters . . . . . . . . . . . . . . . . .
2.4.2.7 Stability Chart . . . . . . . . . . . . . . . . . . . . .
2.4.2.8 Fold Bifurcation . . . . . . . . . . . . . . . . . . . .
2.4.2.9 Limitations of the Model . . . . . . . . . . . . . . . .
Process Damping in Milling . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Process Damping . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Governing Equation with Process Damping . . . . . . . . . . .
2.5.3 Short-Delay Representation . . . . . . . . . . . . . . . . . . .
2.5.4 Periodic Solution . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.5 Linearization and Stability . . . . . . . . . . . . . . . . . . . .
2.5.6 Case Study: Effect of Flank Wear Length on Stability . . . . .
3 Multi-Dimensional Bisection Method
3.1 Introduction . . . . . . . . . . . . . . . . . . . . .
3.2 Bisection Method . . . . . . . . . . . . . . . . . .
3.2.1 Numerous Nearby Roots . . . . . . . . . .
3.2.2 Generalization for Higher Dimensions . . .
3.2.3 n-simplexes and n-cubes . . . . . . . . . .
3.2.4 Selection of Bracketing n-cubes . . . . . .
3.2.4.1 Safe Selection . . . . . . . . . . .
3.2.4.2 Secant Based Selection . . . . . .
3.2.4.3 Fitting Hyper-Planes . . . . . . .
3.2.4.4 Comments on Selection Methods
3.2.5 Continuation at Neighboring n-cubes . . .
3.3 Numerical Implementation . . . . . . . . . . . . .
3.4 Efficiency Number . . . . . . . . . . . . . . . . .
Reference
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
49
51
51
53
57
57
61
63
65
66
67
70
70
70
71
73
78
78
79
80
81
81
82
85
85
86
87
87
88
88
88
89
89
89
89
90
92
97
vi
Outline of the Thesis
The dissertation summarizes our results on the exploration of the correlation between
machined surface quality and the dynamics of cutting processes. The aim of this
exploration is to predict the machined surface properties in hard turning and in
general milling processes where machine tool vibrations are likely to appear at cutting
parameters optimized for high material removal rates (HMRR).
Turning and milling processes are widely used non-abrasive procedures in the
industry. When the HMRR is to be achieved, the cutting parameters are pushed to
the limits. The most critical and the less predictable limit is related to the appearance
of the so-called chatter vibrations. When these high amplitude self-excited oscillations
show up, the machined surface quality gets intolerable, in worst cases even the tool
and the machine tool structure is damaged. This is why the precise prediction of
machine tool vibration is important in practice.
Even if chatter is avoided, the good surface quality is not guaranteed for HMRR
optimized cutting parameters, since these are just in those domains where resonant
forced vibrations occur. While these vibrations are less dangerous, they may still have
a negative effect on surface quality when high precision is required. Special analytical,
simulation and optimization techniques are required to model the machined surface
quality and the dynamic machining process and to find reliable optimal machining
parameters.
The dissertation is divided into Chapters based on three main research topics.
In Chapter 1, the machined workpiece surface is modelled in case of hard turning processes, optimization criteria are defined and optimal cutting parameters are
selected.
Chapter 2 deals with the milling processes in several aspects. First, optimal
machining parameters are determined based on the simulated workpice quality for two
different mechanical models in Sec. 2.2 and Sec. 2.3. Second, the stability properties
of the milling processes are analysed considering large amplitude resonant vibrations
and the so-called process damping phenomenon in Sec. 2.4 and Sec. 2.5, respectively.
Finally, in Chapter 3, we generalize the bisection method to higher dimensions in
order to speed up our computationally intense numerical algorithms, like the construction of stability charts. The method is formulated in a general way to make it usable
in a wide range of engineering problems formed mathematically as root-findings.
1
Chapter 1
Hard Turning Process
1.1
Introduction
In the high precision hard turning processes, the prescribed accuracy and surface
roughness are in the range of micrometers [1]. In case of continuous or uninterrupted
surfaces, these can be reached by choosing the right parameters of the cutting process
[2, 3, 4].
The interruptions of the surface like grooves or holes disturb the stable stationary
cutting process. In high precision hard turning, interruptions raise problems in two
aspects. The first aspect is related to the magnitude of the thrust (passive) component
of the cutting force. In conventional turning processes, the thrust force component has
a negligible magnitude, while this component is dominant in hard turning processes
[5].
The other aspect of the difficulties is that the cutting force decreases to zero at
the edge of the interruption for a short time period and starts cutting again as it
intrudes into the material where the cutting force may reach its maximal value. The
fluctuation of the thrust force creates vibrations perpendicular to the surface, which
causes dimensional and shape error; these errors can be larger than the specified
tolerances. The force fluctuation also creates additional load on the tool-tip, which
reduces the life-span of the tool [6].
Our goals are partly to determine the optimal technological parameters of hard
turning considering the dynamical properties of the machine tool, and partly to predict the machined surface quality in these cases. First, the cutting force is modelled
by a discontinuous (rectangular step) time-function and a simple easy-to-use formula
of the optimal cutting velocities is derived where the surface error is minimal.
The surface pattern is also calculated for a workpiece of arbitrary geometry.
Then, instead of the rectangular step function, an improved time-continuous cutting force model is introduced with the help of appropriately chosen smoothing parameters. The effects of the dimensionless smoothing parameters on the surface waviness
and on the optimal cutting velocity are analysed. Cutting tests are presented for
the validation of the theoretical results by means of the measurements of tool-tip
acceleration.
2
e
groove
F
R
a0
tool
workpiece
vf
0 z
k
f
r
W
m
c
Fp
vf
W
Figure 1.1: Derivation of the single DoF dynamical model and the main technological parameters.
1.2
Dimensionless Dynamical Model
A mechanical model is constructed in order to analyse the properties of the machined
surface cut by the vibrating tool. Accordingly, only the tool motion perpendicular to
the surface is relevant, since this is the vibration that leaves its trace on the machined
surface during the cutting process. Thus, a single degree of freedom (DoF) damped
oscillator is chosen as a dynamical model which is excited by the cutting force which
is periodic with the time period of the workpiece revolution. The parameters of the
model are shown in Fig. 1.1.
The governing equation of the tool-tip motion in the z direction is given by
mz̈(t) + cż(t) + kz(t) = Fp (t),
(1.1)
where Fp is the thrust (passive) force component and the mass m, damping c and the
stiffness k are the modal parameters. The values of these parameters determine the
dominant (undamped) natural frequency ωn , which can be identified by experimental
modal analysis of the given machine tool [7, 8]. To create a general description,
the governing equation is transformed to dimensionless form by using a new general
coordinate
u(t) =
z(t)
,
zstat
(1.2)
where zstat is the static deformation created by a constant (uninterrupted) cutting
force. The dimensionless form of the governing equation is given by
u00 (τ ) + 2ξu0 (τ ) + u(τ ) = Φ(τ ),
(1.3)
where ξ = c/(2mωn ) is the damping ratio, Φ = Fp /(kzstat ) is the dimensionless thrust
force and
p τ = tωn is the dimensionless time. The (undamped) natural frequency
ωn = k/m = 2π/Tn and Tn is the time period of the (undamped) free oscillation.
Note, that the time derivation is transformed into
∂
∂
= ωn .
∂t
∂τ
3
(1.4)
thrust force, Fp
chip width
v
xstat
Wt
tcut
t
tfly
surface
interruption
motion of the cutting
Figure 1.2: Representation of the oscillation of the tool in the cutting and flying phases.
1.3
Rectangular Step Cutting Force Function
According to the so-called surface regeneration effect, the cutting force is a function
of the depth of cut, which, strictly speaking, depends on the current and the previous
position of the tool [9]. This may lead to stability problems. The mathematical
model of the interrupted turning process leads to a delay-differential equation with
time-periodic coefficients. The so-called stability chart of our system were analysed
by the semi-discretization method developed, described and tested by Insperger and
Stépán [10, 11, 12, 13]. Detailed descriptions of the stability of the interrupted turning
process can be found in [14]. In the range of the available technological parameters of
the high precision hard turning tool in question, the theoretical regenerative models
did not predict unstable parameter regions and the chatter vibrations did not show
up during the measurements, either.
Accordingly, if the surface regeneration effect is neglected and the forced vibration
amplitude is much smaller than the prescribed stationer chip thickness (a0 ), then the
(dimensionless) thrust force can be approximated by a constant value (Φ = 1) during
the cutting phase. Of course, the cutting force is zero (Φ = 0) during the flying phase
when the tool has no contact with the material, the oscillation is free. Consequently,
the interrupted turning process can be divided into two phases: cutting and flying
phases where the cutting force switches on and off with the dimensionless durations
τcut and τfly , respectively. These two phases are represented in Fig. 1.2.
The vibration of the tool-tip is superposed to the desired machined surface in the
cutting phase, which leads to surface error. The dimensionless surface waviness Wt
can be computed from the tool motion in the cutting phase as follows
Wt = max(u(τ )) − min(u(τ )),
τfly ≤ τ ≤ τfly + τcut
(1.5)
where the dimensionless time τ is set to zero at the exit point. In this model, the
surface waviness depends only on the damping ratio ξ and the dimensionless flying
time τf ly . The experiments show that the vibrations created by the force fluctuations
die out much sooner than the tool-tip exits the material again [15, 16]. Thus, it
is assumed that there is no oscillation at the end of the cutting phase, that is, the
4
Wt - dimensionless
surface waviness [1]
4
x - damping fac.[1]
0.000
0.010
0.041
0.092
0.163
0.255
0.367
0.500
1.000
3
2
1
0
0
2p
4p
6p
8p
10p
12p
tfly - length of free oscillations [1]
14p
16p
Figure 1.3: Dimensionless surface waviness as a function of the damping ratio ξ and the length of
the free oscillation τfly .
transient motion dies out during the τcut period and the rotational periods can be
analysed with separate initial conditions
u(0) = 1,
u̇(0) = 0,
(1.6)
at the beginning of the flying phase. The motion law of the free oscillation and then
the initial condition of the cutting phase can be calculated in closed form. A similar,
but much more complicated closed form equation can be given for the motion law in
the cutting phase. Thus, it is possible to determine the max(u(τ )), min(u(τ )) and
Wt itself in Eq. (1.5) in closed form, however, the large size of the equation makes it
difficult to check the influence of the relevant parameters.
The solution can also be generated easily by means of numerical integration of
Eq. (1.3) with the initial conditions Eq. (1.6). The effects of the parameters on the
calculated dimensionless surface waviness are plotted in Fig. 1.3.
If the damping is negligible and the dimensionless flying time is
opt
τfly
= j2π
j = 1, 2, 3, . . . ,
(1.7)
that is, the number of free oscillation periods is integer during flying phase, then the
surface waviness is almost zero (see: Fig. 1.3). In this case, Eq. (1.7) defines the
optimal flying time. It is easy to show that, the surface waviness is 4zstat in the worst
case.
As Figure 1.3 shows, increasing damping generally decreases the surface waviness,
that is, improves surface quality. However, this is just the opposite for the specific
optimal lengths of the flying times, especially for the dimensionless flying time at 2π.
In a real machine tool, the damping is usually small and cannot be adjusted much,
but the cutting velocity can easily be set according to the optimal flying time. Since
the length of the interruption is the width e of the groove (see Fig. 1.1), the flying
time is determined as
e ω
eωn
eωn
n
≈
=
,
(1.8)
τfly = 2 arcsin
2r Ω
rΩ
v
5
y [mm]
20
a.)
4
b.)
3
groove
10
2
0
1
-10
vopt=241,48 [m/min]
n=1800 [rpm]
-20
0
-1
-30
-20
0
20
x [mm]
-20
0
z - hight of the surface [mm]
30
20
Figure 1.4: The vibration of the tool and the waviness of the surface. Black lines denote the boundary
of the groove and the work-piece in case of a) constant spindle speed n b) varying spindle speed
providing constant optimal cutting speed vopt . Parameters: m=0.056[kg], c=21[Ns/m], k=20[N/µm],
e=4[mm], R=30 [mm].
where r is the radial position of the cutting edge, v is the cutting speed (the peripheral
velocity) and Ω is the spindle speed. The error of the approximation in Eq. (1.8) is
1.178% in case of r/e=2, while it is only 0.067% in case of r/e=4.
Using this approximation in Eq. (1.8), also Eq. (1.7) and Eq. (1.4), the optimal
cutting speed can be determined as
vopt =
e
jTn
j = 1, 2, 3, . . . .
(1.9)
It can be seen in Eq. (1.9), that in case of a groove where the width of interruption
e is constant along radius r along the workpiece, constant cutting speed can be used
to reach optimal surface quality. To determine the optimal cutting speed, only the
time period of free oscillation Tn must be used, which can be identified by a single
impact test procedure.
In Fig. 1.4, the cut surface is calculated for two cases [17]. First, constant spindle
speed was used (Fig. 1.4a) when the cutting velocity varies during the cutting process
with the radial position of the tool. It can be seen that for some radii, the cutting
speed becomes close to one of the optimal cutting speeds, that is, an integer number
of oscillations takes place in the flying period, which leads to good surface properties.
Between these optimal points along the radial direction, high surface errors take place.
In the second case (Fig. 1.4b), an optimal cutting velocity was used. In this case, the
surface quality is very good at all radii.
1.3.1
Experimental Validation
To validate the identified optimal cutting speeds Eq. (1.9), cutting tests were carried
out on a workpiece shown in Fig. 1.5a. Uddeholm Rigor high speed steel material
was selected for the cylindrical specimens. This material was hardened in its whole
cross section. Experiments were carried out on a Slantbed-Mikroturn 50 CNC lathe
6
a.)
b.)
Figure 1.5: Geometric parameters of the workpiece and the measurements setup.
shown in Fig. 1.5b. The cutting tool was used with a CCMW09T304S-L0-B CBN050C
SECO insert. During, the measurements the vibrations were detected by B&K 4397
accelerometers and analysed by the PULSE data acquisition system. After the double integration of the acceleration signal, the position signal contains a drift error
that was compensated by a filtering method similar to the Zero Velocity Compensation published in [18] and [19]. The frequency response function of the tool-tip was
measured by an impact test [7, 8]. The system has a dominant natural frequency
around ωn =2240 Hz. The damping ration is approximately ξ = 0.05. In the current
measurements, the optimal cutting speeds are vopt =537.6; 268.8; 179.2; 134.4; 107.5;
89.6; ... m/min.
1.3.2
Machining Process Parameters
During the cutting tests, constant feed ratio 25 µm/rev, constant depth of cut 40
µm and a set of constant spindle speeds n = (1232; 1155; 1078; 1001; 924; 821)
rpm were used. Consequently, the cutting speed changed continuously during the
measurements as a function of the radial position of the tool. Based on the results
shown in Fig. 1.4a, varying vibration levels are expected in the measurements during
the cutting phase at different radii. For example: if the spindle speed is n=1001
rpm, smaller vibration is expected around the radii r = 21.6 mm and 28.5 mm (see
Fig. 1.6).
1.3.3
Surface Topography
We assume that the machined surface is created exactly by the vibrating cutting
tool. The machined surface can be represented by transforming the time history of
the tool position (more exactly, its coordinate z) to a spiral in the x − y plane. This
spiral is the path of the motion of the tool relative to the workpiece. The surfaces
generated this way are shown in Fig. 1.7. A photo is also shown where the waviness
of the machined surface can be traced by means of the light reflection. The photo
7
vopt
n
r - radial position of the tool [mm]
Figure 1.6: Optimal cutting speeds (horizontal thick lines) and the applied cutting speeds of the
measurements as a function of the radial position (slanting lines)
and the corresponding generated surface show the same pattern of surface waviness.
Note, that even the small surface error caused by a previous grinding process can be
observed in the upper right corner of the groove at all the generated surfaces.
In Fig. 1.8 a typical signal of the measured tool motion is plotted around the
groove. We can see that the flying phase and the cutting phase are well separated
and have similar pattern as represented in Fig. 1.2. The surface waviness Wt is defined
for the tangential direction as a difference of the maximal and the minimal position
of the tool during the cutting phase (see Fig. 1.8 and Eq. (1.5)), which is, of course,
a function of the radial position of the tool. In Fig. 1.9, these functions are plotted
for the measured and the theoretical surface waviness. To be able to compare the
theoretical and the measured surface waviness, we have to define the static deformation of the tool (see Eq. (1.2)). We were unable to measure this directly, so the static
deformation was chosen to fit to the measurement results (zstat =1.48 µm), which has
no effect for the loci of the radial positions with optimal cutting speeds, that is with
minimal surface waviness. Very good correlation is found between the measurements
and the analytical results. However, in case of n=821 rpm, the optimal points are
slightly shifted to the right. Another slight deviation of theory and measurements can
be observed in Fig. 1.8 where the patterns of the transient vibrations do not present
the exponential decay predicted for the hypothetical rectangular step cutting force
function. Further analysis of these phenomena and a smooth continuous force model
are given in Section 1.5.
8
9
2
0
-2
-4
z - simulated profie hight [mm]
4
-6
Figure 1.7: Simulated surface pattern based on the acceleration measurement and the picture of the
real machined workpiece
10
z - position of the tooltip [mm]
flying phase
highest peak
deepest valley
transient vibration
cutting phase
cutting phase
t - time [ms]
Wt - surface waviness [mm]
Figure 1.8: The measured tool-tip position.
5
4
3
2
1
0
5
4
3
2
1
0
5
4
3
2
1
0
20
22
24
26
28
30
20
22
24
26
28
r - radial position [mm]
Figure 1.9: The theoretical and the simulated surface waviness.
30
Set
M
I1
I2
I3
I4
I5
conditions
|x| < 40 & |y| < 30
x > 45 & |y| < 12.5
x < 30 & |y| < 5
x2 + (y − 20)2 < 102
x2 + y 2 < 52
x2 + (y + 20)2 < 52
Table 1.1: Definition of the set of the initial geometry M (x, y) of the workpiece and the interruption
sets Ij (in (mm) unit).
1.4
Surface Pattern Prediction for General Workpiece Geometry
The optimal cutting speed (Eq. (1.9)) can be used only in symmetrically placed
’groove-like’ interruptions. In other cases, the arc length of the interruption is not
constant at every radii. In case of a cylindrical hole, its size is changing heavily,
and the cutting speed should be changed accordingly as given in Eq. (1.9). Further
problem occurs if the workpiece hase two grooves with different widths. In this case,
the cutting speed should be changed within one rotation, which is hardly possible
considering the high inertia of the spindle in a real machine tool. In case of workpieces
with general geometry like these, the applications of the optimal cutting velocities
are not feasible. Still, the prediction of the machined surface pattern can support the
appropriate choice of the cutting parameters.
The predicted surface pattern can be computed for a general workpiece geometry
by means of the simulation of the governing equation Eq. (1.1) for the whole turning
process. In these tasks, the time history of the thrust force component must be
computed. First, the geometry of the workpiece is described in the x − y plane by
means of Boolean operations. The set W (x, y) of the actual geometry of the workpiece
is given as the initial geometry M (x, y) of the workpiece then subtracting the unions
of interruption sets Ij :
W = M \ (I1 ∪ I2 ∪ . . . In−1 ∪ In )
j = 1, 2, 3, . . . , N.
(1.10)
This type of workpiece representation can easily be updated after any cutting
production step. Table 1.1 describes a rectangular workpiece (shown in Fig. 1.10a)
with holes and groove.
The same assumption is used as in Section 1.3, that is, the thrust force is constant
during the cutting phase and zero otherwise (flying phase):
0
Fp
if (x(t), y(t)) ∈ W
Fp (t) =
(1.11)
0
otherwise,
where (x, y) is the current position of the tool-tip which can be computed from the
current radial position of the tool and the angular position of the workpiece. The
11
12
120
y
W
I3
R5
20
20
10
I4
I2
15
x R5
I1
80
30
R10
25
(a)
I5
40
3
y [mm]
20
2
0
1
-20
0
-40
-1
-80 -60 -40 -20
0
20
x [mm]
40
60
z - hight of the surface [mm]
(b)
80
Figure 1.10: Sketch of workpiece (a) and the simulated workpiece profile and tool motion (b).
Parameters: m=0.034 kg, c=235 Ns/m, k=10 N/µm, Fp0 =40 N, n=5000 rpm.
(a)
(b)
-2
0
2
4
z- hight of the surface [mm]
4
8
Figure 1.11: Picture of the machined crown gear (a) and the simulated workpiece profile (b). Parameters: m=0.632 kg , c=100.5 Ns/m, k=10 N/µm, Fp0 =47.4 N, v=120 m/s.
motion of the tool-tip z(t) is determined by means of the numerical integration of
the governing equation Eq. (1.1), and the surface is generated in the same way as
described in Subsection 1.3.3. The surface waviness, can also be computed from the
tool-motion (see Eq. (1.5)). The predicted surface and the corresponding tool motion
are shown in Fig. 1.10b.
The calculation is tested for a case-study at RABA Automotive Group. Figure 1.11 shows the different views of a crown gear together with its simulated surface
pattern. The modal parameters of the machining system were estimated from the
geometric data of the machined surface pattern and the cutting speed. It can be seen
in Figure 1.11 that the numerically generated and the machined surface patterns are
the same, which validates our model and method.
In the above calculations, a rectangular step force function is used. It was shown,
however, that there are slight deviations in the measured and the computed tool
vibrations (see Fig. 1.9). One reason for this could be the still approximate description
of the thrust force that may start decreasing before it leaves the workpiece material at
a groove (or hole) and which can build up continuously to (temporarily) even higher
values than the stationary thrust force when it enters the workpiece material again. In
the next section an improved continuous force model is introduced by using smoothing
parameters. The effects of the smoothing parameters on the surface waviness and on
the optimal cutting speeds are also given.
1.5
Continuous Force Function
The formula of the optimal cutting speeds Eq. (1.9) is valid only in case of a rectangular step force function. In a real machining process the thrust force in the cutting
phase starts decreasing before reaching the boundary of the material at the interruptions. As the tool reaches the material at the end of the flying phase, the cutting force
13
f(t)
1.4
a
K
1.2
1.0
{
d
c
0.8
0.6
0.4
b
0.2
-2
0
2
4
6
t
Figure 1.12: Representation of the continuous smoothed dimensionless thrust force function. First
phase: decreasing thrust force before exit (a). Second phase: non-cutting or flying (b). Third phase:
increasing cutting force (c) and decay of the increased cutting force coefficient (d) at the intrusion.
Parameters: εexit = 0.4, τfly = 3, εenter = 0.2, K = 0.4, εimp = 1.
starts increasing continuously. The maximum of the thrust force is probably reached
just after the intrusion due to the low temperature of the material at the entrance
point, where the stationary temperature field around the tool-tip has not been build
up yet [20]. In this section, a corresponding model of the cutting force is introduced.
1.5.1
Dimensionless Smoothing Parameters
The improved model of the cutting force function consists of three phases around the
interruptions:
1 − eτ /εexit
0
Φ=
1 − e(τ −τfly )/εenter 1 + Ke(τ −τfly )/εimp
τ ≤0
0 < τ ≤ τf ly
τ > τfly ,
(1.12)
see also Fig. 1.12. The first phase describes the force before the tool leaves the
surface, by means of an exponentially decreasing term with time constant εexit (see
Fig. 1.12a). In the second phase, the cutting force is zero during the free flight, where
no contact takes place between the tool and the workpiece material (see Fig. 1.12b).
The last phase is a product of two exponential functions: the first describes the fast
increase of the thrust force after the intrusion by using the time constant εenter (see
Fig. 1.12c), while the second determines the larger cutting force coefficient due to
the low temperature of the material at the entrance point of the interruption (see
Fig. 1.12d). This subsystem has two parameters: the time constant εimp and the
maximal value of the increase K.
14
1.5.2
Analysis of Smoothing Parameters
The improved model of the force function requires four additional so-called smoothing parameters. Eq. (1.12) describes a rectangular step function if the additional
smoothing parameters are set to zero. According to the preliminary experiments, the
values of these additional parameters should be close to zero, indeed. This is why it is
assumed that the effects of the parameters can be analysed independently except the
parameters K and εimp . The influence of the smoothing parameters on the dimensionless surface waviness Eq. (1.5) and on the optimal cutting speeds Eq. (1.9) are
calculated from the governing equation Eq. (1.3) that includes the improved thrust
force model Eq. (1.12).
First, the effect of the time constant εexit is studied. The optimal flying time can
be computed by means of the solution of Eq. (1.3) in closed form as a function of the
time constant εexit :
εexit
opt
j = 1, 2, 3, . . . .
(1.13)
τfly = 2jπ − arctan
1 + 2ξεexit
If the damping is small enough and εexit << 1, then Eq. (1.13) can be simplified
to
opt
τfly
≈ 2jπ − εexit
j = 1, 2, 3, . . . .
(1.14)
Figure 1.13a) presents surface waviness calculated numerically for different time
constants εexit as a function of the flying time while the damping ratio ξ=0.1 is fixed.
It can be seen that the maximal and minimal surface waviness values move closer
to an average one, and the minima of the surface waviness are shifted to the left as
predicted by the derived formulae Eq. (1.13) and Eq. (1.14).
The time constant enter has similar effect on the optimal flying time
εenter
opt
.
j = 1, 2, 3, . . . ,
(1.15)
τfly = 2jπ − arctan
1 + 2ξεenter
while it improves the surface quality uniformly for any flying time, as it is shown in
Fig. 1.13b.
The effect of the parameters εimp and K cannot be analyzed separately: if any of
them is zero then the other has no effect at all. First, εimp =0.5 is fixed and K is varied
(Fig. 1.13c) then K = 0.5 is fixed and εimp is varied (Fig. 1.13d). The numerical results
show that the increase of the parameter K or εimp leads to the uniform increase of
the dimensionless surface waviness while these parameters have no significant effect
(<2%) on the optimal flying time. If these parameters are dropped then the final
form of the optimal flying time can be constructed from Eq. (1.13), Eq. (1.15) and
Eq. (1.9) as follows:
vopt =
e
Tn j2π − arctan
2π
εenter
1+2ξεenter
15
− arctan
εenter
1+2ξεenter
(1.16)
Wt - dimensionless
surface waviness [1]
1.5
1
0.5
0
Wt - dimensionless
surface waviness [1]
(a)
2
eimp=0, K=0, eenter=0
0
2p
4p
6p
tfly - length of free oscillations [1]
eimp=0, K=0, eenter=0
2p
1.5
1
0.5
2p
4p
6p
tfly - length of free oscillations [1]
6p
K [1]
eimp [1]
0.0
0.5
1.0
1.5
2.0
K=0.5, eexit=0, eenter=0
0
8p
eenter [1]
0.0
0.2
0.4
0.6
0.8
1.0
8p
(d)
eimp=0.5, eexit=0, eenter=0
0
4p
tfly - length of free oscillations [1]
0.0
0.2
0.3
0.6
0.8
1.0
2
0
(b)
0
8p
(c)
2.5
eexit [1]
0.0
0.4
0.8
1.2
1.6
2.0
2.4
2p
4p
6p
tfly - length of free oscillations [1]
8p
Figure 1.13: Dimensionless surface waviness as a function of the dimensionless time-period of the
free oscillation τfly and time constant εexit (a), the time constant εenter (b), the increase parameter
K (c), and the time constant εimp (d). (ξ = 0.1)
vopt ≈
e
2π
Tn j2π − εexit − εenter
j = 1, 2, 3, . . . .
(1.17)
If the dimensionless parameters are measured for a given material and tool geometry, the optimal cutting velocity can be computed for a given geometry of interruption
and for the given dynamic parameters of the machine tool.
1.5.3
Experimental Validation
The same measurement setup is used as described in Sec. 1.3.1 with a tool holder
of smaller tool overhang. The signal noise was reduced by averaging 10 periods of
the measured signals. In the measurements, constant cutting parameters were used:
spindle speed 1000 rpm, depth of cut 40 µm and feed rate 25 µm/rev.
In the validation process, all the system parameters (natural frequency ωn , damping ratio ξ and static deformation zstat ) and all the smoothing parameters of the force
function (flying time τfly , time constants εexit , εenter , εimp and increase parameter K)
were computed by the Least Square Method applied for the time histories of the
measured and the computed tool-tip positions. The resulting parameters are shown
in Table 1.2, and the measured as well as the computed tool positions are plotted in
Fig. 1.14.
The measured tool motion involves an additional frequency, which refers to a
system having at least two degrees of freedom. Since the amplitude of this vibration
is quite small, the 1 DoF model given in Eq. (1.1) is satisfactory. In Table 1, we
can see that the time constants εexit and εenter have relatively large values, so the
16
-7
2 x 10
measured
tool motion
0
0
-0.2 dimensionless
cutting force
-0.4
-2
-4
computed
tool motion
-6
-0.6
-8
-10
0
-0.8
0.01 0.02 0.03 0.04 0.05 0.06
0.026
t - measured time [s]
0.027
0.028
0.029
t - measured time [s]
Figure 1.14: The measured (continuous line) and computed (dotted line) positions and the dimensionless cutting force function (dashed line).
ωn
(Hz)
2729
ξ
(1)
0.184
zstat
(µm)
0.625
τfly
(1)
20.65
εexit
(1)
1.511
εenter
(1)
1.953
εimp
(1)
1.70
K
(1)
0.054
Table 1.2: Parameters determined from the measurements.
improved formula Eq. (1.16) of the cutting speed should be used instead of Eq. (1.9)
derived for rectangular step function. The 0.054 value of parameter K shows that
the entrance cutting force is not significantly higher than the stationary cutting force
of un-interrupted turning, so this parameter can be dropped from the improved force
model, which also means that the parameter εimp has no relevance, either.
While the experimental validation confirmed the use of the improved continuous
cutting force model Eq. (1.12), it also proved that it can be simplified to the form:
τ ≤0
1 − eτ /εexit
0
0 < τ ≤ τf ly
Φ=
(1.18)
(τ −τfly )/εenter
1−e
τ > τfly ,
for hard turning processes and the improved optimal flying time can be given by
Eq. (1.16).
New Results
Thesis 1
Interrupted turning processes are considered on a lathe having a well separated dominant vibration mode with angular natural frequency ωn . When a discontinuous
(rectangular step) cutting force function is assumed, the simple formula
vopt =
e
jTn
j = 1, 2, 3, . . .
is derived for the optimal cutting speeds vopt , where the surface waviness created by
the transient vibration after the interruption is minimized. These optimal cutting
17
speeds depend on the width e of the groove, and on the time period Tn = 2π/ωn .
In case of facing, the machined surface pattern of a workpiece with general interrupted geometry is accurately modelled by the transformation of the numerical
simulation results from the time domain to the tool path relative to the workpiece.
A more accurate model is also constructed to determine the optimal cutting speeds
by means of improved (empirical) continuous cutting force functions. It is shown
that the influence of the continuous decrease and increase of the cutting force at the
entrance and exit of the interruption can be described by the time constants εexit and
εenter , which significantly improve the model accuracy. In the meantime, the effect of
a possible impact-like increase of the cutting force at the entrance has no relevance
regarding the optimal cutting speed, which has the improved formula
vopt =
e
Tn j2π − arctan
2π
εenter
1+2ξεenter
− arctan
εenter
1+2ξεenter
j = 1, 2, 3, . . .
with damping ratio ξ.
All these results are confirmed in an industrial case study and by test measurements for hard turning processes.
Related publications:
[21, 15, 16, 22, 23, 24, 25]
18
Chapter 2
Dynamics of Milling
2.1
Introduction
In milling processes, vibrations are always generated by the rotating milling tool that
provide periodic material removal. This forced vibration is occasionally combined
with the so-called chatter vibration originated in the possible instability of the periodic
milling process. These are the two main sources of vibrations in milling [26].
The levels of the stable periodic forced vibrations are usually small and create
negligible surface errors due to the high dynamic rigidity of the machine tool structure.
However, large periodic forced vibrations are likely to occur just for the most efficient
near-resonant spindle speeds, which are related to low-damped modes (for instance,
in case of thin walled workpieces) [27, 28, 29, 30]. In these cases, the surface integrity
is usually still acceptable, but the surface location error, the surface roughness and
waviness may be large and intolerable. These forced near-resonant vibrations also
create an additional undesired load on the machine tool. These problems result
further limitation to the increase of the material removal rates, and the use of the
cutting parameters in the efficient stable parameter pockets requires additional careful
analysis.
In Chapter 2.2, the effect of these large amplitude forced vibrations are analysed
with a special attention on surface quality. First, three specific parameters are introduced for engineering use to characterize the spatial machined surfaces created by
helical fluted milling tools. Second, it is proved analytically and verified by measurements for helical milling tools that the resonant vibrations caused by the periodic
cutting force can be minimized with the help of appropriately chosen axial immersions.
While in Chapter 2.2, the milling tool is modelled as a rigid body moving in the
plane perpendicular to its axis of rotation, long finger-like milling tools can have large
deformations and the above model is not satisfactory. In Chapter 2.3, an improved
dynamical model of the milling machine is constructed that consists of three main
subsystems. These include the finite element beam model of the finger-like milling
tool, the model of the machine tool structure based on its frequency response measurement, and their experimentally verified dynamic coupling model. The workpiece-side
19
is considered to be rigid. The proposed tool model is able to reconstruct the accurate
shape and roughness of machined surfaces in many industrial applications. For numerical simulations, this machine tool model is subjected to the cutting force which
is measured online during the milling process. The developed method and model
are verified by comparing the numerically reconstructed and the measured surface
topographies for both stationary and transient vibrations.
The second type of vibrations in milling represent an essential limit for the increase
of the material removal rate; this is called chatter caused by the surface regeneration
effect [31, 32, 33, 26]. The corresponding periodic, quasi-periodic or chaotic vibrations
appear when stationary milling process loses its stability [9, 34, 35, 36], in other
words, when the stationary periodic forced vibration becomes unstable. These types
of vibrations destroy the surface integrity and can create heavy load on the machine
tool. The so-called stability chart presents and predicts the chatter-free parameter
domains of the milling process traditionally in the plane of the spindle speed and the
axial immersion.
Most of the mathematical models neglect some essential non-linear effects due to
the fact that they use constant time delay only that is inversely proportional to the
fixed cutting speed. Since the forced vibrations may become very large in the nearresonant parameter domains, the path of the individual cutting edges will strongly
deviate from the approximate circular paths. In [37], an accurate model of the chip
removal of milling is described for these large amplitude vibrations. This model leads
to a time- and state-dependent delay model, which is essentially non-linear even in
the delay parameter. The governing equation of the process and the linearization
around its otherwise unknown periodic solution are described in [37] and this would
leads to a somewhat modified stability chart.
In Chapter 2.4, the governing equation of milling processes is generalized with
the help of the above mentioned accurate chip thickness derivation resulting in a
state-dependent delay model. The periodic motions of this non-linear system and
its stability are calculated. It is observed that the influence of the state-dependent
delay on linear stability can be significant in the near-resonant parameter domains.
The existence of an unusual fold bifurcation is shown where a kind of hysteresis
phenomenon may appear between two different stable periodic motions.
Experimental verification of theoretically predicted stability charts for high spindle
speeds is given in [38] and [39], where good correlation is found between the measured
and the computed stability boundaries. However, the theoretical predictions are too
conservative in the low-speed domain, that is, the experiments present stable cutting
for larger axial immersions compared to the critical ones calculated from the model
[40, 41, 42, 43, 44]. The process damping effect is often used to improve these models.
Additional damping proportional to the time delay is added to the structural damping
to characterize the contact force at the tiny region between the flank face and the
workpiece [40, 45, 46, 41, 44].
In Chapter 2.5, the above mentioned classical process damping is used to improve
the milling model, which leads to a system with a non-smooth (unilateral) damp20
ing. The improved model predicts the expected higher stability limits at low spindle
speeds, but it describes a new intricate pattern of the lobe-structure. It is also shown,
that the mathematically correct description of the widely used classical process damping model is the same as the one based on the so-called short regenerative effect that
contains distributed delay terms [9]. It is proven, that the short regenerative model
linearised with respect to the time delay leads to the classical process damping model,
and some other higher-order approximations can also be derived as its special cases.
21
y
fj
work-piece
Fj
j
j-1
ar
cx
x
kx
v
F rj+1
j+1
W
F
t
j+1
Fj+1
R
ky
cy
tool
Figure 2.1: 2DoF model of the flexible tool and the milling process
2.2
Surface Quality of Milled Surfaces
In this section, surface errors after stable milling operations are investigated numerically in case of helical fluted tools. A detailed two-degree-of-freedom mechanical
model is derived that includes surface regeneration. The governing delay-differential
equation is analysed by means of the semi-discretization method. The surface errors
are predicted based on the (stable) forced motion of the tool. New surface error
parameters are introduced to characterize the properties of the spatial machined surface. The errors are calculated numerically for a given machine tool and workpiece for
different axial depths of cut and spindle speeds. It is shown that both good surface
properties and large material removal rates can be achieved by means of appropriate
axial immersions in case of helical fluted tools even in those spindle speed domains
where the system is near to resonance. This phenomenon is then proved analytically by means of the Fourier transformation of the cutting force. Experiments are
performed with a large industrial milling machine using a flexible test rig with an
essential flexibility in a specified direction. The calculated and the measured forced
vibration signals in the stable domains of the stability chart are in good correlation,
which validates the simple analytical predictions for the non-trivial appropriate axial
immersion parameters. In the last subsection, the values of the non-trivial appropriate axial immersion are computed and numerically tested for a multiple degree of
freedom system.
2.2.1
Mechanical Model
In what follows, a 2 DoF flexible tool and rigid workpiece model is used (see Fig. 2.1),
because the lowest natural frequencies in the x and y modal directions of the cylindrical milling tool are considered to be the relevant ones.
22
In finishing operations, only the stable milling process is acceptable, which can
be identified by means of linear stability analysis. If the cutting force Fj of the jth
tooth is approximated as a linear function of the chip thickness hj , the tangential and
the radial cutting force components can be given as:
Fjt (t) = wKt hj (t),
Fjr (t) = wKr hj (t) = kr Fjt (t),
(2.1)
where w is the chip width, Kt and Kr are the specific cutting force coefficients, while
kr denotes their ratio, kr = Kr /Kt . These coefficients can be determined through
a set of test measurements at different machining process parameters by using force
sensors [47]. The chip thickness consists of two parts: the stationary chip thickness,
and the dynamic one. The stationary part is generated by the constant feed velocity
v, while the dynamic part is generated by the vibration of the tool. If the motion of
the tool centre is described by the general coordinate vector q = [x y]> , then the chip
thickness at the jth tooth assumes the form
τv
hj (t) = cos(η)[sin(φj (t)) cos(φj (t))]
+ q(t) − q(t − τ ) .
(2.2)
0
Here, η = arctan(2πR/p) is the helix angle (see Fig. 2.2), φj (t) is the angular position
of the jth tooth (see Fig. 2.1), the regenerative time delay τ is equal to the toothpassing period
2π
,
(2.3)
ΩN
where N is the number of the flutes and Ω is the angular velocity of the tool. The
term τ v is equal to the feed per tooth. With the help of the modal matrices M
(mass), C (damping) and K (stiffness), the governing equation is given in the form
τv
Mq̈ + Cq̇ + Kq = W(t)
+ q(t) − q(t − τ )
(2.4)
0
τ=
where the right hand side determines the resultant cutting force vector acting on
the tool. This linear non-homogeneous delay-differential equation (DDE) can be
calculated from Eq. (2.1) and Eq. (2.2). The resulting directional force coefficient
matrix W(t) is determined in the following subsections.
2.2.1.1
Straight Fluted Tool
In case of straight fluted tools, the helix angle η is zero, so it is easy to calculate
the directional force coefficient matrix W(t) from Eq. (2.1) and Eq. (2.2) (see, for
example, [38] and [48]). The result is
W(t) =
N
X
ap K t
j=1
2
g(φj (t))T(φj (t)),
23
(2.5)
z
z
W
W
dw
dz
p
p
N
j+1
x
p
N
fj+1(t,z)
j
fj(t,0)
y
R
h
local chip
y
Figure 2.2: Geometry of a helical fluted tool including the helix pitch p
where the axial depth of cut ap is equal to the chip width w. The screen function
g indicates with a value 1 if the jth tooth is in contact with the workpiece, that is,
if the cyclic angle coordinate φj (t) is between the entrance angle φin and the exit
angle φout , while it is 0 otherwise. The projection matrix T is used to project the jth
cutting force vector from the local radial-tangential coordinate system fixed to the
jth tooth to the global (x, y) coordinate system fixed to the tool centre:
−kr + kr cos(φj (t)) − sin(2φj (t)) −1 − kr sin(φj (t)) − cos(2φj (t))
T(t) =
.
1 − kr sin(φj (t)) − cos(2φj (t)) −kr − kr cos(φj (t)) + sin(2φj (t))
(2.6)
2.2.1.2
Helical Fluted Tool
If helical fluted tools are considered in the model (see Fig. 2.2) then the current
angular position of the cutting edge also depends on the axial coordinate z of the
actual tool-segment (see [39]). Thus,
φj (z, t)) = tΩ + j
2π 2πz
−
,
N
p
(2.7)
where p is the helix pitch in Fig. 2.2. Divide the tool into elementary disc segment of
height dz. In a cross section of the tool, the elementary directional force coefficient
matrix dW is calculated in the same way as in the straight fluted model Eq. (2.5)
with an infinitesimally small chip width (dw = dz/ cos(η)):
dW(t) =
N
X
Kt
j=1
2
g(φj (z, t))T(φj (z, t))dz.
(2.8)
Integrate Eq. (2.8) along the z axis from 0 to the axial immersion ap in order to
obtain the resultant directional force coefficient matrix:
24
z
z
fout fin
full period full period
W
p
N
y
p
ap N
ap
x
W
W
p
2p
W
p
2p
W
2p f2(0,z) f1(0,z) 0
fj
dW/dz
dW/df
Figure 2.3: Schematic representation of the variation of the directional force coefficient matrix W
with respect to the vertical coordinate z and also to the angular coordinate φj of the jth cutting
edge.
W(t) =
N Z
X
j=1
ap
0
Kt
g(φj (z, t))T(φj (z, t))dz.
2
(2.9)
The limit of the integration can be changed to the angular coordinate because it is
a linear function of the z coordinate as given in Eq. (2.7), where dz = −(p/(2π))dφj .
The integration by substitution leads to
W(t) =
N Z
X
φj (0,t)
φj (ap ,t)
j=1
pKt
g(φj )T(φj )dφj .
4π
(2.10)
Figure 2.3 shows that at a height z = p/N , the cutting edges’ angles cover a full
period in [0, 2π]. If the axial immersion ap is larger then p/N then the directional
force coefficient matrix in Eq. (2.10) can be separated into two parts:
W(t) = Nf ull W +
N
X
c j (t),
W
(2.11)
j=1
where Nf ull = int(ap N/p) with int denoting the integer part function. The first term
contains those sections where the integration is carried out along full periods of p/N .
This part is constant in time:
Z 2π
Z φout
pKt
pKt
W=
g(φ)T(φ)dφ =
T(φ)dφ.
(2.12)
4π
4π
0
φin
The remainder part of the integration is time dependent:
Z
φj (ap ,t)+ψ
c j (t) =
W
φj (ap ,t)
25
pKt
g(φ)T(φ)dφ,
4π
(2.13)
ap=4mm
80
ap=6mm
ap=0mm8
ap=12mm
ap=14mm
ap=18mm
ap=20mm
W11 [N/mm]
40
0
ap=10mm
80
40
0
80
40
0
ap=16mm
0
p/2
p
3p/2
2p 0
p/2
p
3p/2
2p 0
p/2
p
3p/2
2p
fj, angular position of the tooth [rad]
Figure 2.4: W11 component of the directional force coefficient matrix W for some axial immersion
values (ap ) with p/N = 10mm, ar = 5% radial immersion, R = 4mm Kt = 644 MPa and kr = 0.368.
where the "length" of this integration is
ψ = mod(φj (0, t) − φj (ap , t), 2π).
(2.14)
The structure of the directional force coefficient matrix is important here: it is a
periodic function and its period equals to the time delay τ . However, it can be time
c j (t) is identically zero. This is the case if we have
independent if the periodic part W
to integrate along full periods only, that is, when the axial immersion is just
p
,
l = 1, 2, 3, . . . .
N
This case is represented in Fig. 2.4 by straight horizontal lines.
ap = l
2.2.2
(2.15)
Stationary and Transient Motion of the Milling Tool
It can be shown that the motion of the tool can be decomposed (see, e.g., in [49]) in
the form
x(t)
ξ(t)
+
(2.16)
q(t) = qp (t) + qh (t) =
y(t)
η(t)
where qp (t) = qp (t − τ ) is the forced (chatter-free) stationary periodic motion of the
tool that corresponds to the particular solution of the linear forced DDE (2.4), and
qh (t) is its perturbation that corresponds to the diminishing transient vibration in
case of chatter-free parameters, or it describes the exponentially increasing chatter
in case of unstable machining process parameters, where the subscript h refers to the
homogeneous solution of Eq. (2.4). The forced stationary motion qp (t) satisfies
τv
Mq̈p + Cq̇p + Kqp = W(t)
,
(2.17)
0
which is an ordinary differential equation (ODE). Substituting Eq. (2.16) and Eq. (2.17)
26
into Eq. (2.4), we obtain the linear time-periodic DDE
Mq̈h + Cq̇h + (K − W(t)) qh = −W(t)qh (t − τ ),
(2.18)
which describes the stability of the stationary tool motion qp (t). Equation (2.18)
describes a delayed oscillator subject to parametric excitation. If the qh ≡ 0 trivial
solution of this equation is unstable, then chatter occurs and there is no need to
calculate the surface errors since unstable machining results unacceptable machined
surface. When the trivial solution qh ≡ 0 of Eq. (2.18) is stable, the stationary motion
of the tool centre is given by the solution of Equation (2.17). This (slightly) damped
periodically forced oscillator always has a periodic solution. By using this solution,
we can calculate the stationary motion of the teeth, which generates the machined
surface. This calculation should be carried out in the parameter domains of stable
cutting, which can be determined from Eq. (2.18) by means of the so-called semidiscretization method. This efficient method was developed, described and tested in
details in [50, 10].
2.2.3
Forced Vibration and Machined Surface
In order to calculate a particular solution analytically, we approximate the cutting
force by its truncated Fourier series
NF
2πlt
2πlt
τv ∼ X
al cos
+ bl sin
.
W(t)
=
0
τ
τ
(2.19)
l=0
During the calculations, the use of the first NF = 25 Fourier terms seems to be
satisfactory for realistic engineering applications. Similarly, the particular solution is
considered in the form
qp (t) ∼
=
NF
X
cl cos
l=0
2πlt
τ
+ dl sin
2πlt
τ
.
(2.20)
The coefficients cl and dl can be determined after the substitution of Eq. (2.19) and
Eq. (2.20) into the governing equation Eq. (2.17):
"
cl
K−
=
dl
−
2πl 2
M
τ
2πl
C
τ
2πl
C
τ
2πl 2
− τ
K
#−1
al
bl
M
l = 0, 1, . . . , NF .
(2.21)
Using the solution qp (t), the position vector rj (z, t) of the points of the jth cutting
edge at height z can be described by superposing the forced stationary vibration onto
the cylindrical rotation of the ideally rigid (non-vibratory) tool (see Fig. 2.5):
xp (t)
Rcos(φj (z, t))
rj (z, t) =
+
,
(2.22)
yp (t)
Rsin(φj (z, t))
27
Figure 2.5: Path of the helical edges r(z, t). Dark surface denotes the cutting period, bright surface
denotes the non-cutting period
where R is the radius of the tool and φj comes from Eq. (2.7). In the milling process,
each tooth creates a single surface segment sj (z, t) during one revolution. The finished
surface of the workpiece is formed by the succeeding surface segments (see Fig. 2.6).
In case of stationary cutting conditions, the surface properties can be described by a
single surface segment, which can be determined by the motion rj (z, t) of any of the
cutting edges relative to the workpiece:
vt
sj (z, ϑ(z)) = rj (z, ϑ(z)) +
,
ϑ(z) ∈ [tenter (z), texit (z)].
(2.23)
0
We have to determine the entering and exiting time instants tenter (z) and texit (z) at a
certain height z where the cutting edge just touches the boundary of the surface segment of the workpiece. They can be calculated from the intersection of two succeeding
surface segments:
sj+1 (z, tenter (z)) = sj (z, texit (z)).
(2.24)
Based on Eq. (2.23), Eq. (2.24) can be reformulated as follows:
vtenter (z)
vτ
vtexit (z)
rj (z, tenter (z)) +
+
= rj (z, texit (z)) +
.
0
0
0
(2.25)
From this, tenter (z) and texit (z) are determined by means of a simple gradient based
numerical method. Note, that the results tenter (z) and texit (z) are independent of the
edge number j.
28
workpiece
e
fac
z
des
ire
r
d su
Rtx
y
surface segment
s(z,t)
SLE(z)
texit(z=0)
MSLE
Rtz
t
tenter(z=0)
Figure 2.6: Schematic representation of the surface errors made by a helical tool
2.2.4
Surface Quality Parameters
In case of the straight fluted tool, two surface parameters are commonly in use. The
first is the total height of the profile Rt , which describes the surface roughness. It is
defined by the difference of the maximum peak height and the maximum valley depth
of the surface [51, 28, 52, 29, 53, 54]:
Rt = max(sy (t)) − min(sy (t)),
(2.26)
where the dependence on the coordinate z does not show up for straight fluted tools.
The second parameter is the surface location error SLE [49, 30, 55], which is an offset
error defined as the maximal distance between the desired surface and the machined
surface. SLE can be calculated by
SLE = max(sy (t)) − R.
(2.27)
In case of the helical fluted tool, the total height of the profile and the surface
location error are different at each height z. Consequently, the corresponding parameters Rt (z), SLE(z) are functions of the axial coordinate z (see Fig. 2.6), while the
use of scalar surface quality parameters are preferred in practice. It can be seen from
Eq. (2.22) and Eq. (2.7) that the surface is periodic along the z coordinate, too. The
variation of SLE has a wavelength p/N in the direction z which is usually larger
than the axial depth of cut – this is why this periodicity is not observable in Fig. 2.6,
either. Still, the variation of SLE is not a real roughness but rather a waviness.
Consequently, the surface roughness Rt is to be separated into two parts: the first is
the total height of the profile in feed direction x
29
Rtx = max(Rt (z)),
(2.28)
while the second part is the total height of the profile in axial direction z
Rtz = max(SLE(z)) − min(SLE(z)).
(2.29)
An overall scalar characteristic value of the surface could be the maximum surface
location error
M SLE = max(sy (z, θ(z))) − R = max(SLE(z)).
(2.30)
If we have straight fluted tools, Rtz ≡ 0 and the parameters SLE and M SLE are the
same. For helical fluted tools, the introduced three types of surface errors Rtx , Rtz
and M SLE are represented graphically in Fig. 2.6.
2.2.5
Effect of Machining Process Parameters on Profile
2.2.5.1
Real-case Mechanical Parameters
To represent the effect of the machining process parameters on the surface quality,
we choose a particular machine tool and corresponding workpiece parameters. The
machine tool is described by the same modal parameters as measured in [56]. The
authors determined the frequency response function of the machine tool by a standard
impact test procedure. The machine tool had a flexible tool with small diameter and
large overhang, so the system has one dominant mode in each direction. These modes
correspond to the first bending mode of the tool, so the two natural frequencies of
these modes are very close to each other (ωnx /(2π) = 722.1 Hz and ωny /(2π) = 722.2
Hz) due to the cylindrical symmetry. The mode coupling was also very weak. The
modal parameters like the mass (M), the damping (C) and the stiffness (K) matrices
were identified by means of the curve fitting algorithm of a commercial modal analysis
software:
mx 0
0.01986
0
M=
=
kg,
0 my
0
0.02008
Ns
cx 0
1.6031
0
,
C=
=
0 cy
0
1.1557 m
N
kx 0
408866
0
K=
=
.
0 ky
0
413445 m
During the calculations, we use tools with a standard radius of R = 4 mm. The
material of the workpiece is AlMgSi0.5 aluminium alloy, the corresponding tangential
cutting coefficient Kt and the non-dimensional force ratio kr are [57]:
Kt = 644
N
,
mm2
30
kr = 0.368.
Rt [mm]
SLE [mm]
(a)
200
100
0
w, 5
axi 4
al 3
im
me 2
rsi 1
on
[m 0 0
m]
1000
3000
]
d/s
ed [ra
e
dle sp
in
W, sp
2000
(b)
6
4
2
w, 5
axi 4
al 3
3000
im
me 2
2000
]
d
rsi 1
ra /s
1000
on
eed [
p
s
e
l
[m 0 0
ind
m]
W, sp
Figure 2.7: (a) Surface location error SLE and (b) total height of the profile Rt for down-milling
with straight fluted tool in case of stable cutting parameters (N =4, τ v=0.1 mm, p → ∞, ar = 5%)
2.2.5.2
Effect of Machining Process Parameters
The machined surface parameters should be calculated only in those machining process parameter regions where the cutting process is stable. The surface quality is
especially important in finishing operations, where typically small radial immersions
are used. Calculations using the above structural model and cutting parameters were
carried out to predict surface properties at a fixed 5% radial immersion ar for downmilling. Figure 2.7(a) shows the surface location error (SLE) and Fig. 2.7(b) shows
the total height of the profile (Rt ) for down-milling over the stability chart for a
straight fluted tool. Note, that the total height of the profile in axial direction is
always zero for the straight fluted tool.
We can see that there are large offset errors if we are in the vicinity of those
angular velocities Ω where the ratio of the natural frequency (ωn /(2π)) of the system
and the tooth-passing frequency (1/τ = ΩN/(2π)) is close to an integer. In these
cases, the tool vibration is close to resonance because one of the frequencies of the
higher harmonics of the cutting force is close to the natural frequency.
In Fig. 2.8, we calculated the surface parameters (M SLE, Rtx and Rtz ) for the
same machining conditions as in Fig. 2.7, but we used a helical fluted tool with helix
pitch of 10 mm. We can see that the helix pitch has a large influence on surface quality.
Figure 2.8 shows that both good surface properties and large material removal rate
(that is, large axial immersion at a chosen optimal spindle speed) can be achieved
by the appropriate axial immersion with a helical fluted tool. We found that the
total height of the profile in axial direction (Rtz ) is in the same order of magnitude
as the maximum surface location error (M SLE) (see Fig. 2.8(a) and (c)). A strange
pattern of the surface properties can be viewed by plotting the surface properties
above the unstable regions, too (see Fig. 2.9). Along the dashed line at p/N , when
the directional force coefficient matrix is constant (see Eq. (2.15)), there is a small
but non-zero M SLE. Here the cutting force is constant, hence, there is no forced
vibration and the total height of the profile in axial direction is zero: Rtz = 0, while
M SLE is orders of magnitude smaller than the SLE obtained with a straight fluted
tool.
Along the slanting dotted lines in Fig. 2.9, the surface properties are also very
31
(b)
(a)
6
4
2
Rtx [mm]
MSLE [mm]
32
40
20
0
w5
, ax 4
ial 3
im 2
me
rsi 1
on 0 0
[m
m]
3000
1000
2000
d/s]
ed [ra
e
dle sp
in
W, sp
w, 5
axi 4
al 3
3000
im 2
2000 d/s]
me
rsi 1
1000
ed [ra
on 0 0
e spe
l
d
[m
n
i
m]
W, sp
RRt
tz x[[m
mm]
m]
(c()b)
80
40
0
w, 5
a
xia 4
l im 3
me 2
rsi 1
on
[m 0 0
m]
1000
2000
e
dle sp
in
W, sp
3000
]
d/s
ed [ra
(b) 4
50
p
Z
0
4
w,
m]
[m
on
rsi
me
im
ial
ax
3
2
1
0
0
400
800
1200
ed [rad/s]
spe
W, spindle
50
3 p
Z
MSLE [mm]
(a)
w, axial immersion [mm]
MSLE [mm]
Figure 2.8: (a): Maximum Surface Location Error, (b): total height of the profile in feed direction
and (c): total height of the profile in axial direction for down-milling with helical tool in case of
stable cutting parameters (N = 4, τ v = 0.1 mm, p=10 mm, ar = 5%)
2
1
0
0
0
400
800
1200
W, spindle speed [rad/s]
Figure 2.9: (a): Maximum Surface Location Error above the machining process parameters, and
(b): its top view. Dark areas denote the stable region (N = 4, τ v = 0.1 mm, p = 10 mm, ar = 5%)
good, in spite of the fact that the cutting force itself is not constant there and forced
vibrations exist. This phenomenon is explained below by the Fourier transformation
of the cutting force signal.
2.2.5.3
Fourier transformation of the cutting force
The specific cutting force is distributed along the helical cutting edge and can be
calculated from the right hand side of Eq. (2.4) with Eq. (2.8) for the forced (chatterfree) periodic solution:
f (φ(z, t)) =
N
X
Kt
j=1
τv
g(φj (z, t))T(φj (z, t))
0
2
(2.31)
where, the angular position φ(z, t) of the jth edge is given by Eq. (2.7). This function
f (φ(z, t)) is periodic in time with the tooth-pass period τ due to the trigonometric
functions in the expression of T and g in Subsection 2.2.1.1. Consequently, f (φ), is
also periodic in φ with period 2π/N . By using Eq. (2.4) and Eq. (2.9), the Fourier
transformation of the resultant cutting force is defined by
Z Z ap
2πz
1 τ
f tΩ −
dz e−iωt dt.
(2.32)
Fcut (ω) =
τ 0
p
0
The axial space coordinate z is substituted by the new time variable θ according to
the formula
Ωp
z = (t − θ) ,
(2.33)
2π
and the integration by substitution leads to the convolution integral
1
Fcut (ω) = −
τ
Z
1
=−
τ
Z
τ
Z
!
t
2πa
t− Ωpp
0
τ
Z
f (θΩ)dθ e−iωt dt
(2.34)
∞
−iωt
f (θΩ)u(t − θ)dθ e
dt,
−∞
0
with
u(t) =
1 if 0 < t <
0 otherwise
2πap
pΩ
(2.35)
The Convolution Theorem states that the Fourier Transform transforms a convolution
into a multiplication [58], thus Eq. (2.34) can be written as
Z
Z ∞
1 τ
−iωθ
Fcut (ω) = −
f (θΩ)e
dθ
e−iωt u(t)dt
τ 0
−∞
Z
ωap
1 τ
i
=−
f (θΩ)e−iωθ dθ
1 − ei2π Ωp .
τ 0
ω
33
(2.36)
f (qW)
- 2p
WZ
u (t)
F(w)
Fourier transform
2p
WZ
0
-2WZ -WZ 0
U(w)
Fourier transform
q
2p w
Wp
0
-2Wp
w
t
Wp
- w
WZ
0
2WZ
Wp
w
w
2 Wp
w
w
Figure 2.10: Representation of the functions f (θΩ) and u(t) and their Fourier transforms
Since f is a periodic function, its Fourier transform is a discrete function. It can be
seen from Eq. (2.36) and also in Fig. 2.10, that the resonant Fourier component of
the cutting force Fcut (ωn ) is zero if either the first part of the product in Eq. (2.36)
is zero:
2π
k = kΩN
⇒
τ
or the second part of the product is zero:
ωn 6=
ωn ap
=k
Ωp
⇒
ap =
Ω 6=
ωn
,
kN
Ωp
k,
ωn
k = 1, 2, 3 . . . ,
k = 1, 2, 3 . . . .
(2.37)
(2.38)
The first condition Eq. (2.37) determines the values of the well-known resonant angular velocities. The straight lines ap (Ω) described by the second condition Eq. (2.38)
are the same as the slanting dotted ones in Fig. 2.9, which provide the so-called nontrivial appropriate axial immersions. It is also important to note, that the appropriate
axial immersions do not depend on the force characteristics. The dashed lines at lp/N
(l = 1, 2, 3 . . .) represent minimal surface quality parameters as already explained at
Eq. (2.15), which can also be called as trivial appropriate axial immersions.
The identification of the non-trivial appropriate axial immersions is of great practical importance because just the resonant cutting speeds are preferred to achieve
high material removal rates under stable (chatter-free) machining conditions. If nonappropriate axial immersions are used at these cutting speeds, the machined surface
quality parameters might be extremely poor.
In the next Subsection, our goal is to verify the non-trivial appropriate axial
immersions given in Eq. (2.38) by measurements via constructing an experimental
superchart, which is the combination of the traditional stability chart and the forced
vibration amplitudes which are directly related to some surface quality parameters
(see [30, 59]).
34
milling tool
Steel C45 workpiece
accelerometer
workpiece
holder z
y
vibration
acquisition
system
flexible
parts
Figure 2.11: Measurement setup: test rig flexible in y direction.
p
ωn = k/m
600.6 rad/s
fn = ωn /(2π)
95.6 Hz
ξ = c/(2mωn ) m (kg) c (Ns/m)
0.061%
136
996.63
k (N/µm)
49.069
Table 2.1: Measured modal parameters.
2.2.6
Experimental setup
The measurements are performed with a large industrial milling machine (Soraluce
SV-6000) on which a test rig is mounted with an essential flexibility in y direction (see
Fig. 2.11). The experimental modal analysis performed by impact tests shows that
the system has only one dominant vibration mode, indeed. The verified modal parameters like the natural angular frequency ωn and the damping ratio ξ are presented
in Table 2.1. With the help of the measured mass m of the workpiece holder and
the workpiece, all the mechanical parameters were identified including the stiffness k
in y direction and also the corresponding damping factor c. The tool and the tool
holder together with the whole machine tool structure proved to be rigid relative to
the workpiece holder.
2.2.7
Results and Experimental Validation
In the test case and in the experiments, a two fluted helical cutting tool (N = 2) was
used with 30 mm diameter, 42◦ helix angle and p=104.67 mm helix pitch (Garant
191100 STG 104.67 NC HSS-Co8). The machined material was Steel C45 and its
corresponding cutting coefficients (see Table 2.2) were measured based on a method
explained in details in [32].
During the measurements, constant feed per tooth fz = τ v = 0.15 mm/tooth was
used with radial immersion ar = 0.75 mm (ar /D = 2.5%) in down-milling operation.
Accordingly, the enter and exit angles are φenter = 161.8◦ and φexit = 180◦ , respeck1t (N/mm2 ) k1r (N/mm2 ) k0t (N/mm) k0t (N/mm)
1697
476
58
37
Table 2.2: Specific cutting coefficients.
35
wn
p
2W
w=
w=p/N=52.335 [mm]
W
w=
p wn
800
non-trivial
appropriate
axial
immersions
50
40
600
30
400
20
20
0
20
10
-50
0
50
y(t) [mm/s]
ap -axial immersion [mm]
trivial appropriate axial immersion
200
y(t) [mm]
0
0
50
100
925
150
1387
200
W - spindle speed
250
2500
300
2775 [rpm]
350 [rad/s]
0
max(y(t))-min(y(t)) - amplitude of the forced vibration [mm]
1000
wn
w=
3W
p
60
Figure 2.12: Superchart of milling. Calculated dark gray regions (lobes) refer to chatter (unstable
stationary milling). Calculated amplitudes of the forced vibrations are presented by gray scale in
the stable domains (stable pockets between the lobes). Small panels present measured vibrations in
the phase plane (y, ẏ) at the measurement parameters denoted by tiny circles. A trivial appropriate
axial immersion and the non-trivial ones are denoted by a horizontal dotted line and slanting dashed
lines, respectively.
tively. Our goal was to verify the non-trivial appropriate axial immersion values (see.
Eq. (2.38)). For this reason, the experiments were carried out closed to the first three
resonant spindle speeds n = Ω60/(2π) =2775, 1387, 925 rpm and at a non-resonant
spindle speed n=2500 rpm, for a set of five axial immersions ap =43.6, 34.9, 26.1, 17.4,
8.7 mm. These cutting parameters are denoted by tiny circles in Fig. 2.12.
In this figure, the trivial appropriate axial immersion given by Eq. (2.15) and
the non-trivial ones given by Eq. (2.38) are denoted by a horizontal dotted line and
slanting dashed lines, respectively. During the experiments, the acceleration of the
workpiece was measured by piezo accelerometers then the time signals were integrated
twice numerically and filtered by high-pass filter at the cut off frequency 5 Hz. The
two dimensional projections of the (infinite dimensional) phase space to the plane
(y, ẏ) of the measured signals are plotted in Fig. 2.13 after the transient motions
died out. Note the non-uniform scales of the panels. The same panels show the
corresponding calculated periodic forced vibrations, too, where the static deformation
component is subtracted to be able to compare the results with the measured highpass filtered signals. The same phase plots of the measured vibrations are inserted
to the superchart of Fig. 2.12 where a uniform scale is used. This way, Fig. 2.12
presents a comparison of the measured and theoretically predicted forced vibrations
in the stable cutting parameter domains. Another representation of the measured
data can be found in Fig. 2.14, where the measured and calculated peak-to-peak
amplitudes maxt (y(t))−mint (y(t)) of the periodic forced vibrations are plotted for the
selected spindle speeds denoted by continuous vertical thin lines in Fig. 2.12. All these
36
representations in Fig. 2.12, Fig. 2.13 and Fig. 2.14 show that the measured vibration
amplitudes are significantly smaller in the vicinity of the non-trivial appropriate axial
immersions determined in Eq. (2.38), which provides an experimental validation of
the theoretical predictions.
We measured smaller vibration amplitudes by setting the non-trivial appropriate
axial immersion values in case of resonant spindle speeds, than the vibration amplitudes at any arbitrary axial immersion in case of a non-resonant spindle speed
(compare, for example, the second and the third vertical lines in Fig. 2.12). Note,
that even smaller forced vibrations can be achieved at non-resonant spindle speeds
and at appropriate axial immersions, although the advantage of high efficiency between the instability lobes cannot be exploited in these cases. The main advantage
of the application of non-trivial appropriate axial immersions is that the actual axial immersion can be set to a much smaller value then the helix pitch, which is the
first trivial appropriate axial immersion. The trivial appropriate axial immersions
in Eq. (2.15) often cannot be applied due to certain limits of the machining process
parameters, or simply due to the length of the milling tool as it was the case in our
experiment series, too.
2.2.8
Multiple Degree of Freedom System
In this section, the equations 2.38 and 2.15 are analysed for a given system with
multiple degrees of freedom which is flexible in y direction only (see Fig. 2.15). In
this case, the dynamical properties of the tool can be described by the frequency
response function Hy of the tool tip
Hy (ω) =
F(y)(ω)
,
F(Fy )(ω)
(2.39)
where F refers to the Fourier transform. This can be measured by means of a simple
impact test procedure. After rearrangeing Eq. (2.39) and using Fourier F and inverse
Fourier F −1 transformations, the periodic tool motion can be given as
yp (t) = F −1 Hy (ω)F(Fy )(ω) (t),
(2.40)
where the cutting force function Fy (t) along the flexible direction y can be calculated
from Eq. (2.10) – Eq. (2.17). In practice, the Fourier F and inverse Fourier F −1
transformations are numerically determined by the Fast Fourier Transformation (fft)
and the inverse Fast Fourier Transformation (ifft) routines in Matlab. The following steps of the surface error calculations can be repeated following the procedure
described in Eq. (2.22) – Eq. (2.30).
Test Case
The selected 3 DoF mechanical model is represented by its modal parameters [8]
which are given in Table 2.3 based on the measurements documented in [60].The
corresponding tooltip FRF is represented in Fig. 2.16. The cutting coefficients and
37
38
W - spindle speed
0
0
0
-10
0
10
-0.2
-200
0
200
-0.01
-20
0
20
-0.2
-500
0.2
0.02
0.5
0
0
0
-0.2
-200
0
200
-0.02
-50
0
50
0.01
0.02
0.5
0
0
0
0
-0.1
-200
0
200
0.01
0
-10
0
10 20
-0.01
-20
0
20
-0.02
-50
0
50
-0.5
-500
0.2
0.02
0.5
0
0
0
-0.2
-200
0
200
-0.02
-50
0
50
-0.5
-500
0.1
0.2
0.01
0.2
0
0
0
0
-0.1
-200
0
200
-0.2
-200
y(t) [mm]
0
200
-0.01
-20
y(t) [mm]
0
20
-0.2
-500
y(t) [mm]
500
theory
-0.5
-500
0.1
0
experiment
0
500
0
500
0
500
0
500
26.1 [mm]
y(t) [mm/s]
y(t) [mm/s]
200
43.6 [mm]
0
34.89 [mm]
0
ap - axial immersion
0.2
17.4 [mm]
y(t) [mm/s]
0.01
0
2775 [rpm]
290.60 [rad/s]
0.2
-0.010
-20
y(t) [mm/s]
2500 [rpm]
261.80 [rad/s]
0.1
-0.1
-200
y(t) [mm/s]
1387 [rpm]
145.25 [rad/s]
8.7 [mm]
925 [rpm]
96.86 [rad/s]
y(t) [mm]
Figure 2.13: Measured (gray continuous line) and calculated (dashed black line) phase plots of the
forced periodic vibration of the workpiece. Note changes of scales!
max(y(t))-min(y(t)) - amplitude
of the forced vibration [mm]
100
W=2775 [rpm]
50
0
40
W=2500 [rpm]
20
0
40
W=1387 [rpm]
20
0
40
W=925 [rpm]
20
0
0
10
20
30
40
w - axial immersion [mm]
50
60
Figure 2.14: Amplitudes of the measured and the calculated forced vibrations denoted by gray
continuous lines and circles, and dashed black lines, respectively.
39
v
work-piece
y1
y
x
0
W
tool
y2
0
y3
0
Figure 2.15: Representation of a system with multiple degrees of freedom which is flexible in y
direction only.
wn 2
Hy(w) - frequency response function [mm/N]
0.9
0.8
0.7
0.6
wn 1
0.5
0.4
0.3
wn 3
0.2
0.1
0
0
200
400
600
800 1000 1200 1400 1600 1800
w - frequency [rad/s]
Figure 2.16: Frequency response function of the tool-tip.
mi
(kg)
1 136
2 27.2
3 204
i
ci
(Ns/m)
3986
1275
4584
ki
(N/µm)
49.07
25.12
389.3
ωni
fni
(rad/s) (Hz)
600.7
95.6
961.1 153.0
1381.5 219.9
ξi
(1)
2.44
2.44
0.813
Table 2.3: Selected modal parameters.
the technological parameters of the down milling operation are given in Table 2.4 and
Table 2.5, respectively.
It can be seen in the presented superchart of Fig. 2.17, that the large material
removal rate (large axial immersion and large spindle speed) can be achieved between
the stability lobes where the spindle speeds are close to resonance. In Fig. 2.17
the resonant spindle speeds (ωn1 /N, ωn1 /2N,ωn1 /3N,. . .;ωn2 /N, ωn2 /2N,. . .;ωn3 /N,
ωn3 /2N,. . .) are denoted by vertical dashed lines. Based on Eq. (2.15), the trivial
appropriate axial immersion is plotted by a horizontal continuous line while based on
Eq. (2.38), the non-trivial appropriate axial immersion lines are plotted by dashed
slanting lines for all the natural angular frequencies (ωn1 , ωn2 , ωn3 ).
It can be seen, that in case of resonant spindle speeds, small amplitude vibration
occurs only if the corresponding appropriate axial immersions are used (denoted by
white small circles in Fig. 2.17). Thus, the analytically derived equations 2.15 and
2.38 can effectively be used even for systems with multiple DoF. Note, that the stable
parameter domains reach the trivial appropriate axial immersion (p/N =50 (mm))
k1t (N/mm2 ) k1r (N/mm2 ) k0t (N/mm) k0t (N/mm)
1889.1
775.5
63.1
78.2
Table 2.4: Selected specific cutting coefficients [60].
40
N
(1)
2
D
(mm)
60
p
(mm)
100
fz
(mm/tooth)
0.15
ar
(mm)
3
Table 2.5: Applied technological parameters of a down milling process.
700
ap - axial immersion [mm]
p
N
non-trivial appropriate
axial immersions
600
trivial axial
appropriate
immersion
30
500
400
UNSTABLE
20
UNSTABLE
300
UNSTABLE
200
10
100
STABLE
0
100
. . . w2
3N
200
w2 w1
2N N
STABLE
400
500
w3
w2
2N
N
W - spindle speed [rad/s]
600
STABLE
w3
N
800
max(y1(t))-min(y1(t)) - peak-to-peak amplitude [mm]
50
0
Figure 2.17: Superchart of milling for 3 DoF milling system with an N = 2 fluted tool. Calculated
dark gray regions (lobes) refer to chatter (unstable stationary milling). Calculated amplitudes of the
forced vibrations are presented by gray scale (which have physical meaning in the stable domains
only). A trivial appropriate axial immersion and the non-trivial ones are denoted by a horizontal
continuous line and slanting dashed lines, respectively. The vertical lines denote the resonant spindle
speeds.
only at large spindle speeds, however, these spindle speeds cannot be reached due to
some limitations like maximal cutting speed, maximal power consumption, maximal
spindle speed, size of the workpiece or cutter working length, etc.
New Results
Thesis 2
Peripheral milling tools are considered which can vibrate with two modes normal
to the tool axial direction. It is shown that the helical fluted tools produce complex
spatial machined surface even along the axial direction. Three specific parameters are
introduced to characterize these spatial machined surfaces created by helical fluted
milling tools, which are sufficient to optimally describe the surface and easy-to-use
in engineering practice. The first parameter is the maximal surface location error,
which is a generalization of the classical surface location error defined for straight
41
fluted tools. The second parameter is the total height of the profile in feed direction
that can also be considered as surface roughness. This parameter is shown to be only
slightly affected by the helix angle of the tool. The third introduced parameter is the
total height of the profile in axial direction that can also be considered as a kind of
surface waviness. This is the most relevant parameter in case of milling with helical
fluted tools, which has no classical counterpart.
Related publications:
[61, 62, 63, 64]
Thesis 3
Classical stability charts of milling processes predict that the highest material removal
rates can be achieved when one of the harmonics of the periodic cutting force is close
to one of the natural frequencies ωn of the machine tool structure. The corresponding near-resonant spindle speeds are between the subsequent instability lobes of the
stability charts. Cutting with one of these resonant vibrations does not result selfexcited vibration (chatter) but it still causes large surface waviness in case of helical
tools, which cannot be compensated.
It is proven analytically and verified by measurements for helical milling tools
that the resonant vibrations caused by the periodic cutting force can be minimized
for certain cutting parameters that can be arranged into two classes. The first class
is formed by the trivial appropriate axial immersions ap = jp/N that are equal to the
helix pitch p over the number of teeth N and its integer multipliers j. In this case, the
cutting force is constant, no forced vibration occurs, which is a direct consequence
of the application of helical tools. It is shown that non-trivial appropriate axial
immersions exist, defined by
pΩ
j
j = 1, 2, 3, . . . ,
ωn
which depend linearly on the spindle speed Ω. At these parameter values, the cutting force varies periodically but the harmonic component close to resonance is still
eliminated.
ap =
The above simple formula of the non-trivial appropriate axial immersions can also
be used in case of multi DoF systems if the resonant spindle speeds are isolated.
Related publications: [61, 65, 62, 63, 64, 66, 67, 68]
42
2.3
Surface Reconstruction for Combined FRF and
FEM model
In this section, a method is presented to reconstruct surface topographies of peripheral
milled surfaces based on measured cutting forces in case of flexible finger-like milling
tools.
The dynamic behaviour of the system is described by means of coupling two
different methods. First, the spindle of the machine tool is analysed based on the
frequency response measurements [69, 8]. Second, the dynamic behaviour of a general
tool is modelled using Timoshenko beam elements with axial, torsional, and lateral
dynamics [70, 71]. The tool is fixed in the tool holder which is attached to the spindle,
while the workpiece is assumed to be rigid and rigidly fixed. The overall dynamics of
this system is constructed from the measurement of the spindle and from the Finite
Element Model of the tool.
Then, the measured resultant cutting force is applied to the governing equation of
the coupled system. This way, there is no need to create an advanced/sophisticated
cutting force model which should contain all the known and unknown parameters
of phenomena like non-linear force characteristics, short regenerative effect, process damping, state dependent delay and/or runout. The accuracy of the presented
method depends only on the accuracy of the measured force signals and the coupled
dynamic model of the tool and the spindle (machine tool). The calculated motions
based on the measured forces describe all the periodic forced vibrations, the chatter
vibrations, and even the transient vibrations.
The proposed method is able to reconstruct the accurate shape and roughness of
machined surfaces. It is verified by comparing the reconstructed and the measured
surface topographies. The results demonstrate that the method is able to reconstruct
the surface topography of the machined workpiece from measured resultant cutting
forces and it can be used also for the online monitoring of milling processes.
2.3.1
Development of a Dynamic Tool Model
In order to analyse the forced motion of the coupled tool and machine tool structure
described by the general coordinate vector r, the governing equation of the obtained
system can be determined in the following form (see Eq. (2.4)):
Mr̈(t) + Cṙ(t) + Kr(t) = F(t).
(2.41)
where F(t) is the generalized force originated in the resultant cutting force. Since, the
dynamic model is constructed by means of coupling the above described two parts, the
general coordinate vector and also the mass, the damping and the stiffness matrices
are the results of the concatenating of the two parts as described in the subsequent
Sections.
43
2.3.1.1
Dynamic Model of the Spindle Part
The FRF of the spindle is specific for a machine tool and has to be measured only
once. The modal parameters (natural frequency ωn,l , damping ratio ξl , and the modal
constant Bl ) for the lth mode are determined by fitting rational fraction polynomials
for the FRF, which is a widely used procedure in the industry [69, 8]. It was assumed
that the bending stiffness of the spindle is practically much larger than the bending
stiffness of the finger-like tool, which means that the bending compliance of the spindle
was neglected. Also, we assume that the shaft is torsionally and axially rigid. Then
the dynamic properties of the spindle can be characterized in its modal space given
by the diagonal matrices.
2
ωn,l
2ξl ωn,l
1
s
s
s
, CFRF =
, KFRF =
(2.42)
MFRF =
Bl
Bl
Bl
The coupling of the two systems is defined along the Cartesian coordinates of the
contact points between tool and spindle. These are given by means of the sum of the
modal coordinates qls (t), therefore the general coordinates rs (t) of the spindle part
are redefined as follows:
rls (t) =qls (t)
s
(t) =
rN
m
Nm
X
l = 1, 2, . . . , Nm − 1,
(2.43)
qls (t),
l=1
where Nm is the number of the fitted eigenmodes of the spindle. The inverse transformation of Eq. (2.43) is given by the transformation matrix T̂ as follows:
q1s (t)
1
0 ... 0
r1s (t)
q s (t)
0
s
1 . . . 0
2
r2 (t)
s
s
q (t) = .. = T̂r (t) = ..
(2.44)
.
.. . .
.
.
. ..
.
s
s
qN
(t)
−1 −1 . . . 1 rN
(t)
m
m
The mass, damping and stiffness matrices are transformed into the space of these new
coordinates with the same transformation matrices T̂ :
Msgen = T̂> MsFRF T̂,
Csgen = T̂> CsFRF T̂,
Ksgen = T̂> KsFRF T̂,
(2.45)
This method is used both in x and y directions, separately. If we assume that there
is no cross coupling between the two directions (for example, there are no gyroscopic
effects), then the governing equation of the spindle system can be described by block
diagonal matrices created from the direct sum of the matrices (Eq. (2.45)) in x and
y directions.
44
s,x
s,x s,x
s,x s,x
0
0
Kgen
Cgen
r (t)
ṙ (t)
r̈ (t)
= 0.
+
+
s,y
s,y
rs,y (t)
0
K
0
Cs,y
ṙ
(t)
r̈s,y (t)
gen
gen
(2.46)
Note, that the general force is zero in the right-hand-side due to the spindle having
no external excitation.
0
Ms,x
gen
0
Ms,y
gen
2.3.1.2
Dynamic Model of the Tool Part
The mass, damping and stiffness matrices of the tool could be generated in the same
way as the matrices of the spindle. However, if the mode shapes of the tool have to be
modelled [8, 7], then the measurements have to be made in multiple points along the
tool length and more advanced methods are needed for the evaluation of the system
matrices. Furthermore, several different endmills are used in the industry and the
measurements of each tool step-by-step would take unacceptably long time.
Consequently, the dynamic behaviour of the tool is calculated by finite element
method (FEM). This method can sufficiently predict the dynamic properties of the
tool (like natural frequencies and damping ratios) due to the fact, that the endmill tool
is usually a monolithic part, and does not contain any joints or contact surfaces that
could cause inaccuracies. Note, that the mode shapes of the tool are also presented
by this method. For the finite element model of the tool, a parametric description of
the endmill geometry was developed. The goal is to get a detailed description of the
endmill geometry by means of a minimum of describing parameters. Therefore, the
tool geometry is described by commonly used parameters [72], like the tool radius
R, number of edges N and the helix angle η. In order to calculate the variation of
stiffness along the tool axis, the cross-section zone was modelled for several slices of
the endmill. The boundary curve b of the cross-section of the endmill at height z can
be described by the simplified equation:
R(1 − δ cos(Φ)) sin(θ + φ(z))
b(z) =
θ ∈ [0, 2π],
R(1 − δ cos(Φ)) cos(θ + φ(z))
N
z
Φ = mod(θ, 2π/N );
φ(z) = tan(η);
4
R
(2.47)
where θ is the polar coordinate of the tool, Φ is the tooth pitch angle and φ is the
cutting edge position angle of the first tooth (see Eq. (2.7)). The tool cross-section
is represented by the nominal radius of the tool R, the number of edges N , the chip
space ratio δ (see Fig. 2.18) and the helix angle η. The related sampled boundary
curve b̂ is a polygon defined by the x and y coordinates of the points that are collected
in the matrix b̂ ∈ R2×Nsamp . The generated cross-section simplifies the complex shape
of a real tool cross-section by a minimum of describing parameters (see Fig. 2.18).
The variation of the cross-section geometry along the tool axis z can be described by
the variation of the chip space ratio δ(z). In case of a helical tool, the cutting edge
position angle φ(z) also depends on the axial coordinate. This way, the parametric
45
Figure 2.18: Real cross-section of an endmill a) and simplified parametric cross-section b).
model is able to generate the cross-section geometry of the tool along the cutting edge
length, the tool shaft and even for the tool chuck, where the radius R also depends
on the axial coordinate z.
In order to calculate the stiffness of the tool, the areas and moments of inertia
of the cross-sections are calculated using the trapezoidal method for the polygonal
sampled cross-sections:
l
1X
(yk+1 xk − yk xk+1 ) ,
A(z) =
2 k=1
(2.48)
l
Ixx
1 X 2
=
yk+1 + (yk + yk+1 )yk (yk+1 xk − yk xk+1 ) ,
12 k=1
Iyy
1 X 2
=
xk+1 + (xk + xk+1 )xk (xk+1 yk − xk yk+1 ) ,
12 k=1
l
(2.49)
where xk and yk are the coordinates of the sampled polygon. Note, that Ixy = Iyx = 0
if the tool has at least 2 cutting edges, while Ixx = Iyy for at least 3 cutting edges.
The discretization of the tool along the z axis is represented in Fig. 2.19. The number
of beam elements Nb depends on the desired number of the modelled eigenmodes of
the tool and the required accuracy of the surface generation.
The dominant bending vibration modes of the endmill influence primarily the machined surface topography. Hence, beam theory was chosen for the endmill modelling.
Based on the Timoshenko’s beam theory [73, 74, 75] the bending of the jth node of
the tool’s beam element can be described by four general coordinates referring to 4
degrees of freedom:
rbj = [uj , vj , θjx , θjy ] ,
(2.50)
where the superscript "b" refers to "beam", u and v are the lateral displacements
46
spindle
beam model of the tool
qx,j
vj
x
z
uj
qy,j
y
1
2
3
4 ...
...
Nb-2
Nb-1
Nb
j
Figure 2.19: Representation of the beam model with changing cross-sections and the connection to
the spindle.
in x and y directions, respectively, and θx and θy are the rotations around the x
and y axes, respectively. The mass and stiffness matrices of one beam segment are
determined based on [70] by using the area and the moments of inertia of the current
cross-section (see Eq. (2.48) and Eq. (2.49)), and the material properties of the tool
(density, Poisson ratio and Young’s modulus). The global mass and stiffness matrices
of the beam model of the tool are built up based on the neighbouring elements [73].
The governing equation of the tool dynamics is given by:
Mbglob r̈b (t) + Cbglob ṙb (t) + Kbglob rb (t) = F(t).
(2.51)
where the global proportional damping matrix Cbglob is calculated from the stiffness
matrix: Cbglob = κKbglob . In the computation, the frequency-independent damping
coefficient κ was chosen to create the typical damping ratio ξ1 = 1% for the first
mode. Its value can be calculated by:
2ξ1
,
(2.52)
ωn1
where the first angular natural frequency ωn1 of the tool is computed from the wellknown characteristic polynomial:
κ=
det −ωn2 Mbglob + Kbglob = 0.
2.3.1.3
(2.53)
Global System Matrices
As shown in Fig. 2.19, a flexible connection is modelled between the spindle and
the tool. In Subsection 2.3.1.1, only the general coordinates of the spindle in x and
y directions are defined. Thus, coupling of the matrices can be done between the
s,x
s,y
variables rN
and uNb , and rN
and vNb by considering linear spring (kxce , kyce ) and
m
m
ce
damping (cce
x , cy ) elements between the overlapping coordinates. The stiffness and
47
u1
v1
x
1
y
1
r 2b
b
K glob
r 3b
k yce
-k yce
...
r s,y
uNb
v Nb
x
Nb
y
Nb
r 1s,x
r 2s,x
r 3s,x
...
spindle-x
r s,x
spindle-y
b
r Nb
s,x
r Nm
r 1s,y
r 2s,y
r 3s,y
k xce
-k xce
ce
k qx
ce
k qy
s,x
K gen
s,y
K gen
...
beam model
r 1b
s,x
r Nm
Figure 2.20: Pattern of the coupled symmetric stiffness matrix (white regions denote zero elements).
s,x
uNb ]> coordinates are given
damping matrices of the coupling element for the [rN
m
by:
ce
ce
kx
−kxce
cx −cce
coup
coup
x
Kelem =
,
Celem =
.
(2.54)
cce
−kxce kxce
−cce
x
x
The coupling matrices along the y direction have a similar form. Due to the infinite
bending stiffness of the spindle, the bending coupling can be modelled by torsional
ce
spring kθx
and torsional damping cce
θx elements between the angular coordinate θNb x
and the "ground" (same applies for θNb y ). The final governing equation has the same
form as Eq. (2.41) with general coordinates r = [rb1 rb2 . . . rbNb rs,x rs,y ]> . Figure 2.20
represents the pattern of the coupled stiffness matrix. The specific cutting forces in
x and y directions are applied at the coordinates uj and vj , j = 1, 2, . . . , Nact , where
Nact is the number of "active" elements being in contact with the workpiece. These
coordinates belong to the active cutting length, where the zj coordinate of the jth
cross-section is smaller than the depth of cut ap .
2.3.1.4
Forced Motion
The main goal of the developed method is to apply measured cutting forces on the
endmill model and to estimate the surface quality based on the forced motion. During
the cutting process, the resultant cutting forces in the x, y and z directions are
provided by a force dynamometer. The specific cutting forces acting on the cutting
edges may have non-trivial variation along the axial coordinate z, but for the forced
motion it is supposed that the measured cutting force is distributed uniformly along
the active engagement zone. This assumption is illustrated in Fig. 2.21.
The generation of the forced motion based on the full matrices in Eq. (2.41) would
48
force
measurement
F
PC
force
acquisition
system
Fxyz
forced motion simulation
surfcae modelling
theoretical force
modelled force
workpiece
Fx
W
ae
Fy
ap
1 2 ...
vf
j
Nact
Figure 2.21: The measurement setup and the modelled uniform force distribution.
take too much computational effort due to the large sizes of the matrices. Therefore,
the system is reduced by means of the modal transformation [8, 7]. The system is
transformed into the modal space from the space of the general coordinates (Tq = r)
as follows:
T> MTq̈(t) + T> CTq̇(t) + T> KTq(t) = T> F(t),
(2.55)
where the transformation matrix T contains the mode shapes T = [A1 , A2 , . . . ,
ANb +NFRFx +NFRFy ] calculated by means of the eigenvalue problem:
−ωn2 Mbglob + Kbglob Al = 0.
(2.56)
The size of the system can be reduced by selecting the relevant modes only. In our
case, only the first Nmod = 50 modes were modelled, and the truncated transformation
matrix T = [A1 , A2 , . . . , ANmod ] was used. Equation (2.55) could be solved in the time
domain, however, if periodic cutting processes are considered, the frequency domain
solution is more appropriate (see: Subsection 2.2.8). The periodic forced motions of
the selected modes are given by means of the inverse Fourier transform:
(2.57)
qp (t) = F −1 (T> (−ω 2 M + iωC + K)T)−1 F T> F (ω) (t),
For the surface generation, only the general coordinates of the active segments
(j = 1, 2, . . . , Nact , see Fig. 2.21) of the tool are needed as explained in the subsequent
section.
2.3.2
Generation of Surface Topography in Peripheral Milling
To be able to reconstruct the surface topography for peripheral milling processes
with a flexible finger-like tool, the material removal model (see Subsec. 2.2.3) must
be evaluated for each active segments of the tool. Based on Eq. (2.22) and Eq. (2.23),
the edge path of the jth tooth at the axial coordinate zk (k = 1, 2, . . .) relative to the
49
Figure 2.22: Principle of the material removal Boolean calculation and the blank workpiece denoted
by w0 .
workpiece is given by
ukp (t)
Rcos(φj (zk , t))
v t
rj (zk , t) =
+
+ f ,
vkp (t)
Rsin(φj (zk , t))
0
(2.58)
where the lateral displacements of the forced motions ukp and vkp are determined
from Eq. (2.57), φj (zk , t) is the angular position the jth edge (from Eq. (2.7)) at
axial coordinate zk and vf t describes the feed motion. This equation describes the
relevant motion of the cutting edges along the tool path in the x − y plane. In the
same way as the endmill has been described along the tool axis, the workpiece can be
described as a number of polygons over the height of the workpiece. These polygons
give the geometry of the blank part wk0 in the kth level (see Fig. 2.22). It is now
possible to describe the workpiece contour at one time step t as the subtraction of
the subsequent tool paths and the previous workpiece contour. This can be done by
Boolean operations like intersections and subtractions between the sampled polygons
of the subsequent tooth paths (see Subsec. 1.4). The resulting geometry wk in the
kth level contains the contour of the resulting workpiece (see Fig. 2.22).
The shape of the machined workpiece, including the surface roughness of the
machined contour, can be obtained by means of a level curve representation of the
workpiece [76]. The material removal procedure is done separately for each level curve
of the workpiece. The generated level curves of the machined workpiece incorporate
the geometric profile of the machined surface.
Such profile contains surface roughness caused by the feed movement of the endmill (Fig. 2.22) and it also incorporates the surface shape due to the occurring tool
vibrations. The workpiece model contains a 2.5-dimensional representation of the
workpiece volume (see Fig. 2.23). Thereby, it presents discrete information along the
axial direction z.
50
Figure 2.23: Generation of the 2.5-dimensional workpiece model (a) and meshed surface topography
of the machined workpiece (b).
The finer the resolution of the discretization of the tool and the workpiece geometry is, the better the information on the vertical shape of the machined workpiece
will be. To be able to examine these surface characteristics, the machined workpiece
shape was meshed and interpolated to receive a smooth mesh. Then the meshed
surface can be compared to the requested surface shape and thereby, quality indicators like surface roughness, maximal surface loaction error (MSLE) and further shape
deviation parameters can be inspected.
2.3.3
Experimental Verification
In order to evaluate the developed method, the reconstruction of machined surfaces
is carried out for measured cutting forces and the reconstructed surface topographies
are compared to the measured ones as prescribed in the previous section.
2.3.3.1
Description of the Experimental Setup
To implement the dynamic behaviour of the machine tool into the global system
matrix (see Subsec. 2.3.1), the dynamic response of the machine tool spindle was
analysed while the structure was excited. This was done by a laser vibrometer using
the impact hammer method. The response was captured in two orthogonal directions
x and y (see Fig. 2.19). Thereby, it was found, that the frequency response function of
the spindle system is almost cylindrically symmetric, so the cross effects are neglected
in this model.
In the finite element beam model of the tool, the geometric data and the material
data of the endmill are given in Table 2.6. During the computation, Nb = 250 FEM
nodes were used with the system coupling parameters in [70] (see Table 2.7). The
51
Diameter
Number of edges
Helix angle
Helix angle
Young’s modulus
10 mm
3
45◦
10.47 mm
200 GPa
Cutting length
50 mm
Total length
115 mm
Chip space ratio
50%
Density
7.8 kg/dm3
Poisson number
0.33
Table 2.6: Milling tool geometry and material properties.
Stiffness
Damping
Latheral
Bending (x, y)
ce
ce
300 N/µm kθx , kθy 1.6 × 106 Nm/rad
ce
182 Ns/m cce
1240 Nms/rad
θx , cθy
kxce , kyce
ce
cce
x , cy
Table 2.7: System coupling parameters [70].
first relevant natural frequencies and mode shapes are plotted in Fig. 2.24 by using a
meshed level curve representation of the tool.
For the examination of tool dynamics in milling processes, an experimental setup
was designed. The setup is illustrated in Fig. 2.21. During a peripheral milling
process, the cutting forces were measured by a three component force dynamometer
(Kistler 9257 B). The experimental parameters were chosen for the variation of spindle
speed and constant engagement conditions Table 2.8.
By varying the cutting velocity vc , the tooth engagement frequency fte was changed
in ten steps. Thereby, the excitation of the machine tool was influenced. Ten different
milling processes were performed by measuring the cutting forces in the workpiece
coordinates. Afterwards, the surface topographies of the machined surfaces were captured by a confocal laser microscope. The experimental results are compared to the
reconstructed surface topographies in the next section.
wn1= 1467 [Hz]
wn2= 1467 [Hz]
wn10= 2335 [Hz]
wn26= 6683 [Hz]
wn47= 21490 [Hz]
Figure 2.24: Representation of the magnified vibration modes of the coupled tool-spindle system.
52
Nr.
ap (mm)
ar (mm)
fz (mm)
vc (m/min)
fte (Hz)
n (rpm)
1
40
63.7
1273
2
60
95.5
1910
3
80
127.3
2546
4
5
100
159.2
3183
6
16
1
0.05
120
140
191.0 222.8
3820 4456
7
8
9
10
160
254.6
5093
180
286.5
5730
200
318.3
6366
220
350.1
7003
Table 2.8: Process parameters for milling.
Figure 2.25: Measured surface topography (a) and reconstructed surface based on the measured
cutting forces (b) and reconstructed surface based on the theoretical cutting forces (c) for parameter
set 1 in Table 2.8
2.3.3.2
Reconstruction of Surface Topographies
In the first experimental run, moderate process conditions were chosen to be able to
validate the influence of tool deflection and surface generation mechanisms. The first
milling process was performed using the parameter set 1 in Table 2.8. The milling
process was found to be stable.
The cutting forces from the milling process were recorded and applied to the
dynamic tool model. The machined surface was reconstructed by applying the tool
motion to the material removal model. The measured (a) and the reconstructed
surface topography (b) are compared in Fig. 2.25.
As it can be seen from the measured surface, that the topography contains a
characteristic shape in the axial direction. Surface roughness can also be obtained
53
from the measured surface, but it is affected by measurement noise. The surface
roughness of the reconstructed surface (see Fig. 2.25b) includes some additional noise
from the measured cutting forces. The agreement of the characteristic shape of the
machined surface can be obtained from the profiles of the topographies (a) and (b).
Referring to this validation for moderate process conditions, the developed method
for the reconstruction of surface topographies shows very good agreement between
measured and reconstructed surface shapes.
In order to verify the developed model even further, a third test case was performed by using theoretically calculated cutting forces (see Sec. 2.2.1). The simulated
cutting forces do not include any process vibration or measurement noise, so its different effects on surface generation can be separated. Computed cutting forces were
applied to the tool model and the resulting deflections were used for the simulation
of the theoretical surface topography (see Fig. 2.25c). It can be seen, that surface
roughness is very small due to uniform engagement conditions and the absence of
measurement noise. The nominal profile shape of the measured surface topography
could not be predicted well, because the cutting force model was implemented in open
loop, thus, the material removal and cutting force generation is not influenced by the
tool deflection. In this case, the developed method for the reconstruction of surface
topographies from measured cutting forces delivers qualitative information for the
validation in milling processes.
Furthermore, the developed method is able to reconstruct transient changes in
the cutting conditions. Figure 2.26 shows a comparison between measured and reconstructed topography for an emerging cut. It can be seen that the shapes of the
two surfaces fit well both in the axial and in the feed direction. The slope of the
surface in the feed direction at the end of the process is also predicted properly.
Additionally, the reconstruction of surface topographies was verified for unstable
process conditions. The measured cutting forces can easily be applied in the model in
the presence of chatter. The measured and reconstructed surfaces for a milling process
(Table 2.8 No. 8) with chatter vibration is shown in Fig. 2.27. Due to the necessary
simplifications of the model of system dynamics, some deviations can be seen in the
occurrence of dynamic surface marks. The reconstructed surface roughness is larger
than the measured one, but the characteristics of the chatter marks are similar.
In case of chatter vibration, the presented method provides rather a qualitative
prediction for the machined surface. Depending on the application of the developed
method, however, the reconstruction of workpiece quality for unstable process conditions will not be relevant as long as process monitoring strategies will identify instability, anyway. Especially for the task of process monitoring, the developed method will
aim the understanding of process signals and will provide online evaluation criteria
for quality related information from manufacturing processes.
54
55
profile height [mm]
(a)
]
axial dir
dd
ection [m
fee
m]
n
ctio
[mm
ire
profile height [mm]
(b)
m]
c
axial dir
d
fee
ection [m
m]
(c)
reconstructed
profile height [mm]
profile height [mm]
[m
tion
dire
(d)
reconstructed
microscopy
microscopy
axial direction [mm]
feed direction [mm]
Figure 2.26: Measured (a) and reconstructed (b) surfaces created by a transient vibration for parameter set 1 in Table 2.8. The surface profiles along the tool axis (c) and the surface profile along
the feed direction (d) are plotted along the black lines in (a) and (b) for both cases.
ial
ax
d ir
ec
tio
m]
n
[m
m
d
fee
]
e
dir
ctio
n
[m
(c)
microscopy
reconstructed
ial
d ir
ec
tio
m]
n
[m
m
n[
m
ctio
]
d
fee
e
dir
reconstructed
profile height [mm]
ax
profile height [mm]
profile height [mm]
(b)
(a)
axial direction [mm]
(d)
microscopy
feed direction [mm]
Figure 2.27: The measured (a) and the reconstructed (b) surfaces created by a slight chatter vibration
for parameter set 8 in Table 2.8. The surface profile along the tool axis (c) and the surface profile
along the feed direction (d) are plotted along the black lines in (a) and (b) for both cases.
New Results
Thesis 4
A combined dynamical model of the milling machine is constructed that consists of
a finite element beam model of the finger-like milling tool, a model of the machine
tool structure based on its frequency response measurement, and an experimentally
verified dynamic coupling model between them, while the workpiece-side is considered
to be rigid. This dynamical model is subject to the cutting force which is measured
online during the milling process, and the machined surface profile is generated with
high accuracy numerical simulations of the cutting edge path and by means of Boolean
operations. The advantage of this method is that the surface quality of the workpiece can be predicted online during the cutting process by means of cutting force
measurements.
The method is verified experimentally for the case of forced vibrations of milling
processes for both stationary and transient vibrations: the accuracy of the predicted
surface quality shows good quantitative agreement with the offline direct measurement
of the machined surface. In case of self-excited vibrations, the experiments showed
only qualitative agreement with the predicted structure of the machined surface,
which might still be useful in chatter identification.
Related publication: [77]
56
k
j
hj
vj
Figure 2.28: Stationary chip thickness definition.
2.4
2.4.1
Effects of Large Amplitude Vibrations on Milling
Stability
Velocity Dependent Dynamic Chip Thickness in Milling
Processes
The classical mechanical models developed for stability predictions of milling processes use the assumptions that the amplitude of the stationary vibration caused by
the stationary rotation of the milling tool is negligible relative to the peripheral motion of the cutting edges, and that the feed per tooth can be neglected relative to the
tool radius. Hence, the edge path is modelled by a circular motion.
However, when the periodic forcing drives the system close to resonance, and in
the same time, we use large axial immersion and feed per tooth, then the forced
vibrations can reach a level where the above assumptions are not valid anymore. In
this case, a so-called velocity dependent chip thickness has important quantitative
effect on the onset of chatter (dynamic instability) and more surprisingly, it also has
an unusual qualitative effect on the static stability of the stationary cutting.
A key element of our approach in this subsection is the use of a simple mechanical
model that considers the variation of the chip thickness due to the stationary forced
vibration in the milling process. The same 2 DoF model with straight edged tool
is used in the orthogonal cutting model of milling as described in Sec. 2.2.1 and
Subsec. 2.2.1.1.
To analyse the effect of the forced vibration, the chip thickness is re-defined here,
and we take into account the direction of the velocity of the tool-tip. Usually, the
chip thickness is defined as the size of the chip normal to the velocity vector of the jth
tool-tip vj . The corresponding rake angle is denoted by κ in Fig. 2.28. The cutting
coefficients depend on this angle in an intricate way [78], but our model focuses on
the effect of the chip thickness variation only, and the cutting coefficient variation
on the rake angle is neglected here. While we take into account the variation of the
direction of the velocity vector, we suppose that this variation is not ’too large’, that
is, the chip thickness does not change too much relative to the size of the contact
zone. An even more accurate model will be discussed later in Sec. 2.4.2.
A schematic representation of the geometry of dynamic chip forming at the jth
cutting edge is shown in Fig. 2.29. The displacement vector f̃ of the tool-tip over the τ
57
ideal path
h
RW cos(fj )+y(t)
fj
fz
vj
y(t) x(t)
x(t-t)
y(t-t)
~
f
RW sin(fj )+x(t)
nvj
current path
previouse path
Figure 2.29: Dynamic chip thickness h and the feed-per-tooth fz . Double-line vectors denote velocities, where vj is the tool-tip velocity.
tooth pass period is defined by the difference of the present and the past intersections
of the rake faces and the tool paths:
fz + x(t) − x(t − τ )
f̃ =
,
(2.59)
y(t) − y(t − τ )
where fz is the feed per tooth. Using the tool-tip velocity
RΩ cos(φj (t)) + ẋ(t)
vj =
,
RΩ sin(φj (t)) + ẏ(t)
of the jth tool-tip, the vector nj normal to the tool-tip velocity is given by:
RΩ sin(φj (t)) + ẏ(t)
nj =
.
−RΩ cos(φj (t)) − ẋ(t)
(2.60)
(2.61)
In Eq. (2.60) and Eq. (2.61), the feed velocity N fz Ω/(2π) is neglected since it
is small compared to the peripheral cutting velocity RΩ. In case of real machining process parameters, their ratio is less than 1%. This feed velocity is also small
compared to the vibration velocity components (ẋ and ẏ) if the system is close to
resonance and large forced vibrations occur. Using Eq. (2.59) and Eq. (2.61), the
dynamic chip thickness can be derived by projecting the tool-tip displacement over
tooth pass period to the unit vector normal to the actual tool-tip velocity:
nj
=
|nj |
(fz + x(t) − x(t − τ ))(RΩ sin(φj (t)) + ẏ(t)) − (y(t) − y(t − τ ))(RΩ cos(φj (t)) + ẋ(t))
p
.
(RΩ sin(φj (t)) + ẏ(t))2 + (RΩ cos(φj (t)) + ẋ(t))2
(2.62)
hj = f̃ >
This way, the governing equation can be built, which leads to a nonlinear delaydifferential equation, where a new kind of nonlinearity appears in accordance with the
expression 2.62. Note, that if the velocity [ẋ(t) ẏ(t)]> of the tool centre is negligible,
then we obtain the classical chip thickness expression (see: Eq. (2.2)):
58
nj
≈
|nj |
(fz + x(t) − x(t − τ )) sin(φj (t)) − (y(t) − y(t − τ )) cos(φj (t)),
hj = f̃ >
(2.63)
which is commonly used in traditional chatter models.
Stability Chart
To describe the effect of the cutting parameters, a particular machine-tool/workpiece
configuration was chosen with a standard straight fluted milling tool given in Subsec. 2.2.5.1. During the calculations, 10% radial immersion was used and a tool was
considered with N = 1 tooth only.
We use DDE-BIFTOOL [79] to find the solution of the periodically forced system,
and to check the stability of the resulting periodic solution. In this case, the stability
charts also depend on the feed per tooth parameter in spite of the fact that the cutting
force is considered to be a linear function of the chip thickness. In other words, if
we neglect the velocity of the tool centre and use a linear cutting force model then
the stability of the system does not depend on the feed per tooth. Consequently,
for small values fz of feed per tooth, we obtain the classical stability chart with the
conventional instability lobes in it (see the white domains of Fig. 2.30). The new
type of chip thickness calculation is important only in the case when large amplitude
forced vibrations occur. Such large vibrations can be found easily at resonant spindle
speeds (like around 2253 rad/s in Fig. 2.30). For large feed per tooth, we find a new
type of stability loss leading to a fold bifurcation [80]: an eigenvalue of the Floquet
transition matrix crosses the unit circle at +1 (see Fig. 2.31).
With DDE-BIFTOOL, we also calculated the periodic solutions for the constant
spindle speed Ω=2253 rad/s. A bifurcation diagram was created by plotting the amplitude of the vibration against the axial immersion ap (see Fig. 2.32). The emerging
fold bifurcations give the boundaries of a new type of bistable regions, where two
stable periodic motions are separated by an unstable one. All the three have the
same periodicity as the forcing has. These vibrations have nothing to do with chatter caused by the regenerative effect, since the fold bifurcation points are located
in the shaded stable region of the stability chart of Fig. 2.31, where stable region
means stability against chatter. The continuous lines of Fig. 2.31 present the first
fold bifurcation points only.
If the chip thickness were not velocity dependent, the fold bifurcation could never
occur; if it is velocity dependent then fold bifurcation appears even for small radial
immersion values like 10%. This also means that the velocity dependent chip thickness
model must be used in the region of resonant spindle speeds, that is, where the basic
natural frequency of the structure is just an integer times the spindle speed. On the
other hand, there are certain problems with the application of the velocity dependent
dynamic chip thickness model.
The first problem is originated in the simplified cutting force model. If the amplitude of the forced vibration is large, then the change of the rake angle κ may become
59
ap -
60
Figure 2.30: Stability chart for different feed per tooth parameters with 10% radial immersion.
Shaded region denotes the stable region for the linear model. Continuous lines denote the static loss
of stability in the velocity dependent dynamic chip thickness model.
Im
i
m
1
Re
Figure 2.31: Fold bifurcation: the largest characteristic multiplicator crosses the unit circle at +1.
ap -
Figure 2.32: Bifurcation diagram at Ω=2253 [rad/s], fz =0.5 [mm]. Solid blue lines denote stable
periodic solution, dashed red line denotes unstable periodic solution.
Figure 2.33: Machined workpiece surface created during large amplitude near-resonant vibrations.
significant, resulting in inconsistency with our assumption of constant cutting coefficients (Kt , kr ).
Secondly, for especially large vibrations, the chip thickness can become negative,
which would create an unrealistic negative cutting force in the above model. The
model may seem to be easily adjustable to this effect in the same way as the classical
loss of contact effect is handled in non-linear machine tool vibrations [36], but there
are even more complicated situations to appear here, like the undercut. In this case,
the tooth leaves the surface and can make contact again during the same cutting
phase, and the surface is not formed from the current tooth path only.
A third problem of the model is that the projection in the chip thickness calculation in Fig. 2.29 is correct only for a tooth paths with small curvatures.
To describe the surface correctly, we need a mathematical model that involves
state-dependent delays as discussed in detail in the subsequent section.
2.4.2
State-Dependent Delay Model
In the paper [37], an accurate model of the chip removal in milling is described for large
amplitude vibrations. This model leads to a time- and state-dependent delay model,
which is essentially non-linear even in the delay parameter. The governing equation
and the process of the linearisation around its otherwise analytically unknown periodic
solution are presented and this leads to a somewhat modified stability chart. Still,
the conclusion in [37] is that changes caused by the state-dependency of the delay are
negligible in the practically relevant range of machining parameters. The results of
the previous Section, also manufactured workpiece surfaces like the ones in [53] or the
one in Fig. 2.33 which was generated during our experiments with large amplitude
near-resonant vibrations imply that more detailed analysis of the state-dependent
delay effect might be necessary.
We assume that the system can be described by a flexible tool / rigid workpiece
model. The motion of the centre of the tool in the x − y plane is described by
the general coordinates q = [x y]T (see Fig. 2.34). The differential equation of the
corresponding two degrees of freedom oscillator has the form:
M q̈ + C q̇ + Kq = F .
61
(2.64)
R
y
j–1
tool
x
t
Fj
a Fj
W
j+1
j
v
r
Fj
fj
workpiece
Figure 2.34: Representation of the mechanical model.
Table 2.9: Type of the cutting force functions
type of approximation
F j (hj )
linear
w (k0 + k1 hj )
cubic polynomial power (q ≈ 0.75)
w k1 hj + k2 h2j + k3 h3j
wkhqj
F denotes the resultant of the cutting forces F j , which depend on the machining
parameters. Among these parameters, we pay special attention to the chip width w
(equal to the axial immersion for orthogonal cutting) and chip thickness h (related
to the feed velocity v). In the local coordinate system given by the tangential and
radial directions of the milling tool at the jth cutting edge, the cutting force is given
T
in the form F j = Fjt Fjr .
Table 2.9 summarizes the most frequently used cutting force expressions. They
all describe the cutting force as a linear function of the chip width w and there are
different linear [78], power law [81, 31] and cubic polynomial [82, 36] functions with
respect to the chip thickness hj . These functions are also displayed in Fig. 2.35.
The total cutting force acting on the tool centre is the resultant of the forces
acting on each tooth
F (hj , t) =
N
X
T (φj (t)) F j (hj ),
(2.65)
j=1
where
T (φj (t)) =
cos (φj (t)) sin (φj (t))
− sin (φj (t)) cos (φj (t))
(2.66)
is the projection matrix from the local (cutting-edge-fixed) coordinate system to the
62
Fj
(b)
(a)
(c)
hj
Figure 2.35: Approximations of the cutting force function: a) linear-, b) power-, and c) polynomial
function.
y
Pj-1(t)
tv
R
v
O(t)
r
hj
Pj-1(t-tj)
O(t-tj)
W
fj-1(t-tj)
Pj(t)
fj(t)
x
t
workpiece
Figure 2.36: Accurate model of chip thickness calculation.
x − y coordinate system, N is the number of teeth and φj (t) = tΩ − (j − 1)2π/N is
the angular position of the jth tooth (j = 1, 2, ...N ), where the milling tool angular
velocity is Ω.
2.4.2.1
Improved Chip Thickness Modelling
The chip thickness of the jth tooth is determined by the current tool-tip position and
the previously machined surface which is generated by the previous (j − 1)st tooth
path. Figure 2.36 shows an accurate model of the chip thickness calculation. It can
be seen that during the machining, the chip thickness is defined by the position of
the jth tool-tip and the intersection point of the jth tooth rake face and the (j − 1)st
edge path. The tool-tip position P j (t) is determined by the angular position of the
tooth and the position of the tool centre O(t) = q(t) = [x(t) y(t)]T according to
x(t) − R sin (φj (t))
P j (t) =
.
(2.67)
y(t) − R cos (φj (t))
The above mentioned previous intersection point depends also on the feed displacement vτj as given by
63
P j−1 (t − τj ) =
x(t − τj ) − R sin (φj−1 (t − τj )) + vτj
y(t − τj ) − R cos (φj−1 (t − τj ))
,
(2.68)
where the delay τj is the time needed by the milling tool to get from shadowing the
tool-tip position P j−1 to P j . This is not exactly the time (T = 2π/N Ω) between
two consecutive cuts of the jth and (j − 1)st teeth used in conventional regenerative
theory. The intersection point P j−1 (t−τj ) can also be reached from the tool-tip point
P j (t) by using the chip thickness hj
P j−1 (t − τj ) = P j (t) + hj
sin (φj (t))
cos (φj (t))
=
x(t) − (R − hj ) sin (φj (t))
y(t) − (R − hj ) cos (φj (t))
. (2.69)
The combination of (2.68) and (2.69) provides the relation
x(t − τj ) − R sin (φj−1 (t − τj )) + vτj
y(t − τj ) − R cos (φj−1 (t − τj ))
=
x(t) − (R − hj ) sin (φj (t))
y(t) − (R − hj ) cos (φj (t))
, (2.70)
which are two non-linear scalar equations with respect to two scalar unknowns, the
actual chip thickness hj and the time delay τj . Equation (2.70) can be projected to
the local tangential – radial coordinate system by means of the inverse T −1 (φj (t))
of the projection matrix. The chip thickness is obtained from the equation of radial
projection after a straightforward calculation:
hj = R 1 − cos (τj Ω − 2π/N ) +
(2.71)
vτj + x(t − τj ) − x(t) sin (φj (t)) + y(t − τj ) − y(t) cos (φj (t)) .
The other equation of tangential projection leads to the following implicit equation
with respect to the unknown τj :
R sin τj Ω − 2π/N + vτj + x(t − τj ) − x(t) cos (φj (t))
− y(t − τj ) − y(t) sin (φj (t)) = 0.
(2.72)
This equation expresses that the time delay τj is time and state dependent, due to
the fact that it depends on the current and delayed positions of the tool centre, i.e.,
on x(t), y(t) and x(t − τj ), y(t − τj ), respectively. If it were possible to solve this
equation in closed form, we would obtain the time delay as a function τj (t, xt , yt ),
where xt (s) = x(t + s), yt (s) = y(t + s), s ∈ [−τmax , 0], and τmax ∈ R+ (see also
[37]) represent the present and the past states as a ’tail’ of their time history, that
is, as a function with respect to the past-time parameter s. The time delay also
depends on the feed velocity, however, its effect is observable only in case of small
radial immersion (<10%) and not relevant for higher values of radial immersions as
it is shown in [83].
64
y
gr=0
workpiece
gr=1
ymax
x
gr=0
ymin
hj>0
gh=1
hj<0
gh=0
Figure 2.37: Values of the screen function gr and gh
2.4.2.2
Improved Screen Function Modelling
The jth chip thickness can be calculated from (2.71) and (2.72) for the cutting phase
only, while it is zero for the non-cutting phase. This phenomenon is often described
by using a screen function that depends on the tool actual angular position only
[84, 85]:
1 if φenter < φj (t) < φexit
g (φj (t)) =
,
(2.73)
0
otherwise
where the constant parameters φenter and φexit are the angles where the jth tooth
enters and exits the workpiece material. Clearly, this screen function is approximate in the sense that it is valid for stationary cutting and small-amplitude vibrations only. The improved screen function consists of two parts. The first part
describes tool-workpiece contact period with respect to the radial immersion. The
corresponding screen function gr depends on the tool tip P j vertical position yj (t) =
y(t) − R sin (φj (t)) obtained from: (2.67) relative to the workpiece vertical boundaries ymin and ymax (see: Fig. 2.37). Clearly, the (theoretical) radial immersion is just
R + ymax if ymin < −R. Thus the first part of the screen function has the form
1 if ymin < y(t) − R sin (φj (t)) < ymax
gr (y(t), t) =
.
(2.74)
0
otherwise
The second part gh of the screen function is used to express the fact that toolworkpiece contact appears only in case of positive chip thickness:
1 if hj > 0
gh (hj ) =
.
(2.75)
0 otherwise
The governing equation is formed from (2.64), (2.65), (2.74), and (2.75):
M q̈(t) + C q̇(t) + Kq(t) =
N
X
gr y(t), t gh hj T (φj (t)) F j (hj ),
(2.76)
j=1
where the chip thickness can be expressed in an abstract form hj (t, xt , yt ) or with an
65
obscure classical notation:
hj t, x t , y t ,
x t − τj (t, x(t), y(t), x (t − τj (...)) , y (t − τj (...)) , τj (...)) ,
y t − τj (t, x(t), y(t), x (t − τj (...)) , y (t − τj (...)) , τj (...)) ,
τj t, x(t), y(t), x (t − τj (...)) , y (t − τj (...)) , τj (...)
involving infinite sequences of embedded functions originated in (2.71), and in the implicit definition (2.72) of τj . Since the time delay τj depends on the state coordinates,
the governing equation (2.76) is a non-autonomous state-dependent delay differential
equation (SD-DDE) with time-periodic coefficients of period T = 2π/N Ω and with
time-periodic excitation of the same period T .
2.4.2.3
Periodic Solution
We apply a shooting method to determine the chatter-free T -periodic solutions of the
periodically forced SD-DDE (2.76). The method is a generalization of the approach
used in [86, 87] for non-linear ordinary differential equations (ODEs).
The time integration in the time delayed system can be carried out by selecting a
suitable initial function ϕ(t), t ∈ [−τmax , 0]. We define the time integration over one
principal period T as a non-linear mapping z of the function-space C[0, T ; R2 ] into
itself, where C[0, T ; R2 ] is the space of continuous R2 valued functions on the time
interval [0, T ]. The periodic solution q̄, is the fix point of this non-linear mapping,
i.e.,
q̄ − z(q̄) = 0.
(2.77)
We obtain successive approximate solutions of (2.77) using the following procedure:
1. We take a finite dimensional subspace, C n [0, T, R2 ] of C[0, T, R2 ] where the
time interval [0, T ] is discretized by n time-steps of length ∆t = T /n, and pick an
initial guess q̄ n0 for q̄ in the form
q̄ n0
=
x̄(0) x̄(∆t) . . . x̄(n∆t)
ȳ(0) ȳ(∆t) . . . ȳ(n∆t)
T
.
(2.78)
2. We compute z(q̄ n0 ) (the details of the numerical integration procedure is explained below) and define the error function, or the so-called ’shooting function’
H 2 (q̄ n0 ) =
n
X
n,i
q̄ n,i
0 − z(q̄ 0 )
i=1
66
2
.
(2.79)
3. If H 2 (q̄ n0 ) < then we use q̄ n0 as an approximation of q̄, otherwise we find
= 1, 2... by minimizing H(q̄ n0 )2 iteratively over the subspace C n [0, T, R2 ].
In our calculations, all the codes are built in Matlab. For the time integration
in step 2, a 3rd order Runge-Kutta method with fixed time step is used. In each
integration step we solve (2.72) consecutively for all τj (j = 1, ..., N ), then (2.71) for
all the chip thicknesses hj (j = 1, ..., N ), then (2.65) for the resultant cutting force F .
Having the righ-hand-side of (2.76), we implement the Runge-Kutta method to make a
time step ahead to obtain q̄ n,n+1
and the procedure is repeated again
to obtain q̄ n,n+2
0
0
,
n,n+1
n,n+2
n,2n
n,2n
n
q̄ 0
... q̄ 0
.
and so on, n times, till we obtain q̄ 0 . All these serve z(q̄ 0 ) = q̄ 0
In our calculations in step 3, the fixed point is generated by using a built-in Matlab
function ’fsolve’, which finds the root of the system of nonlinear equations in (2.79),
i.e., solves H 2 (q̄ n0 ) = 0. Function values between time steps are obtained by spline
interpolation; these are needed during the solution of the implicit equation (2.72) and
also in (2.71).
q̄ nj , j
2.4.2.4
Linearization Around the Periodic Solution
In order to investigate the stability of the above numerically determined periodic
solutions q̄(t), the SD-DDE (2.76) is linearized around these periodic solutions with
respect to the new "small" variable
u(t) = q(t) − q̄(t).
(2.80)
Since the system is piecewise continuous due to the screen functions, we can
distinguish a cutting phase where gr = gh = 1, and a flying phase where gr = 0
and/or gh = 0. We start the linearization with the latter case. In the flying phase,
the cutting force is zero, hence the right-hand-side of (2.76) is zero, and the SD-DDE
is simplified to a linear ODE. Thus, we are left with the linearized equation
M ü(t) + C u̇(t) + Ku(t) = 0,
(2.81)
of the two DoF damped oscillator.
In the cutting phase, the system is non-linear due to the state-dependency. The
linearization procedure follows the one presented in [88]. The associated linear system
is given by
M ü(t) + C u̇(t) + Ku(t) =
N
X
T (φj ) D1 F j (h̄j ) D2 h̄j u(t)
(2.82)
j=1
+ D3 h̄j u t − τj (t, q̄ t ) − q̄˙ t − τj (t, q̄ t ) D2 τj (t, q̄ t )ut + D4 h̄j D2 τj (t, q̄ t )ut ,
where Di denotes the derivative with respect to the ith argument, and the arguments
that are not displayed in the above linear form have the following sequence in the
chip thickness function:
67
Table 2.10: Derivatives of the cutting force functions
type of approximation
F j (hj )
D1 F j h̄j
linear
w (k0 + k1 hj )
wk1
cubic polynomial
w k1 hj + k2 h2j + k3 h3j
w k1 + k2 2h̄j + k3 3h̄2j
power (q ≈ 0.75)
wkhqj
wkq h̄q−1
j
hj t, q(t), q (t − τj (t, q t )) , τj (t, q t )
(2.83)
and the overbar refers to its value at the periodic solution determined above.
The derivative
D1 F j (h̄j ) =
∂F j
(h̄j ),
∂hj
(2.84)
of the force function can easily be derived for the commonly used cutting force functions (see Table 2.10).
The partial derivatives Di hj of the chip thickness (2.71) are gradients with respect
to the 2nd and 3rd arguments (i = 2, 3):
∂hj
t, q̄(t), q̄ (t − τj (t, q̄ t )) , τj (t, q̄ t )
∂q(.)
h
i
=
− sin (φj (t)) − cos (φj (t)) ,
D2 h̄j =
∂hj
t, q̄(t), q̄ (t − τj (t, q̄ t )) , τj (t, q̄ t )
∂q(. − τj )
h
i
= sin (φj (t)) cos (φj (t)) = −D2 h̄j ,
D3 h̄j =
(2.85)
(2.86)
while D4 hj is simply the derivative with respect to its 4th argument that is the time
delay:
∂hj
t, q̄(t), q̄ (t − τj (t, q̄ t )) , τj (t, q̄ t ) = RΩ sin (τj Ω − 2π/N ) + v sin (φj (t)) .
∂τj
(2.87)
The partial derivative D2 τj of the time delay τj with respect to the function qt is a
so-called Frechet-derivative [88], which is the infinite dimensional counterpart of the
gradient
D4 h̄j =
D2 τj =
∂τj
.
∂q t
(2.88)
However, the term D2 τj (t, q̄ t )ut can only be computed from the variational system
of the implicit equation for the time delay τj in (2.71). With the state vector q(t),
68
this equation can be rewritten as
R sin (τj Ω − 2π/N ) + vτj cos (φj (t))
h
i
+ cos (φj (t)) − sin (φj (t)) (q(t − τj ) − q(t)) = 0.
(2.89)
Its implicit differentiation at the periodic solution leads to:
RΩ cos (τj (t, q̄ t )Ω − 2π/N ) D2 τj (t, q̄ t )ut
h
i
+v cos (φj (t)) D2 τj (t, q̄ t )ut + cos (φj (t)) − sin (φj (t))
× (u (t − τj (t, q̄ t )) − q̄˙ (t − τj (t, q̄ t )) D2 τj (t, q̄ t )ut − u(t)) = 0.
(2.90)
(2.91)
Since it is a linear equation with respect to the term D2 τj (t, q̄ t )ut in question, we
can determine its solution in the form:
D2 τj (t, q̄ t )ut = pj (t) u (t − τj (t, q̄ t )) − u(t) ,
(2.92)
where the T -periodic coefficient can be expressed as
pj (t) =
(2.93)
h
i
− cos (φj (t)) − sin (φj (t))
h
i
.
˙ − τj (t, q̄ t ))
RΩ cos (τj (t, q̄ t )Ω − 2π/N ) + v cos (φj (t)) − cos (φj (t)) − sin (φj (t)) q̄(t
Combining (2.82) with (2.85), (2.86), (2.87) and (2.92), and collecting the variables, the linear equation of motion is given by
M ü(t) + C u̇(t) + Ku(t) =
(2.94)
N
X
˙ − τj (t, q̄ t ))
T (φj ) D1 F j (h̄j ) −D2 h̄j + pj (t) D4 h̄j + D2 h̄j q̄(t
j=1
× u(t − τj (t, q̄ t )) − u(t) .
This linear DDE has T -periodic coefficients and also T -periodic time delay but
it does not contain state-dependent delay anymore. Note, that it is valid for the
cutting phase only. Whether the linear equation of the cutting phase (2.94) or that
of the flying phase (2.81) is to be used depends on the forced periodic solution q̄(t)
substituted into the screen functions (2.75) and (2.74). All these are expressed by the
form
69
M ü(t) + C u̇(t) +
K+
N
X
!
W j (t) u(t) =
j=1
N
X
W j (t)u (t − τj (t, q̄ t )) ,
(2.95)
j=1
where the jth T -periodic directional force coefficient matrix [61] is
W j (t) = gr (q̄(t), t) gh h̄j T (φj ) D1 F j h̄j
−D2 h̄j + pj (t) D4 h̄j + D2 h̄j q̄˙ (t − τj (t, q̄ t )) .
(2.96)
This linearized form of the governing equation is a generalization of the one derived
in [37].
2.4.2.5
Stability by the Semi-Discretization
The variational system (2.95) is a parametrically forced linear DDE with time periodic
delays. The stability of such systems can be investigated by the semi-discretization
method. The theory is developed and described in details in [13, 89], while some
calculations and results for machine tool vibrations can be found in [90, 91]. Note,
that the forced periodic motion q̄(t) must be calculated in advance by means of the
algorithm described in Section 2.4.2.3 in order to perform the stability analysis by
means of the semi-discretization method.
2.4.2.6
Real-Case Parameters
We carry out the calculations using the real-case modal parameters in the mechanical
model given in Subsec. 2.2.5.1, where the natural frequencies of the two first bending modes of the tool are very close to each other (ωnx = 4537.3 rad/s and ωny =
4537.6 rad/s) due to cylindrical symmetry. During the calculations, we consider a
straight fluted tool of radius R = 4 [mm] and of number of teeth N = 2. The workpiece material is AlMgSi0.5 aluminium alloy and the cutting force vector function
F j (hj ) is chosen to be linear (see 2nd column of Table 2.10). The corresponding
cutting coefficients k0 and k1 were identified experimentally by the authors of [57]
(see: Subsec. 2.2.5.1)
N
N
644 000
, k1 =
.
k0 = 0
236 992
m
m2
The workpiece vertical boundaries ymin = −∞ and ymax = −1.2 mm (see Eq.
(2.74)) describe an up-milling process with 35% radial immersion. The feed-rate, the
axial depth of cut and the cutting speed are varied during the stability calculations.
2.4.2.7
Stability Chart
First, we check our model in the limit case when the effect of the state-dependent
delay is negligible. The stability chart is calculated in two different ways. On one
70
w - axial immersion [mm]
2
1.5
unstable
1
0.5
stable
0
1000 wn
2N
wn 2500
N
W - angular velocity [rad/s]
1500
2000
3000
Figure 2.38: Comparison the stability charts of the standard (shaded region with continuous line) and
the new state-dependent delay model (dashed line with crosses). Parameters: N = 2, ymin = −∞,
ymax = −1.2 [mm], v = 1 [µm/s]
hand, we calculate the stability boundaries of the standard model as described in [80]
when the forced vibration due to periodic change of the cutting force during milling
assumed to be small (see shaded area with the continuous boundary lines in Fig. 2.38).
On the other hand, we also calculate the stability boundaries (see dashed line with
crosses in Fig. 2.38) with our improved model containing the state-dependent delay;
but the feed velocity v = 1 µm/s is set to such an (unrealistically) small value that
the chip thickness is very small just like the cutting force is; this leads to a very small
forced oscillation, consequently, the state-dependency and the time-periodicity of the
time delay are negligible. The two stability charts are to be almost identical for small
feed rates. We can check numerically the difference between the two charts: Fig. 2.38
shows tiny gaps between the boundary lines, especially in the thin resonant regions
(e.g. at ωn /N ≈2000 rad/s), which proves that the resonant regions are very sensitive
for the state-dependency of the delay even in case of very small feed rates.
Stability charts for non-negligible state-dependent delays are presented in the next
subsection, where additional bifurcations will also be detected.
2.4.2.8
Fold Bifurcation
Now, we check the results of our model for a set of larger realistic feed velocities
v = 20, 30, 40, 50 mm/s. In case of constant time delay τ = T = 2π/N Ω, the
commonly used feed-per-tooth parameter is varying from 0.079 to 0.020 mm/tooth
in case of v = 20 mm/s and from 0.196 to 0.051 mm/tooth in case of v = 50 mm/s
while the angular velocity is in the range Ω = 800 − 3100 rad/s.
The state-dependent delay model (2.76) is essentially non-linear, and the stability
boundaries can be calculated from its linearized form (2.95), but the forced periodic
solution in it will change as the cutting parameters are varied in the stability chart.
We solve this task by using a continuation method similar to the one described in
[86, 87].
During the calculations, the angular velocity Ω is fixed, and the chip width (axial
immersion) w is used as a continuation parameter. The calculation is initialized at
71
w - axial immersion [mm]
2
1.5
v - feed velocity
20 mm/s
30 mm/s
40 mm/s
50 mm/s
1
0.5
0
wn
2400
N
W - angular velocity [rad/s]
1800 1900 2000 2100 2200
2500 2600
Figure 2.39: Stability chart computed by the semi-discretization method (shaded region with continuous thin line) and the limit of first bifurcation points wb of the state-dependent delay model
(thick lines). Parameters: N = 2, ymin = −∞, ymax = −1.2[mm].
the small (unrealistic) chip width w =0.01 mm and with zero initial guess matrix q̄ n0
(for the periodic solution q̄(t), see Section 2.4.2.3). Due to the small amplitude, the
periodic solution can be found in a few iterations. Then we increase the chip width
by ∆w the value of which depends on the number of iterations we had to use to find
the forced periodic solution; this varies in the range of 0.002–0.04 mm. At the next
chip width, the next initial guess vector is calculated by a predictor-corrector method,
then we iterate to the new slightly modified forced periodic solution. This periodic
solution is substituted into the variational system (2.95) and the semi-discretization
method provides its stability. Following this procedure for increasing chip widths,
the stability boundary can be defined at the first chip width point wb (Ω) where the
instability of the periodic solution is detected for the fixed angular velocity. The
subscript b refers to "bifurcation" since the periodic solution will bifurcate at the
stability boundary. We will also refer to this stability boundary as "limit of first
bifurcation points".
The stability boundaries of the standard model and the limit of first bifurcation
points in our improved model are plotted in Fig. 2.39 for different feed velocities v.
It can be seen, that in case of non-resonant angular velocities, the stability limits are
almost the same for the two models, the state-dependent delay has no effect at all.
However, when the feed velocity v is increased at resonant angular velocities (Ω =
ωn /k, k = 1, 2...), that is, for parameters within the pockets between the standard
lobes, the limit of first bifurcation points appear at much smaller chip width values
than in case of the standard model. This means that the standard model overestimates
the region of stability between the stability lobes.
In Fig. 2.40, the bifurcation diagram of the forced periodic motion is plotted
for Ω = 2230 [rad/s] and v = 50 mm/s, where the commonly used feed-per-tooth
parameter is 0.07 mm/tooth. Clearly, fold bifurcation occurs at the limit of first
72
bifurcation point wb = 0.385 mm. At this point, an eigenvalue of the transition
matrix calculated in the semi-discretization method crosses the unit circle at +1.
The shaded region of Fig. 2.40 between wb − wc also shows that there is a region of
the chip width w, where the system has at least one stable and one unstable forced
periodic solution. In this shaded region, the stable forced periodic motion is stable
only locally. This means, that an unexpected perturbation may push the system to
the "other side" of the identified unstable periodic motion and it cannot get back to
the stable stationary cutting again. This is why this zone is also called “unsafe zone”;
it is unsafe in the sense that certain perturbations may result in global instability. It
is still unknown for us what kind of vibration occurs outside of the unstable branch.
Similar phenomenon occurs in case of orhogonal cutting processes at the stability
limits where sub-critical Hopf bifurcations occur (see the Hopf bifurcation point in
the straight line in Fig. 2.40), and an unstable periodic motion appears around the
stable stationary cutting. "Outside" this unstable motion, the so-called "fly-over
effect" takes place: the tool leaves the workpiece and then re-enters it. This is a nonsmooth non-linearity which generates a non-smooth bifurcation in the branch of the
unstable periodic solution (see details in [92]). The authors found that the unstable
periodic motion separates the domain of attractivity of the stable stationary cutting
and the domain of chaotic large amplitude vibrations. In our model we expect similar
behavior: the "fly-over effect" also exists here in the sense that one cutting edge of
the milling tool leaves and re-enters the workpiece within one "cutting phase" (see
Fig. 2.41).
If we choose the axial immersion below the first bifurcation point (w < wb ), it
might be critical if it falls in the unsafe zone (wc < w < wb , see Fig. 2.40), due to the
above mentioned phenomena. If we want to stay in the parameter domain of globally
stable stationary cutting, we should know the width of this zone along the w axis, that
is, we should calculate how far the unstable forced periodic solution folds back. The
lower limit wc of the unsafe zone also denoted in Fig. 2.40 gives an upper estimation
for this. At this parameter point our continuation method loses its validity, due to
new types of non-linearities showing up in the milling process, which are discussed
in the next Subsection. One of these is the above described fly-over effect. In these
situations the above described equations (2.71) and (2.72) provide non-realistic chip
thicknesses and cutting forces.
Still, as it is also presented in Fig. 2.40 as “hypothetic region”, a stable vibration
generated by the fly-over effect is likely to exist outside the unstable branch of periodic
motions. If it is so, the unsafe zone could also be a so-called bi-stable zone where
stable stationary cutting and stable chaotic-like oscillation co-exists.
2.4.2.9
Limitations of the Model
In case of the so-called fly over effect the amplitude of the vibration is so large that
the tool leaves the material during the cutting phase. This is a kind of self-regulated
interrupted cutting as opposed to the parametrically interrupted cutting caused by
the rotation of the milling tool. The chip thickness calculation (2.71) is not valid
in that case, and multiple previous tool tip paths should be considered in a correct
73
74
chaotic vibration
with fly-over effect
Hopf bifurcation
Fold
bifurcation
wb
branch for constant
delay model
non-smooth
bifurcation
6
0.3
2
-2
0
-4
-6
-6
-2
-4
-2
0
2
6
4
6
-4
4
-6
-6
2
4
6
0
-6
-6
-4
-2
0
2
4
-2
6
6
-4
4
-6
-6
2
-2
0
2
4
6
4
2
-4
0.2
-4
6
0
-2
0
-6
-6
-4
-2
0
2
4
6
-2
-4
6
-6
-6
4
-4
-2
0
2
4
6
branch for state
dependent delay model
2
0
-2
-4
-6
-6
-4
-2
0
2
4
hypothetic region of
large amplitude vibration
with fly-over effect
6
6
0.1
2200 wn
2400
N
W - angular velocity
[rad/s]
0
2
-4
0.2
0
-2
4
-2
wc
wb
-4
6
0
0.3
0
4
0
2
1
0.5
6
2
1.5
w - axial immersion [mm]
w - axial immersion [mm]
4
unsafe zone
2
4
2
0.1
0
-2
-4
-6
-6
-4
-2
0
0
2
4
6
2
4
6
6
4
2
0
-2
-4
-6
-6
0
0
-4
0.2
-2
0.4
0.6
0.8
1
1.2
||x,y|| - average amplitude [mm]
1.4
1.6
Figure 2.40: The bifurcation diagram of the milling process and the path of the cutting edges. Black
circles denote stable processes, white circles denote unstable processes. The straight line represents
the system with constant time-delay, and its Hopf bifurcation. Continuous lines represent stable
branches, dashed lines represent unstable branches. Dark clouds denote the hypothetic chaotic
region. Parameters: N = 2, Ω = 2270 [rad/s], ymin = −∞, ymax = −1.2 [mm], v = 50 [mm/s],
τ v = 0.07 [mm/tooth].
computational error
computed
chip thickness
real chip thickness
Figure 2.41: Representation of the fly-over effect.
path of the
tool-tip
P(t-t)
dt
>1
dt
P(t)
workpiece
Figure 2.42: Computational error of the chip thickness in case of dτj /dt > 1
model (see, e.g., [93], and also Fig. 2.41).
Equation (2.71) does not give the proper chip thickness either when dτj /dt > 1,
i.e., t − τj is decreasing – similarly to undercut. In this case, the computed chip
thicknesses overlap (see the black region in Fig. 2.42).
There is a third problem, when the implicit equation for the state-dependent time
delay (2.72) has multiple solutions. This situation is represented in Fig. 2.43 when
the solution of (2.76) is not unique – this is also called "multiple chip thicknesses". In
mechanical sense, the problem is that two distinct contact regions may appear on the
rake face between the chip and the tool tip. Our traditional cutting force calculations
lose their sense and further mathematical conditions would be needed to guarantee
the uniqueness of the solution.
In this study we managed to follow the unstable forced periodic solution till one
of the above mentioned three conditions are violated. This limit point provides an
upper estimation of the possible turning point of the periodic solution, which would
be the lower boundary of the unsafe zone. In Fig. 2.44 the first bifurcation points
and the limit points are plotted around the 1st resonant cutting speed.
75
76
tool-tip
multiple
chip
workpiece
tj
multiple
solution
T
t
Figure 2.43: Representation of the multiple solutions of the time delay τj in space and in time
domain.
w - axial immersion [mm]
2
1.5
fly over effect
1
0.5
multiple chip
thicknesses
0
1800
1900
wn
2400
N
W - angular velocity [rad/s]
2000
2100
2500
2600
Figure 2.44: Limit points around the 1st resonant cutting speed (white circles: fly over effect, black
circles: multiple chip thickness), the first bifurcation point wb (dashed line with crosses) and the
stability chart computed by the semi-discretization method (shaded region with continuous line).
Parameters: N = 2, ymin = −∞, ymax = −1.2 [mm], v = 50 [mm/s].
New Results
Thesis 5
In case of milling with near-resonant spindle speeds, the amplitudes of the forced
vibrations can be very large relative to the tool diameter. It is shown, that the stable
parameter regions predicted by linear time delayed mechanical models do not provide
conservative estimates for the selection of the milling parameters due to the existence
of an unusual bistable parameter domain within the near-resonant zone. The corresponding fold-type bifurcation is identified by two improved non-linear mechanical
models.
The first model uses the chip thickness variation projected to the velocity vector
of the tooltip in contrast to the conventional approximation of the chip thickness
variation by means of the approximated circular path of the tooltip. The second
model is the most accurate by describing the exact path of the tooltip leading to a
state-dependent delay model with an implicit algebraic expression of the delay.
A numerical method was developed for state-dependent delay differential equations to determine the periodic forced vibration and to follow it by continuation as
one bifurcation parameter varies. Both models predict the first fold bifurcation of
stationary milling in the near resonant zone.
Related publications:
[94, 95, 96, 97, 98]
77
2.5
Process Damping in Milling
Accurate experimental verifications of the stability charts for the high speed milling
are given in e.g. [99, 56, 39]. However, the experimentally determined stability
boundaries are usually underestimated by the theory of high-speed milling at low
spindle speeds [41, 42, 43, 44]. This is the inclusions of the process damping effect,
which is often used to improve these models: an additional damping proportional to
the time delay is added to the structural damping to characterize the contact force
at the tiny region between the flank face and the workpiece ([40, 45, 46, 41, 44]).
In this section, the classical idea of the process damping is used to improve the
milling model, which leads to a system with a non-smooth (unilateral) damping.
2.5.1
Process Damping
The investigated mathematical model of the milling process is the same as the one
described in Sec. 2.2.1 and Subsec. 2.2.1.1 but an additional force component is also
applied to describe the so-called process damping.
In the presence of flank wear, the process damping is explained as a contact force
between the flank face and the workpiece (see: Fig. 2.45). As the tooth penetrates
the workpiece, the material approaching the cutting edge is separated into two parts.
In the upper part, the material is removed as a part of the chip while in the lower
part, the material that cannot move upward, is compressed under the flank face of
the tooth. This elastic-plastic deformation creates a pressure p under the worn flank
face leading to a resultant contact force Fproc
having radial component which comes
j
from the pressure and tangential component due to the friction between the material
and the worn flank face. Following [45, 46], suppose that the force Fproc
is linearly
j
proportional to the total volume V of the displaced material, so the process damping
force is given as
µ
proc
Fj = gmat (φj (t)) T (φj (t))
k V,
(2.97)
1 sp
where ksp is a specific volumetric coefficient [46], µ is the dynamic friction coefficient
and the screen function gmat is the same as the one described after Eq. (2.5). The
above mentioned contact force occurs only in case of penetrating tooth, that is when
the radial component of the tool-tip velocity
vr = e>
r (φj (t)) q̇(t)
(2.98)
is less than the radial component of the feed velocity v, as it is demonstrated in
Fig. 2.45: er is the radial unit vector. This behaviour is described in Eq. (2.97) by
an additional velocity-dependent screen function
1
vr < v sin (φj (t))
gproc (vr ) =
(2.99)
0
otherwise.
If the vibration velocity q̇(t) is small compared to the peripheral component of the
78
vr v
vr
RW
v
lw
r
proc
Fj
t
p
RW
r
V
t
Figure 2.45: Process damping is created by the contact force under the worn flank face. The force
appears only in case of penetrating tool (left) and it is zero otherwise (right).
cutting velocity RΩ and its wave length is much larger than the characteristic size lw
of the flank wear, then the volume V of the displaced material can be approximated
as (see [46])
2
l vr − v sin (φj (t))
gproc (vr ) ,
(2.100)
V ∼
= −w w
2
RΩ
where −(vr − v sin (φj (t)))/(RΩ) provides the tangent of the penetration angle of
the triangular section of the displaced material (see Fig. 2.45). The negative sign
shows that the volume is positive in case of negative penetration velocity, it is zero
otherwise.
The resulting cutting force acting on the tool is the sum of process damping forces
Fproc
and the previously calculated cutting forces Fcut
derived from the dynamic chip
j
j
thickness in Sec. 2.2.1:
F(t) =
N
X
proc
Fcut
(t) .
j (t) + Fj
(2.101)
j=1
2.5.2
Governing Equation with Process Damping
The combination of the equations (2.1 - 2.6) and (2.97-2.101) leads to the governing
equation
Mq̈(t) + (C + Cproc (q̇(t), t)) q̇(t) + (K + W(t)) q(t) = Fst (t) + W(t)q(t − τ ) (2.102)
where the process damping matrix is expressed as
Cproc (q̇(t), t) = w
N
X
j=1
l2
µ >
gmat (φj (t))gproc (vr )T(φj (t))
er (φj (t)) ksp w , (2.103)
1
2RΩ
and the stationary part of the resultant cutting force is also modified by the process
damping:
79
τv
v
Fst (t) = W(t)
+ Cproc (q̇(t), t)
.
0
0
(2.104)
The process damping matrix is time-periodic and velocity-dependent. It is a piecewise constant function of the velocity due to the screen function gproc (vr ), so the
governing equation is non-linear, similarly to the case of the turning models with
process damping [45, 46]. In contrast to the turning process, however, the system is
periodically excited by the τ -periodic stationary cutting force Fst (t), which generates
a forced periodic motion. Clearly, this forced periodic motion will also contribute to
the stationary cutting force through the process damping (see Eq. (2.104)).
2.5.3
Short-Delay Representation
The volume of the displaced material can also be computed without considering long
wave lengths for the vibration as Eq. (2.100) does in accordance with [45, 46] by
considering a triangular shaped compressed material volume. Following [100, 9], this
volume can be expressed by means of an integral of the variation of the radial position
of the tool-tip along the flank wear in the form
Z lw
v sin(φj (t))
s
−
q(t)
+s
ds
(2.105)
e>
(φ
(t))
q
t
−
V =w
j
r
RΩ
RΩ
0
if we consider continuous penetration for small vibrations, that is, gproc (vr ) ≡ 1. The
time instant t − s/RΩ defines the time when the tool-tip was at the position s, where
s is the arc-length along the flank face starting from the tool-tip. In this approach, we
have a term with distributed delay instead of the process damping term Cproc q̇. The
form of the corresponding governing equation is similar to the one, which is derived
in [101], where another short-regenerative effect in the milling model is applied with
the chip-flow on the rake face.
The time delayed term in Eq. (2.105) can be approximated by its Taylor series for
variable s as follows:
s 2
s
s
= q(t) − q̇(t)
+ q̈(t)
+ h.o.t..
(2.106)
q t−
RΩ
RΩ
RΩ
Assume that the flank wear is small, that is, the angular position φj (t) can be
assumed to be a fixed value for the whole flank face, then the compressed material
volume is then given by
Z lw
s 2
s
lw2 v sin(φj (t))
>
+ wer (φj (t))
−q̇(t)
+ q̈(t)
+ h.o.t. ds
V =w
2
RΩ
RΩ
RΩ
0
lw2 v sin(φj (t))
lw2
lw3
>
>
=w
− wer (φj (t))q̇(t)
+ wer (φj (t))q̈(t)
± h.o.t..
2
RΩ
2RΩ
3(RΩ)2
(2.107)
It can be seen that the first order approximation of the compressed volume is iden80
tical to Eq. (2.100), which means that the traditional process damping effect can be
treated as a first order approximation of a mathematically more precise short-delay
representation Eq. (2.105).
If we also consider the second order term in Eq. (2.107), then the coefficient of
the acceleration will modify the mass matrix proportionally to the square of the time
delay, that is, inversely proportionally to the square of the spindle speed Ω. Similar
idea can be found in [41] where the modelling of process damping is improved by
adding a correction not only to the ’velocity’ term, but also to the ’inertial’ term.
Note, that if we use the approximated compressed volume and apply the third
order terms or even higher order ones, then third and higher order derivatives of the
general coordinate will appear in the governing equation, respectively. It is shown
in [9], that this kind of approximation is divergent and cannot be used for stability
calculations for higher than second order. However, the stability of a short-delay
model can also be determined by the semi-discretization method and this model
contains the accurate description of the compressed volume.
2.5.4
Periodic Solution
First, the chatter-free forced τ -periodic solution has to be computed in order to
perform its stability analysis. It is computed by the shooting method developed in
Subsec. 2.4.2.3. If we substitute the τ -periodic solution q̄(t) in Eq. (2.102), then the
time-delayed terms cancel, and a modified governing equation can be provided for the
periodic solution in the form
τ
v
v
¨ (t) + C + Cproc (q̄(t),
˙
˙ + Kq̄(t) = W(t)
˙
Mq̄
t) q̄(t)
+ Cproc (q̄(t),
t)
. (2.108)
0
0
This is an non-linear ordinary differential equation. Instead of simulating Eq. (2.108)
˙
for a long time with an arbitrary initial conditions q̄(0), q̄(0),
we determine the
˙
˙ ).
initial conditions corresponding to the periodicity: q̄(0) = q̄(τ ) and q̄(0)
= q̄(τ
Based on these conditions, an iterative method is applied leading to the numerical
determination of the corresponding forced periodic solution (see [94, 86]).
2.5.5
Linearization and Stability
In order to investigate the stability of the above determined periodic solution q̄(t),
Eq. (2.102) is linearized around q̄(t) with respect to the small perturbation variable
u(t) = q(t) − q̄(t)
(2.109)
leading to
˙
Mü(t) + C + Cproc (q̄(t),
t) u̇(t) + (K + W(t)) u(t) = W(t)u(t − τ )
(2.110)
This is a system of linear DDEs with time-periodic coefficients, since the directional
81
force coefficient matrix W and the process damping matrix Cproc are both τ -periodic:
W(t) = W(t + τ )
˙
˙ + τ ), t + τ ).
Cproc (q̄(t),
t) = Cproc (q̄(t
(2.111)
The stability of such systems can be investigated by the semi-discretization method
[80].
2.5.6
Case Study: Effect of Flank Wear Length on Stability
Calculations are performed for real-case modal parameters identified experimentally
in [102]. Two dominant modes are selected, where the natural frequencies are 995.21
Hz and 922.78 Hz, the damping factors are 3.4% and 2.4%, and the modal stiffnesses
are 12.8 N/µm and 19.6 N/µm in the principal x and y directions, respectively. The
modal matrices are
N
Ns
12.83
0
328.2
0
141.1
0
, K=
.
M=
g, C =
0
19.62 µm
0
583.4
0
162.4 m
During the calculations, we consider a straight fluted tool of radius R = 10 mm
with number of teeth N = 2. The workpiece material is AlMgSi0.5 aluminium
alloy and the corresponding cutting coefficients are Kt = 777.9N/µm and Kr =
109.2N/µm, which were identified experimentally in [60]. To characterize the process
damping, the specific volumetric coefficient is ksp = 17670N/m3 , which was measured for aluminium alloy in [46]. A typical value µ = 0.3 of the dynamic friction is
used [43]. The entrance angle φenter = 0◦ and the exit angle φenter = 90◦ describe an
up-milling process with 50% radial immersion. A constant feed per tooth fz = 0.1
mm is chosen, which determines the feed velocity v = fz /τ .
Figure 2.46 shows the stability charts, where the stability boundaries are plotted
as a function of the spindle speed Ω and the axial immersion w for different lengths
of flank wear lw . Cutting speed is represented inversely in the horizontal axis in order
to enlarge the deformation of the lobe-structure at the low cutting speed domain. In
Fig. 2.46a, the velocity-dependency is neglected, that is gproc ≡ 1 if the velocity of
the tool vibration is small. Similar results can be found in [43] where even the timeperiodicity of the process damping coefficients were averaged and process damping
was considered to be constant.
On contrary, the stability charts in Fig. 2.46b are constructed by the time-periodic
and velocity-dependent process damping matrix Eq. (2.103). In case of zero flank
wear, we obtain the standard lobe-structure without process damping even in the low
cutting speed domain. Increasing the flank-wear lw slightly, the stability boundaries
move upwards just as it is in the case of the time-periodic but velocity-independent
process damping matrices. For large lw , however, the pattern of the lobe-structure
changes radically for the time-periodic velocity-dependent process damping matrix.
While the critical axial immersion values at the stability limits are about the same,
the lobes are shifted in an intricate way with respect to the spindle speed parameter
(see Fig. 2.46b).
In case of tools with helical cutting edges, appropriate axial immersion values can
82
w [m] - axial immersion
0.05
lw - 0.0 [mm]
lw - 0.2 [mm]
lw - 0.3 [mm]
0.04
C proc (t )
a)
0.03
0.02
0.01
0
5000
2000
1300 1000
800 700
600
500
450
400 375 350
325
Ω[rad/s] - spindle speed
w [m] - axial immersion
0.05
lw - 0.00 [mm]
lw - 0.10 [mm]
lw - 0.15 [mm]
0.04
lw - 0.2 [mm]
lw - 0.3 [mm]
lw - 0.5 [mm]
Cproc (p(t), t )
b)
0.03
0.02
0.01
0
5000
2000
1300 1000
800 700
600
500
450
400 375 350
325
Ω[rad/s] - spindle speed
Figure 2.46: Stability boundaries for different size of flank wear. b) time-periodic velocity-dependent
process damping a) time-periodic process damping
be chosen as shown in Sec. 2.2. Close to these appropriate axial immersion values, the
forced vibrations are very small and the velocity-dependent screen function is identical
to 1 during the whole cutting phase, that is, gproc ≡ 1 according to Eq. (2.99), and
consequently, the stability condition is likely to be similar to that of Fig. 2.46a. If
the axial immersion values are far from the appropriate ones, the forced vibrations
become significant and the time-periodic velocity-dependent process damping matrix
results stability charts similar to that of Fig. 2.46b.
83
New Results
Thesis 6
An improved model of the process damping effect is developed for milling which was
originally introduced to explain the improved stability properties of cutting at low
cutting speeds. This model takes into account the ratio of the feed velocity and the
cutting speed, which is traditionally neglected. This leads to a non-linear phenomenon
due to the fact that the process damping effect at the flank face is switched off when
the stationary forced vibration results in large vibration velocities of the tool-tip while
moving outwards of the workpiece material. This way, an accurate stationary forced
vibration can be calculated and its linear stability is different from the traditional
stability boundary of stationary milling. The intricate structure of the stability chart
for milling models with the improved process damping description is explored.
It is shown, that the mathematically correct description of the widely used process damping model is equivalent to the description of the so-called short regenerative
effect that contains distributed delay terms. It is proven, that the short regenerative
model linearised with respect to the time delay leads to the traditional process damping model, and some other higher-order approximations can also be derived as its
special cases.
Related publication:
[103]
84
Chapter 3
Multi-Dimensional Bisection Method
3.1
Introduction
During the stability calculations of the milling process the semi-discretization method
was used, in which the parameter set of the two main machining parameters – the
spindle speed Ω and the axial depth of cut ap – were determined for which the largest
Floquet multiplier µ equals to one. Mathematically, it is a root-finding problem of a
scalar non-linear function in a two dimensional parameter space of (Ω, ap ):
f (Ω, ap ) := max |µi (Ω, ap )| − 1 = 0.
i
(3.1)
Another well known stability calculation method is the so-called multi-frequency
solution. In this case, the stability boundary lines are calculated in a three dimensional parameter space – the spindle speed Ω, the axial depth of cut ap and the
chatter frequency ωc – by means of the real and the imaginary parts of a truncated
Hill’s infinite matrix [104]:
f1 (Ω, ap , ωc )) := Re (det(G(Ω, ap , ωc ))) = 0,
f2 (Ω, ap , ωc )) := Im (det(G(Ω, ap , ωc ))) = 0.
(3.2)
The computational time of these methods can radically be decreased if we evaluate
the functions Eq. (3.1) or Eq. (3.2) only around the boundary lines. In order to achieve
this goal, we generalize the bisection method for higher dimensions.
The bisection method – or the so-called interval halving method – is one of the
simplest root-finding algorithms which is used to find zeros of continuous non-linear
functions. This method is very robust and it always tends to the solution if the
signs of the function values are different at the borders of the chosen initial interval.
Unfortunately, it has only linear convergence, more precisely, it doubles the accuracy
with each iteration, which is relatively slow. Hence, it is usually used to find only
a proper initial value for alternative root finding methods which have better rate of
convergence (e.g.: the Newton’s method). A big advantage of the bisection method
is, however, that it can be applied for non-differentiable continuous functions, too.
85
f(x)
x b3
x b2
x b1
x a2 x a3
x
}
}
x a1
}
}}
refined intervals
Figure 3.1: Schematic representation of the iteration by bisection.
Geometrically, root-finding algorithms of f (x) = 0 find one intersection point
of the graph of the function and the axis of the independent variable. In many
applications, this 1-dimensional intersection problem must be extended to higher
dimensions, e.g.: intersections of surfaces in a 3D space (volume) or the stability
calculation as explained above, which can be formulated as the solution of a system
of non-linear equations. In higher dimensions, the existence of multiple solutions
become very important, since the intersections of two surfaces can create multiple
intersection lines.
In this chapter, we give a description and categorization of the root-finding problems in higher dimensions and we create a generalized version of the bisection method.
In Sec. 3.2, the bisection method is defined and generalized for multiple solutions. In
Sec. 3.2.2, the basics of the generalization are detailed and the possible solutions of
the occurring problems are analyzed. The steps of the Multi-Dimensional Bisection
Method are presented in the form of a flowchart (Fig. 3.4). Some hints for the numerical method are given in Sec. 3.3. The computational efficiency of the method is
analyzed for different problems in Sec. 3.4, based on numerical results.
3.2
Bisection Method
The bisection method is used to find the root of a nonlinear real-valued scalar function
f (x) = 0 (f : R → R). The method is initialized with the limits of an interval
[xa , xb ], where the function is defined (see Fig. 3.1). The signs of f (xa ) and f (xb )
must be opposite. In this case xa and xb are bracketing at least one root, since,
by the Intermediate Value Theorem, the function f must have at least one zero in
the interval (xa , xb ). At the beginning of the iteration, the midpoint of the interval
xc = (xa + xb )/2 is computed and f (xc ) is evaluated. If the sign of f (xc ) is the same
as the sign of f (xa ), then xc is set as a new value instead of xa , otherwise, xc is set as
a new value instead of xb and the iteration is repeated. If xc happens to be a root of
f (i.e., f (xc ) = 0) than the iteration is completed in finite number of steps. In this
process, the length of the interval is reduced by 50% in each step, leading to a strictly
monotone linear convergence. The drawback of this process is that it can find only
one intersection point. In the next subsection, a generalized version of this method
is described, which is able to find numerous roots.
86
}
}
}}
}}
}
x
}
}
}
}
}
}
}}
unrecognized
roots
initial mesh
}
f(x)
refined intervals
Figure 3.2: Numerous roots computed by bisection method, evaluated on an initial mesh.
3.2.1
Numerous Nearby Roots
The original form of the bisection method can easily be extended to find numerous
roots of a non-linear equation in a given interval. If the function values are computed
in an initial mesh on the examined interval, then the original method can be used for
each neighboring points where the signs of the function values are different. This way,
some roots may be skipped if the initial mesh is not fine enough and even number of
roots are placed inside a single interval. An example is shown in Fig. 3.2.
Note, that the number of the "useful" intervals that are bracketing the roots, is
the same as the number of the detected roots, so if we use a relatively fine initial
mesh, then the detection of all roots can be ensured without the dramatic increase
of the computational time. It is also important that the bracketing intervals are
mapped into themselves after one iteration step. In case of Newton’s method, e.g.,
the monotone convergence cannot always be guaranteed.
3.2.2
Generalization for Higher Dimensions
In many applications, the roots of a system of non-linear equations f (x) = 0 must
be computed. Let f : RDS → RDC with RDS representing the DS -dimensional vector
space and DC is the codimension (or relative dimension) of the subspace or submanifold of the solution. The topological dimension of the submanifold is DT = DS − DC .
For example: a surface (DT = 2) can be defined by one equation DC = 1 in a 3D space
(DS = 3). Note, that in case of fractal type submanifolds, the topological dimension
DT and the fractal dimension DF of the object are not the same (DT ≤ DF ≤ DS ).
The fractal dimension has also been characterized as a measure of the space-filling
capacity of a pattern that tells how a fractal scales differently than the space it is embedded in; a fractal dimension does not have to be an integer [105, 106]. For example:
Julia set for c = −1 is defined by the periodic points of map g(x) = (x1 + ix2 )2 − 1
in the complex plane resulting DT = 1 and DF = 1.2683 (DS = 2, DC = 1) [107].
Solution methods already exist and are well developed for the most typical cases.
In case DF = DT = 0, the Newton’s method can easily be generalized for higher dimensions by using the gradient of f . Numerical solutions can be found by computer algebraic software for DC = 1 and DS = 2, 3 (see Wolfram Mathematica: ContourP lot
and ContourP lot3D, Matlab: contour and isosurf ace, etc.) [108] [109]. For large
values of DS , the continuation methods (like [110] [111] [112]) are used usually. These
are typically applied for DF = DT = 1, however, multiparameter continuation meth87
ods also exist [113]. The main drawbacks of these continuation methods are that the
proper initial value is needed and the closed submanifolds (compact submanifolds or
isolas) can hardly be detected. From the engineering point of view, it is important
to find all the roots in a given domain automatically, which could be a set of closed
submanifolds. Our main goal is to create a robust algorithm based on the bisection
method, which can determine a properly dense point cloud of all the submanifolds
with a given accuracy.
In the following subsections, the details of the suggested method and its main
steps are described and analyzed.
3.2.3
n-simplexes and n-cubes
The generalized version of the initial interval of the bisection method must be a higher
dimensional object, which can be divided into multiple self-similar objects. The first
obvious choice could be the n-simplex as presented in [114], which is a generalization
of a triangle or a tetrahedron to higher dimensions. The second obvious choice is the
n-cube (or hypercube), which is a generalization of a square or a cube. In a numerical
implementation, the n-cube has many advantages. The data corresponding to a node
(or vertex) of an n-cube can easily be stored in multilevel arrays (or hypermatrices),
the coordinates of the neighbours of an n-cube can easily be determined. Furthermore,
the halving of an n-cube along each dimension is trivial, while halving of an n-simplex
is complicated in higher dimensions.
3.2.4
Selection of Bracketing n-cubes
In higher dimensions, an initial mesh must be defined to ensure the detection of
multiple submanifolds similarly as it was shown in Sec. 3.2.1. First, two limit values
xaj and xbj and the size of the mesh Nj must be defined for each dimension (j ∈
{1, 2, ...DS }). The edge length of the initialized n-cube in the j th dimension is ∆xj =
(xbj − xaj )/(Nj − 1). In the next step, the function value must be computed in all the
QDS
j=1 Nj nodes. Before the refinement, the bracketing n-cubes which contain parts of
the submanifold must be selected from all the initialized n-cubes. This step also has
to be modified in case DC > 1 according to the given problem.
3.2.4.1
Safe Selection
If the signs of fk are different at any two nodes of an n-cube for all k = 1, 2, ...DC ,
then it is possible, that the n-cube brackets a part of a submanifold. Unfortunately,
this is just a necessary condition, however, this condition can be used to define an
n-cube as a bracketing one. In this case, non-bracketing n-cubes will also be selected
and refined, which does not lead to computational error if further iteration is applied,
because for smaller refined n-cubes, a ghost bracketing will disappear. The only minor
problem is that the unnecessary selections will increase the computational time.
88
3.2.4.2
Secant Based Selection
A more precise selection of the bracketing n-cubes is based on the secant method. In
this case, the submanifold inside the n-cube is to be found by linear interpolations.
The iteration starts with connecting all the nodes with positive f1 values to all the
nodes with negative f1 values, then the roots of f1 along these lines are linearly approximated by the secant method. Then,f2 , f3 , . . . , fDC are also linearly interpolated
at these points.
The above process is then repeated for f2 , then for f3 , ... and fDC . If in any stage
of the iteration, the signs of fk are the same in all interpolated points, then the next
step of the iteration can not be carried out, which means that the current n-cube is
non-bracketing. If an exact 0 is found, then it is a bracketing n-cube.
3.2.4.3
Fitting Hyper-Planes
Another possible solution is to fit a hyper-plane to the function values of the nodes at
the current n-cube. The normal vectors of the hyper-plane can be calculated as well
as the hyper-plane point closest to the center of the n-cube. If this point is "close
enough" to the n-cube then the n-cube can be selected as a bracketing one.
3.2.4.4
Comments on Selection Methods
An advantage of the selection methods presented in the above subsection is that the
additional evaluation of f is not needed. Other types of selection methods could be
created if we compute f at additional points.
The appropriate method of the selection depends on the given problem. The
safe selection is very fast, but the computational time increases if the codimension-1
surfaces (fj (x) = 0) are almost parallel at the intersection points. Still, it is advised
to use this if the evaluation of f is fast enough. The selection of the bracketing n-cubes
can be the bottleneck of the computation.
If the evaluation of f is slow, then the selection based on the secant method or
the fitting of hyper-plane provide better performance. Based on our experiments, if
DS < 4 then the secant based selection is faster than fitting a hyper-plane. These
two methods give poor results if the curvatures of the surfaces are relatively large
compared to the size of the n-cube.
If DT ≥ 1, the missing parts of the submanifold can be found by using some kind
of continuation method at the very end of the iteration.
3.2.5
Continuation at Neighboring n-cubes
At the end of the iterations, all the neighbors of the bracketing n-cubes must be
examined. If additional bracketing n-cubes are found then a missing part of the submanifold is detected. These tests should be repeated till all the neighboring n-cubes
of the bracketing n-cubes are non-bracketing. In this case, the whole submanifold is
found in the selected region, and the iteration is closed. This method is similar to a
multiparameter continuation method, because if a small part of the submanifold is
89
1
initail
mesh
0.8
0.6
0.4
x2
0.2
0
-0.2
-0.4
0.1
0.05
-0.6
0
-0.8
-0.05
0.15
-1
-1
-0.8
-0.6
-0.4
-0.2
0
x1
0.2
0.2
0.4
0.25
0.3
0.6
0.35
0.8
0.4
1
Figure 3.3: Numerical result with continuation algorithm. The black crosses denote the points where
f is evaluated, the dense bright gray dots show up as a line and denote the roots of f computed by
the secant method in the final bracketing n-cubes. Parameters are DS = 2, DC = 1, DT = DF = 1,
xa1 = xa2 = −1, xb1 = xb2 = 1.1, N1 = N2 = 3.
detected, then it is extended during the iterations until it is closed into itself (in case
of an isola) or until the boundary of the selected range is reached.
This continuation method is efficient, because the evaluation of f is carried out
only along the submanifold.
An example is shown in Fig. 3.3, where a unit circle in 1/3 metric f (x) = kxk1/3 −1
1/p
P
p
DS
|x
|
is evaluated in a plane. We can clearly see, that the thin
with kxkp =
i
i=1
parts of the object was not detected by the bracketing n-cubes, but even these parts
are discovered efficiently by the continuation of the bracketing n-cubes.
The final step of the method is the interpolation of the roots inside the bracketing
n-cubes, which can be done by the secant method, the hyper-plane fitting or other
advanced methods.
All the necessary steps of the proposed algorithm are described and collected in
proper sequence on the flowchart in Fig. 3.4.
3.3
Numerical Implementation
The proposed method is implemented and tested in Matlab. To create a fast algorithm, indexing of the n-cubes and the nodes should be optimized. It should be stored
which n-cubes are bracketing or non-bracketing and which nodes had already been
computed, in order to avoid the redundant computation and analysis of a node or an
n-cube.
The memory problems can be eliminated in case of large DS if spares matrices or
linearly indexed arrays are used to store the values of f which are typically computed
90
91
Input:
f(x)=0
xja,xjb,Nj j=1,...,DS
Number of iteration: I
Create the initail mesh
and define the n-cubes
Evaluate f(x) at the new
nodes of the new n-cubes
Select the bracketing cubes
End of iteration?
Yes
No Refine the
bracketing
n-cubes
Select the neighbouring cubes
Evaluate f(x) in the new
nodes of the new n-cubes
Are
Yes
there new bracketing
cubes?
No
Interpolation of the roots in the
bracketing n-cubes
Figure 3.4: Flowchart of the proposed Multi-Dimensional Bisection Method.
Figure 3.5: Application of the Multi-Dimensional Bisection Method: point cloud and the fitted
surface of the stability boundary of the shimmy problem [115]. Gray scale map refers to the 4th
dimension ω.
very sparsely in the selected intervals. Still, memory problems could occur due to the
large number of computed points if the fractal dimension of the object DF is large
and many iterations are used.
The result of the proposed method is a dense point cloud. In some cases, it
is necessary to connect the neighboring points and determine the edges and/or the
surfaces of the detected submanifold (e.g., for plotting). This additional computation
can be done by means of a simplex-based algorithm, the details of which are out of
the focus of this study.
In order to demonstrate the efficiency of the algorithm, we present the computation
of the point cloud that corresponds to the stability boundary surface of a shimmying
wheel, which is described and analyzed in details in [115]. The surface was determined
by a simplex based method. Equation (4.5) in [115] defines the stability boundary
surface of a shimmy problem for a 4-parameter space (DS = 4, DC = 2, DT = DF =
2). Our final result in Fig. 3.5 fits perfectly to Fig. 4.5 of the dissertation [115], which
was calculated by a basically different semi-analytical method.
3.4
Efficiency Number
In order to analyze the efficiency of the proposed method, the number NP of the points
is determined where the function f is evaluated. The efficiency of the computation
time can be measured similarly to the definition of the Hausdorff (or fractal) dimension
[116] [117]. This dimension is closely related to the box-counting dimension D0 , which
considers the number of boxes NB of size ε that contains the submanifold (attractor)
92
and is given by
log NB (ε)
.
ε→0 log(1/ε)
D0 = lim
(3.3)
This definition can be applied directly to our method where n-cubes are refined.
Thus, define
log Nx (I)
Dx = lim
,
(3.4)
I→∞
log 2I
where I is the number of the iteration steps of the refinements of the n-cubes. If
we substitute Nx with the number of bracketing cubes then the fractal dimension
of the submanifold DF is determined, but if the number of function evaluation NP
is substituted, then we get the fractal dimension of the evaluated points DP . This
number shows, how dense the evaluated points are in the space. We define the
efficiency number E of the computation method as follows
1
.
(3.5)
DP − DF
First, look at the so-called "brute force method", where the function is evaluated
at all points in the final mesh for which NP is given by
E=
NP =
DS
Y
I
I DS
(Nj − 1)2 + 1 ≈ 2
j=1
DS
Y
Nj .
(3.6)
j=1
It is clear, that for the brute force method, the fractal dimension of the evaluated
points DP tends rapidly to the dimension of the whole vector space DS , which leads
to a small value of the efficiency number: E = 1/(DS − DF ). Note, that E = 1/DC in
case of non-fractal type manifolds (recall that DC is the codimension of the subspace
or submanifold of the solution, see Sec. 3.2.2).
The limit of expression (3.4) cannot be determined by a numerical method, however the tendency can be analyzed by increasing the number I of iterations. In
the next example, we determine the surface of a unit sphere in Euclidean norm by
f (x) = kxk2 − 1 in a volume defined by xaj = −2, xbj = 2, Nj = 3 and j ∈ {1, 2, 3}
(DS = 3, DC = 1, DT = DF = 2).
It can be seen in Fig. 3.6, that the fractal dimensions converge almost exponentially to a given value, which can be approximated by the limit of the fitted
exponential curve. The efficiency number of the afore mentioned brute force method
tends to 0.9989, which is very close to the correct value E = 1. The efficiency
number of the proposed Multi-Dimensional Bisection Method tends to E = 3.7821
(DP = 2.2644). An ideal perfect method, which computes the dense pointcloud only
on the submanifold would have DP = DF and E → ∞.
The examined unit sphere does not have fractal shape, so the fractal and the
topological dimensions must be the same (DT = DF = 2). The fitted continuous
curve in Fig. 3.6 tends to DF = 2.0328 which is very close to the topological dimension
DT = 2.
93
3.5
DF ,DP
3
2.5
2
1.5
3
4
5
6
7
I - number of iteration
8
Figure 3.6: Convergence of the fractal dimensions as a number of iterations. Dots denote the
evaluated box-counting dimension for the integer values of I. Lines denote the fitted exponential
curve in a form D = A + Be−λI . Dashed line denotes the fractal dimension of the submanifold
DS . The fractal dimension DP of the evaluated points for the brute force method and for the
proposed Multi-Dimensional Bisection Method are denoted by the dotted and the continuous lines,
respectively.
0.8
0.6
0.4
x2
0.2
0
-0.2
-0.4
-0.6
-0.8
-2
-1.5
-1
-0.5
0
x1
0.5
1
1.5
2
Figure 3.7: Julia set computed by Multi-Dimensional Bisection Method for c = −1.
The same analysis was performed for the afore mentioned Julia set [107]) with an
initial mesh xaj = −2, xbj = 2, Nj = 3 and j ∈ {1, 2} (max(I) = 12), see Fig. 3.7.
During the numerical evaluation, the fractal dimension of the bracketing n-cubes
tends to DF = 1.2753, which has only ≈ 0.5% error relative to the fractal dimension
of the Julia set. Thus, as a side result, the proposed Multi-Dimensional Bisection
Method can be also used to determine the fractal dimension of an object.
The efficiency number was also computed for other non-fractal type problems and
the results are collected in Table 3.1. The dimension of the vectorspace DS and the
number of codimensions DC were varied from 1 to 4. The initial mesh was defined
by xaj = −2, xbj = 2, Nj = 3 and j ∈ {1, 2, ..., DS }. The corresponding objects are
defined by the equations
f1 (x) = kxk2 − 1,
fk (x) = xk − xk−1
for k > 1 and k ∈ {1, 2, . . . , DC },
(3.7)
(3.8)
which define a unit n-sphere in Euclidean space and codimension-1 hyper-planes,
94
Table 3.1: The efficiency number of the Multi-Dimensional Bisection Method for Eqs. (3.7, 3.8).
E
DS =1
DS =2
DS =3
DS =4
Brute Force
DC =1
2.9878
5.4795
3.7821
7.7821
1
DC =2 DC =3
×
×
2.0500
×
1.6603 1.1860
1.3978 0.9936
0.5
0.333
DC =4
×
×
×
0.8103
0.25
axial depth of cut (mm)
100
ZOA
UNSTABLE
80
MFS
Instability
gradient
2
EMFS
SD
60
Number of unstable
Floquet multipliers
(EMFS)
40
1
2
20
1
STABLE
1
0
0
0
2000
4000
6000
8000
10000
spindle speed (rpm)
Figure 3.8: Application of the Multi-Dimensional Bisection Method: Stability chart of milling with
varying helix pitch [104]. The chart is computed by means of the Zeroth Order Approximation (ZOA,
black thin dashed lines), the traditional Multi Frequency Solution (MFS, thick continuous lines),
the Extended Multi Frequency Solution (EMFS, fuzzy lines) and the Semi-Discretization method
(SD, shaded area). Instability gradient vectors are also visualized in case of EMFS. Note, that these
vectors overlap at low spindle speeds due to the scaling.
respectively.
Table 3.1 shows that the computational efficiency varies depending on the given
problem. For the ideal method, these numbers would tend to infinity. In all cases,
the efficiency numbers are much better than the ones of the brute force method. A
linear surface can be fitted well on the computed values of the fractal dimension of
the evaluated points DP = 0.9764DS − 0.6692DC ≈ Ds − 2/3DC , which leads to
E ≈ 3/DC ; this shows that the efficiency is about 3 times larger than that of the
brute force method.
The developed Multi-Dimensional Bisection Method is especially efficient and useful when intricate stability charts are to be determined with the extended multi
frequency solution in cases of complex milling processes where fractal-like stability
boundaries and/or stable and unstable islands occur like the ones in Fig. 3.8 [104].
95
New Results
Thesis 7
In industrial applications of chatter predictions for cutting, the fast, possibly online
reconstruction of the stability charts is essential even in cases of complex real-world
machine tool dynamics. To achieve this goal, the bisection method is generalized to
higher dimensions for stability chart reconstruction. A fast numerical algorithm of
this ’Multi-Dimensional Bisection Method’ is built in Matlab, which is able to find the
submanifolds formed by the sets of the roots of a system of nonlinear equations, where
the number of equations, that is the codimension DC of the submanifolds, is smaller
or equal than the number of the independent variables, which gives the dimension
DS of the parameter space. The developed method can be used for the automatic
determination of disjunct solution sets in selected regions for optional dimensions and
codimensions.
An efficiency number is introduced to characterize the performance of the corresponding numerical methods. It is based on the box-counting dimension DP of the
evaluated points in the space of the parameters, and it is defined as
1
,
DP − DF
where DF is the fractal dimension of the submanifolds. It is shown that the efficiency
number of the so-called brute force method, that checks a uniform net of the selected
region of the space (DP = DS ), is the small value E = 1/(DS − DF ), which is
E = 1/DC in case of non-fractal type submanifolds. An ideal perfect method, which
computes the dense pointcloud on the submanifold only, would have DP = DF and
E → ∞.
Based on numerical tests, the efficiency number of the proposed ’Multi-Dimensional
Bisection Method’ is approximated by E ≈ 3/DC . This means that the developed
numerical method is efficient enough to present stability charts online for milling processes in the DS = 3 dimensional parameter space of spindle speed, axial depth of cut
and chatter frequency and for the DF = 1 dimensional instability lobes, even when
the lobes have fractal-like structure and/or they form closed stability or instability
islands.
E=
Related publications:
[118, 104, 119]
96
Bibliography
[1] Schuman, W., 2007. “Das bearbeiten harter eisenwerkstoffe mit geometrisch
bestimmter schneide”. Industrie Diamanten Rundschau, III(07), pp. 52–71.
[2] Mészáros, I., and Szepesi, D., 2006. “Optimization of finish hard turning”. In
EUSPEN 6th International Conference, Vienna, pp. 100–103.
[3] Takács, M., and Verõ, B., 2003. “Material structural aspects of micro scaled
chip removal”. Materials Science Forum, 414 - 415, pp. 337–342.
[4] Meszáros, I., and Szepesi, D., 2005. “Edzett acélok nagypontosságú esztergálása
i-iv.”. Gépgyártás, XLV.
[5] Kratz, H., 2006. “Belastungsoptimierte werkzeuge in wichtigen anwendungsgebieten der cbn schneidstoffe”. Industrie Diamanten Rundschau, III/06,
pp. 62–67.
[6] Chou, K., and Evans, C., 1999. “Cubic boron nitride tool wear in interrupted
hard cutting”. Wear, 225-229, pp. 234–245.
[7] Ewins, D. J., 2000. Modal Testing: theory, practice and application. Research
Studies Press. John Wiley and Sons Inc., England.
[8] Dombovari, Z., and Stépán, G., 2012. “Marószerszámok dinamikai tulajdonságai
és azok hatása a megmunkálás stabilitására”. GÉP, 65(12), pp. 163–174.
[9] Stépán, G., 1989. Retarded Dynamical Systems. Longman, Harlow.
[10] Insperger, T., and Stépán, G., 2011. Semi-Discretization for Time-Delay Systems, 1 ed. Applied Mathematical Sciences. Springer.
[11] Insperger, T., Stepan, G., and Turi, J., 2008. “On the higher-order semidiscretizations for periodic delayed systems”. Journal of Sound and Vibration,
313, pp. 334–341.
[12] Insperger, T., and Stépán, G., 2004. “Updated semi-discretization method for
periodic delay-differential equations with discrete delay”. International Journal
of Numerical Methods in Engineering, 61(1), pp. 117–141.
97
[13] Insperger, T., and Stépán, G., 2001. “Semi-discretization of delayed dynamical
systems”. In Proceedings of the ASME 2001 Design Engineering Technical
Conferences, no. DETC2001/VIB-21446 (CD-ROM), ASME, p. 6.
[14] Szalai, R., and Stepan, G., 2006. “Lobes and lenses in the stability chart of interrupted turning”. Journal of Computational and Nonlinear Dynamics, 1(3),
pp. 205–212.
[15] Bachrathy, D., and Meszaros, I., 2009. “Dynamical problems in interrupted high
precision hard turning”. In LAMDAMAP 2009: 9th International Conference
and Exhibition on Laser metrology, machine tool, CMM and robotic performance. London, England, 06.30-07.02, no. ISBN: 978-0-9553082-7-7, pp. 357–
367.
[16] Bachrathy, D., and Meszaros, I., 2009. “Dynamical analysis of high precision
hard turning processes for interrupted machining (nagypontossagu kemeny esztergalas dinamikai vizsgalata megszakitott feluletek eseten) (in: Hungarian)”.
Gepgyartas, XLIX(4-5), pp. 47–53.
[17] Bachrathy, D., and Meszaros, I., 2009. “Nagypontosságú kemény esztergálás dinamikai vizsgálata megszakított felületek esetén”. Gépgyártás, XLIX, pp. 47–
53.
[18] Lim, J., Sohn, Y., and Kwon, D., 2007. “Real-time accelerometer signal processing of end point detection and feature extraction for motion detection”. In 10th
IFAC/ IFIP/ IFORS/ IEA Symposium on Analysis, Design, and Evaluation of
Human-Machine Systems.
[19] Yang, J., Chang, W., Bang, W., Choi, E., Kang, K., Cho, S., and Kim, D.,
2004. “Analysis and compensation of errors in the input device based on inertial
sensors”. In Int. Conf. on Information Technology: Coding and Computing
(ITCC’04), Vol. 2.
[20] Y. Huang, S. Y. L., 2003. “Cutting forces modeling considering the effect of
tool thermal property-application to cbn hard turning”. International Journal
of Machine Tools and Manufacture, 43.
[21] Meszaros, I., Bachrathy, D., and Farkas, B., 2008. “Dynamical problems in
high precision hard cutting”. In Biannual 19th International Conference on
Manufacturing 2008. Budapest, Hungary, 11.06-11.07, no. ISBN: 978-963-905824-8, pp. 18–24.
[22] Bachrathy, D., and Meszaros, I., 2010. “Optimal cutting speeds, entrance and
exit force in interrupted high precision hard turning”. In CIRP ICME Š10 7th CIRP International Conference: Intelligent Computation in Manufacturing Engineering. Capri, Italy, 06.23-06.25, Vol. B4/3, ISBN: 978-88-95028-65-1,
pp. 1–4.
98
[23] Bachrathy, D., Reith, M. J., and Meszaros, I., 2010. “Optimal cutting speeds
and surface prediction in interrupted high precision hard turning”. In Gepeszet
2010, Proceedings of the Seventh Conference on Mechanical Engineering. Budapest, Hungary, 05.25-05.26, no. 036, ISBN: 978-963-313-007-0. Budapest University of Technology and Economics, Budapest University of Technology and
Economics, pp. 241–246.
[24] Bachrathy, D., and Meszaros, I., 2010. “Surface modeling, optimal cutting
speeds and entrance and exit force in interrupted high precision hard turning”. In 4th CIRP International Conference on High Performance Cutting.
Gifu, Japan, 10.24-10.26, no. C16, pp. 1–6.
[25] Reith, M. J., Bachrathy, D., and Meszaros, I., 2010. “Smoothed force model for
interrupted high precision hard turning”. In Manufacturing 2010: The XX. Conference of GTE on Manufacturing and related technologies. Budapest, Hungary,
10.20-10.21, no. 17, ISBN: 978-963-9058-31-6, pp. 1–8.
[26] Tobias, S. A., 1965. Machine Tool Vibration. Blackie and Son, Ltd., London.
[27] Insperger, T., Gadisek, J., Kalveram, M., G., S., Weinert, K., and Goverkar,
E., 2004. “Machine tool chatter and surface quality in milling processes”. In
Proceedings of ASME International Mechanical Engineering Conference and
Exposition, Anaheim CA, paper no. IMECE2004-59692 (CD-ROM).
[28] Montgomery, D., and Altintas, Y., 1991. “Mechanism of cutting force and
surface generation in dynamic milling”. Transactions of ASME Journal of Engineering for Industry, 113, pp. 160–168.
[29] Shirase, K., and Altintas, Y., 1991. “Cutting force and dimensional surface
error generation in peripheral milling with variable pitch helical end mills”.
International Journal of Machine Tools and Manufacture, 113, pp. 160–168.
[30] Mann, B. P., Young, K. A., Schmitz, T. L., and Diley, D. N., 2005. “Simultaneous stability and surface location error predictions in milling”. Journal of
Manufacturing Science and Engineering, 127, pp. 446–453.
[31] Tlusty, J., and Spacek, L., 1954. Self-excited vibrations on machine tools. Nakl.
CSAV, Prague. in Czech.
[32] Tobias, S., and Fishwick, W., 1958. “The chatter of lathe tools under orthogonal
cutting conditions”. Transactions of the ASME, 80, pp. 1079–1088.
[33] Tlusty, J., and Poláček, M., 1963. “The stability of machine-tool against selfexcited vibration in machining”. Proceedings of the International Research in
Production Engineering, American Society of Mechanical Engineers (ASME),
p. 465.
[34] Szalai, R., Stépán, G., and Hogan, S., 2004. “Global dynamics of low immersion
high-speed milling”. Chaos, 14(4), pp. 1069–1077.
99
[35] Stepan, G., Szalai, R., Mann, B. P., Bayly, P. V., Insperger, T., Gradisek, J.,
and Govekar, E., 2005. “Nonlinear dynamics of high-speed milling - analysis,
numerics, and experiments”. ASME Journal of Vibration and Acoustics, 127,
pp. 197–203.
[36] Dombóvári, Z., Wilson, R., and Stépán, G., 2008. “Estimates of the bistable
region in metal cutting”. Proceedings of the Royal Society A: Mathematical,
Physical and Engineering Science, 464(2100), pp. 3255–3271.
[37] Insperger, T., Stépán, G., Hartung, F., and Turi, J., 2005. “State dependent
regenerative delay in milling processes”. In Proceedings of ASME International
Design Engineering Technical Conferences, no. DETC2005-85282 (CD-ROM),
ASME, p. 10.
[38] Gradisek, J., Govekar, E., Grabec, I., Kalveram, M., Weinert, K., Insperger,
T., and Stepan, G., 2005. “On stability prediction for low radial immersion
milling”. Machining Science and Technology, 9(1), pp. 117–130.
[39] Zatarain, M., Muñoa, J., Peigné, G., and Insperger, T., 2006. “Analysis of the
influence of mill helix angle on chatter stability”. CIRP Annals - Manufacturing
Technology, 55(1), pp. 365–368.
[40] Das, M., and Tobias, S., 1967. “The relation between the static and the dynamic
cutting of metals”. International Journal of Machine Tool Design and Research,
7, pp. 63–89.
[41] Altintas, Y., Eynian, M., and Onozuka, H., 2008. “Identification of dynamic
cutting force coefficients and chatter stability with process damping”. CIRP
Annals - Manufacturing Technology, 57(1), pp. 371–374.
[42] Budak, E., and Tunc, L. T., 2009. “A new method for identification and modeling of process damping in machining”. Journal of Manufacturing Science and
Engineering, 131(5), p. 10.
[43] Budak, E., and Tunc, L. T., 2010. “Identification and modeling of process damping in turning and milling using a new approach”. CIRP Annals Manufacturing
Technology, 59(1), p. 6.
[44] Eynian, M., and Altintas, Y., 2009. “Chatter stability of general turning operations with process damping”. J. Manuf. Sci. Eng., 131(4), pp. 1–10.
[45] Wu, D. W., 1988. “Application of a compreshensive dynamic cutting force model
to orthogonal wave-generating processes”. Mechanical Sciences, 30(8), p. 581.
[46] Chiou, Y. S., Chung, E. S., and Liang, S. Y., 1995. “Analysis of tool wear effect
on chatter stability in turning”. International Journal of Mechanical Sciences,
37(4), pp. 391–404.
100
[47] Budak, E., Altintas, Y., and Armarego, E. J. A., 1996. “Prediction of milling
force coefficients from orthogonal cutting data”. J. Manuf. Sci. Eng, 118,
pp. 216–225.
[48] Insperger, T., Mann, B. P., Stépán, G., and Bayly, P. V., 2003. “Stability of upmilling and down-milling, part 1: alternative analytical methods”. International
Journal of Machine Tools and Manufacture, 43(1), pp. 25 – 34.
[49] Insperger, T., Gadisek, J., Kalveram, M., G., S., Weinert, K., and Goverkar,
E., 2006. “Machine tool chatter and surface location error in milling processes”.
Journal of Manufacturing Science and Engineering, 128, pp. 913–920.
[50] Insperger, T., and Stépán, G., 2002. “Semi-discretization method for delayed
systems”. International Journal for Numerical Methods in Engineering, 55,
pp. 503–518.
[51] Kline, W., 1982. “The prediction ofcutting forces and surface accuracy for
the end milling process”. Dissertation, Ph.D. Thesis, University of Illinois at
Urbana-Champaign, 223 pages; AAT 8209596.
[52] Peigné, G., Paris, H., and Brissaud, D., 2004. “Surface shape prediction in high
speed milling”. International Journal of Machine Tool and Manufacture, 44,
p. 1567Ű1576.
[53] Surmann, T., 2006. “Geometric model of the surface structure resulting from
the dynamic milling process”. In Proceedings of the 9th CIRP International
Conference on Modeling of Machining Operations, no. ISBN 961-6536-06-0,
CIRP, pp. 187–192.
[54] Xu, A.-P., Qu, Y.-X. a. Z. D.-W., and Huang, T., 2003. “Simulation and
experimental investigation of the end milling process considering the cutter
flexibility”. International Journal of Machine Tools and Manufacture, 43,
p. 283Ű292.
[55] Schmitz, T., and Ziegert, J., 1999. “Examination of surface location error due
to phasing of cutter vibrations”. Precision Engineering, 23(1), pp. 51 – 62.
[56] Gradišek, J., Kalveram, M., Insperger, T., Weinert, K., Stépán, G., Govekar,
E., and Grabec, I., 2005. “Stability prediction for milling”. International Journal
of Machine Tools and Manufacture, 45(7-8), pp. 769–781.
[57] Gradišek, J., Kalveram, M., and Weinert, K., 2004. “Mechanistic identification
of specific force coefficients for a general end mill”. International Journal of
Machine Tools and Manufacture, 44(5), pp. 401–414.
[58] Randall, R. B., 1987. Frequency Analysis, 3rd edition 1st print ed. Brüel &
Kjaer: Naerum, Denmark.
101
[59] Schmitz, T., and Smith, K. S., 2009. Machining Dynamics: Frequency Response
to Improved Productivity. Springer, New York, NY.
[60] Munoa, J., 2007. “Desarrollo de un modelo general para la predicción de la
estabilidad del proceso de fresado. aplicación al fresado periférico, al planeado
convencional y a la caracterización de la estabilidad dinámica de fresadoras
universales”. Mondragon University, PhD Thesis (in Spanish).
[61] Bachrathy, D., Insperger, T., and Stepan, G., 2009. “Surface properties of
the machined workpiece for helical mills”. Machining Science and Technology,
13(2), pp. 227–245.
[62] Bachrathy, D., Homer, M., Insperger, T. I., and Stepan, G., 2007. “Surface
location error for helical mills”. In 6th International Conference on High Speed
Machining. San Sebastian, Spain, 03.21-03.22, no. Paper C100., pp. 379–384.
[63] Bachrathy, D., and Stepan, G., 2007. “Surface error for helical mills”. In 6th
International Congress on Industrial and Applied Mathematics, Zurich, Switzerland, 07.16-07.20, no. Paper 1268, p. 131.
[64] Bachrathy, D., Insperger, T., and Stepan, G., 2007. “Computation of surface
quality in case of helical milling tools (felületi minõség számítása csavart élû
szerszámmal történõ marás során) (in hungarian)”. In X. Magyar Mechanikai
Konferencia, Miskolc, Hungary, 08.27-08.29, pp. 1–4.
[65] Bachrathy, D., and Stepan, G., 2010. “Optimal axial immersion for helical
milling tools based on frequency response function (optimalis axialis fogasmelyseg csavart elu maroszerszamra frekvencia atviteli fuggveny alkalmazasaval (in
hungarian)”. Gep, LXI(9-10), pp. 3–6.
[66] Bachrathy, D., and Stepan, G., 2007. “Good surface properties at efficient
technological parameters in milling process”. In 24th Danubia-Adria: Symposium on Developments in Experimental Mechanics. Sibiu, Romania, 09.19-09.22,
no. ISBN: 978-973-739-456-9, pp. 45–46.
[67] Bachrathy, D., and Stepan, G., 2008. “Experimental setup for fast stability
chart reconstruction of milling processes”. In Proceedings of Sixth Conference
on Mechanical Engineering. Budapest, Hungary, 05.29-05.30, no. G-2008-H-17,
ISBN: 978-963-420-947-8, Budapest University of Technology and Economics,
pp. 1–7.
[68] Bachrathy, D., and Stepan, G., 2008. “Efficient experimental detection of milling
stability boundary and the optimal axial immersion for helical mills”. In International Multi-Conference on Engineering and Technological Innovation: IMETI
2008: International Symposium on Manufacturing Systems and Technologies:
ISMST 2008. Orlando, USA, 06.29-07.02, Vol. 1, pp. 7–11.
102
[69] Richardson, H. M., and Formenti, D. L., 1982. “Parameter estimation from
frequency response measurements using rational fraction polynomials”. In 1st
IMAC Conference, Orlando, FL.
[70] Schmitz, T. L., and Duncan, G. S., 2005. “Three-component receptance coupling substructure analysis for tool point dynamics prediction”. Journal of
Manufacturing Science and Engineering, 127(4), November, pp. 781–790.
[71] Okwudire, C. E., and Altintas, Y., 2009. “Hybrid modeling of ball screw drives
with coupled axial, torsional, and lateral dynamics”. Journal of Mechanical
Design, 131(7), p. 071002.
[72] Engin, S., and Altintas, Y., 2001. “Mechanics and dynamics of general milling
cutters.: Part i: helical end mills”. International Journal of Machine Tools and
Manufacture, 41(15), pp. 2195 – 2212.
[73] William Weaver, Stephen Timoshenko, D. H. Y., 1990. Vibration problems in
engineering, 5 ed. John Wiley & Sons Inc.
[74] Davis, R., Henshell, R., and Warburton, G., 1972. “A timoshenko beam element”. Journal of Sound and Vibration, 22(4), pp. 475 – 487.
[75] Yokoyama, T., 1990. “Vibrations of a hanging timoshenko beam under gravity”.
Journal of Sound and Vibration, 141(2), pp. 245 – 258.
[76] Schmidt, C., 2006. “Level curve based geometry representation for cutting force
prediction”. In DAAAM International Symposium.
[77] Denkena, B., Kruger, M., Daniel, B., and Gabor, S., 2012. “Model based reconstruction of milled surface topography from measured cutting forces”. International Journal of Machine Tools and Manufacture, 54-55, pp. 25–33.
[78] Budak, E., and Altintas, Y., 1996. “Prediction of milling force coefficients from
orthogonal cutting data”. Journal of Manufacturing Science and Engineering,
118(2), pp. 216–225.
[79] Engelborghs, K., Luzyanina, T., and Samaey, G., 2001. DDE-BIFTOOL v.
2.00: a Matlab package for bifurcation analysis of delay differential equations.
Tech. Rep. TW330, Department of Computer Science, K.U. Leuven, Leuven,
Belgium.
[80] Insperger, T., and Stepan, G., 2011. Semi-Discretization for Time-Delay Systems: Stability and Engineering Applications. Springer, New York.
[81] Taylor, F. W., 1907. On the art of cutting metals. Transactions of the American
Society of Mechanical Engineers.
[82] Shi, H. M., and A., S., 1984. “Theory of finite amplitude machine tool instability”. International Journal of Machine Tool Design and Research, 24(1),
pp. 45–69.
103
[83] Faassen, R., van de Wouw, N., Oosterling, H., and Nijmeijer, H., 2005. “Updated
tool path modelling with periodic delay for chatter predition in milling”. In 5th
EUROMECH Nonlinear Oscillations Conference (ENOC), D. van Campen, ed.,
no. CDROM.
[84] Insperger, T., and Stépán, G., 2000. “Stability of high-speed milling”. In Proceedings of ASME International Mechanical Engineering Congress and Exposition, no. 241, ASME, pp. 119–123.
[85] Laczik, B., 1986. “Vibration monitoring of cutting processes”. PhD thesis, BME.
[86] Kerschen, G., Peeters, M., Golinval, J., and Vakakis, A., 2009. “Nonlinear
normal modes, part i: A useful framework for the structural dynamicist”. Mechanical Systems and Signal Processing, 23(1), January, pp. 170–194.
[87] Peeters, M., Viguié, R., Sérandour, G., Kerschen, G., and Golinval, J.-C., 2009.
“Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques”. Mechanical Systems and Signal Processing,
23(1), January, pp. 195–216.
[88] Hartung, F., 2005. “Linearized stability in periodic functional differential equations with state-dependent delays”. Journal of Computational and Applied
Mathematics, 174(2), February, pp. 201–211.
[89] Balachandran, B., Kalmár-Nagy, T., and Gilsinn, D. E., 2009. Delay Differential
Equations: Recent Advances and New Directions. Spinger.
[90] Hartung, F., Insperger, T., Stépán, G., and Turi, J., 2006. “Approximate stability charts for milling processes using semi-discretization”. Applied Mathematics
and Computations, 174(1), pp. 51–73.
[91] Insperger, T., and Stépán, G., 2004. “Stability analysis of turning with periodic
spindle speed modulation via semi-discretization”. Journal of Vibration and
Control, 10, pp. 1835–1855.
[92] Dombovari, Z., Barton, D. A., Wilson, R. E., and Stepan, G., 2010. “On the
global dynamics of chatter in the orthogonal cutting model”. International
Journal of Non-Linear Mechanics, 46, pp. 330–338.
[93] Balachandran, B., 2001. “Nonlinear dynamics of milling processes”. Philosophical Transactions of the Royal Society A, 359, pp. 793–819.
[94] Bachrathy, D., Stepan, G., and Turi, J., 2011. “State dependent regenerative
effect in milling processes”. Journal of Computational and Nonlinear Dynamics,
6(4), p. 9.
[95] Bachrathy, D., and Stepan, G., 2009. “Bistable parameter region caused by velocity dependent chip thickness in milling process”. In 12th CIRP Conference on
Modelling of Machining Operations. San Sebastian, Spain, 05.07-05.08, no. 14,
ISBN: 978-84-608-0866-4, pp. 867–871.
104
[96] Bachrathy, D., Turi, J., and Stepan, G., 2009. “Analysis of the state dependent
regenerative delay model of the milling process”. In SICON CF: Nonlinear
dynamics, stability, identification and control of systems and structures. Roma,
Italy, 09.21-09.25, pp. 1–3.
[97] Bachrathy, D., and Stepan, G., 2011. “State dependent regenerative effect in
milling processes”. In Proceedings of the 7th European Nonlinear Dynamics
Conference (ENOC 2011): Systems with Time Delay (MS-11). Rome, Italy,
07.24-07.29, no. MS11-21, ISBN: 978-88-906234-2-4, pp. 1–2.
[98] Bachrathy, D., and Stepan, G., 2011. “Fold bifurcation in the state-dependent
delay model of milling”. In ASME 2011 International Design Engineering Technical Conferences (IDETC) and Computers and Information in Engineering
Conference (CIE): 8th International Conference on Multibody Systems, Nonlinear Dynamics, and Control (MSNDC). Washington DC, USA, 08.28-08.31,
DETC2011-48300, pp. 1–7.
[99] Mann, B., Inspergeer, T., Bayly, P., and Stepan, G., 2003. “Stability of upmilling and down-milling, part 2: Experimental verification”. International
Journal of Machine Tools and Manufacture, 43, pp. 35–40.
[100] Stépán, G., 2001. “Modelling nonlinear regenerative effects in metal cutting”.
Phil. Trans. R. Soc. Lond. A, 359, pp. 739–757.
[101] Khasawneh, F., Mann, B., Insperger, T., and Stépán, G., 2009. “Increased
stability of low-speed turning through a distributed force and continuous delay model”. Journal of Comupational and Nonlinear Mechanics, 4(4), p. 12.
No.041003.
[102] Dombovari, Z., Altintas, Y., and Stepan, G., 2010. “The effect of serration on
mechanics and stability of milling cutters”. International Journal of Machine
Tools and Manufacture, 50(6), pp. 511 – 520.
[103] Bachrathy, D., and Stepan, G., 2010. “Time-periodic velocity-dependent process
damping in milling processes”. In 2nd CIRP International Conference on Process
Machine Interactions. Vancouver, Canada, 06.10-06.11, no. M09, ISBN: 978-09866331-0-2, pp. 1–12.
[104] Bachrathy, D., and Stepan, G., 2013. “Improved prediction of stability lobes
with extended multi frequency solution”. CIRP Annals - Manufacturing Technology, 62(1), pp. 411 – 414.
[105] Mandelbrot, B., 1967. “How long is the coast of britain? statistical self-similarity
and fractional dimension”. Science, 156, pp. 636–638.
[106] Vicsek, T., 1992. Fractal growth phenomena. No. ISBN 978-981-02-0668-0. New
Jersey: World Scientific, Singapore.
105
[107] Saupe, D., 1987. “Efficient computation of julia sets and their fractal dimension”.
Physica D: Nonlinear Phenomena, 28(3), pp. 358 – 370.
[108] Wolfram Research, I., 2007. Mathematica, version 6.0 ed. Wolfram Research,
Inc., Champaign, Illinois.
[109] MathWorks, I., 2011. Matlab, r2011b ed. MathWorks, Inc., Natick, Massachusetts, United States.
[110] Doedel, E., Champneys, A., Fairgrieve, T., Kuznetsov, Y., Sandstede, B., and
Wang, X., 1997. Auto97: Continuation and bifurcation software for ordinary
differential equations (with homcont). Tech. rep., Concordia University.
[111] Engelborghs, K., Luzyanina, T., and Roose, D., 2002. “Numerical bifurcation
analysis of delay differential equations using dde-biftool”. ACM Transactions
on Mathematical Software, 28(1), pp. 1–21.
[112] Szalai, R., 2009. “Knut: A continuation and bifurcation software for delaydifferential equations”.
[113] Henderson, M. E., 2002. “Multiple parameter continuation: Computing implicitly defined k-manifolds”. International Journal of Bifurcation and Chaos,
12(3), pp. 451–476.
[114] Dombovari, Z., Iglesias, A., Zatarain, M., and Insperger, T., 2011. “Prediction
of multiple dominant chatter frequencies in milling processes”. International
Journal of Machine Tools and Manufacture, 51, pp. 457–464.
[115] Takács, D., 2011. “Dynamics of towed wheels - nonlinear theory and experiments”. PhD thesis, Budapest University of Technology and Economics, Department of Applied Mechanics.
[116] Mandelbrot, B. B., 1982. The fractal geometry of nature. No. ISBN 0-71-6711869. W. H. Freeman Press.
[117] Tél, T., and Gruiz, M., 2006. Chaotic dynamics: an introduction based on
classical mechanics. Cambridge University Press.
[118] Bachrathy, D., and Stepan, G., 2012. “Bisection method in higher dimensions
and the efficiency number”. Periodica polytechnica. Mechanical engineering,
56(2), pp. 81–86.
[119] Bachrathy, D., and Stepan, G., 2013. “Efficient stability chart computation
for general delayed linear time periodic systems”. In ASME 2013 International
Design Engineering Technical Conferences (IDETC) and Computers and Information in Engineering Conference (CIE): 9th International Conference on
Multibody Systems, Nonlinear Dynamics, and Control (MSNDC). Portland,
Oregon, USA, 08.04-08.07, DETC2013-13660, pp. 1–9.
106