Machining Sounds - Part C: Machining¶
Part C: Machining¶
A FAB26 workshop lesson, companion to Machining Dynamics: Jupyter Notebook Edition (DOI: forthcoming)
Authors: Tony L. Schmitz and Michael F. Gomez
"Chatter is the most obscure and delicate of all problems facing the machinist." F. W. Taylor, On the Art of Cutting Metals (1907)
Part A provided the vocabulary of vibrating systems and Part B described the algorithms for sound analysis. Part C applies both to a milling machine. We describe the milling process following the treatment of Chap. 4 of Schmitz and Smith [1], predict the sound for a stable cut using Part B, and investigate regenerative chatter, a feedback loop that occurs in milling and adds new frequency content to the spectrum for a stable cut. The stability map that separates milling parameters that are stable from those that produce chatter is then developed twice, using two analytical solutions presented side by side, and is tested against the Sect. 4.4 time domain simulation from [1], which we provide with an explanation of the input parameters. Part C closes with a practical understanding of milling stability that combines Parts A, B, and C: listen to a cut, identify the chatter frequency if it exists, and choose a better spindle speed when chatter occurs.
The purposes of Part C are to:
- Describe the milling operation: kinematics, chip thickness, cutting force, and the tooth passing frequency.
- Predict and synthesize the sound of a stable cut based on the operating parameters.
- Explain regeneration, the mechanism of self-excited vibration in milling.
- Compute stability maps by two analytical solutions, the average tooth angle approach and the Fourier series (zero-order) approach, which are described together with their assumptions.
- Apply the
p_4_10_1.mmilling time domain simulation from [1] using all input parameters. - Build a stability map directly from time domain simulations and compare it with the analytical solutions.
- Detect chatter in the frequency content of a milling signal and select a better spindle speed based on the chatter frequency.
Knowledge from Parts A and B is assumed. Equations taken from [1] are provided with their original numbers.
Nomenclature¶
The following symbols are used consistently throughout the lesson, matching the notation of [1]:
| Symbol | Meaning | Units |
|---|---|---|
| $\Omega$ | spindle speed | rpm |
| $N_t$ | number of teeth | dimensionless |
| $f_{tooth}$ | tooth passing frequency, $\Omega N_t / 60$ | Hz |
| $\tau$ | tooth period, $60 / (\Omega N_t)$ | s |
| $\phi$ | cutter rotation angle | deg |
| $\phi_s$, $\phi_e$ | cut start and exit angles | deg |
| $f_t$ | feed per tooth | mm/tooth |
| $h$ | instantaneous chip thickness | mm |
| $b$ | axial depth of cut (chip width) | mm |
| $a$ | radial depth of cut | mm |
| $r$, $d$ | tool radius and diameter | mm |
| $k_t$, $k_n$ | tangential and normal cutting force coefficients | N/mm² |
| $K_s$, $\beta$ | specific force, $\sqrt{k_t^2 + k_n^2}$, and force angle | N/mm², deg |
| $F_t$, $F_n$ | tangential and normal force components | N |
| $F_x$, $F_y$ | fixed-frame force components | N |
| $n$ | vibration along the surface normal | m |
| $f_n$, $k$, $\zeta$ | modal parameters per direction (Part A) | Hz, N/m, dimensionless |
| $\mu_x$, $\mu_y$ | directional orientation factors | dimensionless |
| $b_{lim}$ | limiting axial depth of cut | mm |
| $N$ | stability lobe number | dimensionless |
| $\varepsilon$ | tooth-to-tooth undulation phase | rad |
| $f_c$ | chatter frequency | Hz |
In code, $\Omega$ appears as omega_rpm, $\phi$ as phi, and the modal symbols keep their Part A names.
Setup: Install and Import Libraries¶
Run the install cell once. If this is the first install on your machine, restart the kernel afterward (Kernel menu, Restart), then continue from the import cell.
%pip install -q numpy matplotlib plotly
Note: you may need to restart the kernel to use updated packages.
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from IPython.display import Audio, display
AUDIO_RATE = 48000 # Hz, sample rate for all synthesized sound
print('Libraries loaded.')
Libraries loaded.
The utility functions carried over from Parts A and B are defined in the collapsed cell below and documented in the Appendix. The analytical stability solutions and the time domain simulation are not utilities; they are the content of Sects. C.4 and C.5 and appear there in full view.
# --- Utility functions used throughout the lesson (documented in the Appendix) ---
def fade(x, ms=15, fs=AUDIO_RATE):
"""Apply a short raised-cosine fade-in and fade-out so audio clips start
and stop without clicks."""
n = int(ms * 1e-3 * fs)
ramp = 0.5 * (1 - np.cos(np.pi * np.arange(n) / n))
y = x.copy()
y[:n] *= ramp
y[-n:] *= ramp[::-1]
return y
def compute_spectrum(x, fs, fmax=4000):
"""Amplitude spectrum of signal x on a linear scale, normalized to a peak
value of one. Returns the frequency vector (Hz, limited to fmax) and the
normalized spectrum."""
x = np.asarray(x, dtype=float)
w = np.hanning(len(x))
X = np.fft.rfft((x - x.mean()) * w)
f = np.fft.rfftfreq(len(x), 1 / fs)
amp = 2 * np.abs(X) / np.sum(w)
amp = amp / amp.max()
sel = f <= fmax
return f[sel], amp[sel]
def compute_spectrogram(x, fs, nfft=2048, hop=512, fmax=4000):
"""Short-time amplitude spectrogram on a linear scale, normalized to a
peak value of one. Returns (t_frames, f, S) with one column per frame."""
w = np.hanning(nfft)
n_frames = 1 + max(0, (len(x) - nfft) // hop)
f = np.fft.rfftfreq(nfft, 1 / fs)
sel = f <= fmax
S = np.empty((int(sel.sum()), n_frames))
for j in range(n_frames):
S[:, j] = np.abs(np.fft.rfft(x[j * hop:j * hop + nfft] * w))[sel]
S /= S.max()
t_frames = (np.arange(n_frames) * hop + nfft / 2) / fs
return t_frames, f[sel], S
def decimate_for_plot(t, x, max_points=None):
"""Pass plotting data through unchanged by default (max_points=None).
To thin a long record for display, set max_points to an integer: each
output bin keeps its local minimum and maximum, preserving the visual
envelope without aliasing. Computations always use the full-rate data."""
n = len(x)
if max_points is None or n <= max_points:
return np.asarray(t), np.asarray(x)
bins = max_points // 2
edge = (n // bins) * bins
xb = np.asarray(x)[:edge].reshape(bins, -1)
tb = np.asarray(t)[:edge].reshape(bins, -1)
tt = np.repeat(tb.mean(axis=1), 2)
xx = np.empty(2 * bins)
xx[0::2] = xb.min(axis=1)
xx[1::2] = xb.max(axis=1)
tt = np.append(tt, np.asarray(t)[-1])
xx = np.append(xx, np.asarray(x)[-1])
return tt, xx
def resample_to_audio(x, fs_in, fs_out=AUDIO_RATE):
"""Linearly interpolate a simulated record onto the audio sample grid so
it can be played."""
t_in = np.arange(len(x)) / fs_in
t_out = np.arange(int(t_in[-1] * fs_out)) / fs_out
return np.interp(t_out, t_in, x)
C.1 Milling Description¶
In milling, a rotating tool with defined cutting edges is moved relative to a workpiece to remove material and obtain the desired geometry and dimensions [1]. The tool is mounted in a holder attached to the spindle, and the spindle provides the tool's rotational speed, torque, and power. At minimum, three mutually perpendicular linear axes manipulate the tool-holder-spindle relative to the workpiece, traditionally labeled x, y, and z, with z indicating the tool axis. Machines with additional rotational axes provide contouring capability for non-prismatic parts. Milling machines may be manual or computer numerically controlled (CNC) and both vertical and horizontal spindle configurations are available.
Cutting tools come in many varieties tailored to peripheral, end, contour, and face milling. For our analysis here, we focus on peripheral and end milling, although the concepts extend to other operations. Endmills are loosely categorized by their free-end geometry as square, ball nose, or bull nose, and the teeth are either ground into the body (typically high-speed steel or sintered carbide) or provided as replaceable inserts clamped to a steel body. The cutting edge is not usually parallel to the tool axis; it is inclined with a helix angle $\gamma$ so the chip is spread over an increased edge length and the cutting pressure is reduced. The chip width $\tilde{b}$ is then related to the axial depth of cut $b$ by
$$\tilde{b} = \frac{b}{\cos\gamma} \tag{4.1}$$
For the straight-teeth analysis of this lesson we take $\gamma = 0$, so the chip width and axial depth are equal.
Under the circular tool path approximation, the instantaneous chip thickness varies with the cutter's rotation angle $\phi$:
$$h(\phi) = f_t \sin(\phi) \tag{4.2}$$
where $f_t$ is the feed per tooth, related to the linear feed $f$ (mm/min), the spindle speed $\Omega$ (rpm), and the tooth count $N_t$ by
$$f_t = \frac{f}{N_t \, \Omega} \tag{4.3}$$
The chip thickness is zero at $\phi = 0$ and 180 deg and reaches its maximum of $f_t$ at 90 deg. In conventional, or up, milling the chip thickness increases as the tooth rotates through the cut, while in climb, or down, milling it decreases. The decrease in force as it rotates toward the finished surface is one reason down milling is often selected for finishing passes. The range of engagement angles depends on the radial depth of cut $a$ and tool radius $r$. For up milling, the start angle is $\phi_s = 0$ and the exit angle is
$$\phi_e = \cos^{-1}\!\left(\frac{r - a}{r}\right) \tag{4.4}$$
while for down milling the exit angle is $\phi_e = 180$ deg and the start angle follows from the same geometry. A cut with $a$ equal to half the diameter is called 50 percent radial immersion, and a slot is 100 percent.
Neither the chip thickness nor the chip width is fixed in milling. In a straight slotting cut the chip thickness encountered by each tooth varies continuously as each tooth enters and exits the cut. In partial immersion pocket milling the radial depth may change as well. In sculptured surface milling, even the axial depth varies. The motion of the rotating tool is simple to visualize, but the exact engagement conditions of the cutting edge can be surprisingly complicated.
The cutting force on a tooth is proportional to the chip area through the specific force. The tangential and normal components are
$$F_t = k_t \, b \, h(\phi) \tag{4.7}$$
$$F_n = k_n \, b \, h(\phi) \tag{4.8}$$
where $k_t$ and $k_n$ are the tangential and normal cutting force coefficients for the tool-material pair, equivalently expressed through the specific force $K_s = \sqrt{k_t^2 + k_n^2}$ and force angle $\beta = \tan^{-1}(k_t / k_n)$. These components rotate with the tool. However, measurement and analysis are more convenient in a coordinate system fixed to the machine tool's axes, so they are projected using the cutter angle:
$$F_x = F_t \cos(\phi) + F_n \sin(\phi) \tag{4.9}$$
$$F_y = F_t \sin(\phi) - F_n \cos(\phi) \tag{4.10}$$
When more than one tooth can be engaged, the fixed-frame forces sum over the teeth, with each tooth's angle $\phi_j$ checked against the range of engagement angles $[\phi_s, \phi_e]$ before its contribution is added:
$$F_x = \sum_{j=1}^{N_t} \left[F_t(\phi_j)\cos(\phi_j) + F_n(\phi_j)\sin(\phi_j)\right] \tag{4.15}$$
$$F_y = \sum_{j=1}^{N_t} \left[F_t(\phi_j)\sin(\phi_j) - F_n(\phi_j)\cos(\phi_j)\right] \tag{4.16}$$
The following explorer shows the chip thickness for each of the four teeth and the resulting single-revolution force profiles for a set of cut geometries, computed with the Example 4.2 parameters ($k_t = 750$ and $k_n = 250$ N/mm², $b = 5$ mm, $f_t = 0.1$ mm/tooth, four teeth).
# --- Fig. C.1: chip thickness and forces across cut geometries ---
kt_c, kn_c = 750e6, 250e6 # N/m^2
b_c, ft_c = 5e-3, 0.1e-3 # m, m/tooth
Nt_c = 4
GEOMS = [('25% radial immersion, up', 0.0, 60.0),
('50% radial immersion, up', 0.0, 90.0),
('100% radial immersion (slot)', 0.0, 180.0),
('50% radial immersion, down', 90.0, 180.0),
('25% radial immersion, down', 120.0, 180.0)]
START_G = 2 # slider starts at the slot
phi_deg = np.linspace(0, 360, 1441)
phi_rad = np.deg2rad(phi_deg)
fig = make_subplots(rows=1, cols=2, column_widths=[0.42, 0.58], horizontal_spacing=0.10,
subplot_titles=('Chip thickness, all four teeth (Eq. 4.2)',
'One revolution of force (Eqs. 4.15 and 4.16)'))
for gi, (name, ps, pe) in enumerate(GEOMS):
Fx = np.zeros_like(phi_deg)
Fy = np.zeros_like(phi_deg)
vis = (gi == START_G)
for j in range(Nt_c):
pj_deg = (phi_deg + j * 360 / Nt_c) % 360
pj = np.deg2rad(pj_deg)
eng = (pj_deg >= ps) & (pj_deg <= pe)
hj = np.where(eng, ft_c * np.sin(pj), 0.0)
# Tooth 1 is drawn in full colour and is Eq. 4.2 exactly; teeth 2 to 4 are
# drawn faintly so the sum of Eqs. 4.15 and 4.16 can be read across panels.
style = (dict(color='#1f77b4', width=2) if j == 0
else dict(color='#aec7e8', width=1))
fig.add_trace(go.Scatter(x=phi_deg, y=hj * 1e3, mode='lines',
name=f'tooth {j + 1}', line=style,
showlegend=False, visible=vis), row=1, col=1)
Ft = kt_c * b_c * hj
Fn = kn_c * b_c * hj
Fx += Ft * np.cos(pj) + Fn * np.sin(pj)
Fy += Ft * np.sin(pj) - Fn * np.cos(pj)
fig.add_trace(go.Scatter(x=phi_deg, y=Fx, mode='lines', name='F<sub>x</sub>',
line=dict(color='#1f77b4'), visible=vis), row=1, col=2)
fig.add_trace(go.Scatter(x=phi_deg, y=Fy, mode='lines', name='F<sub>y</sub>',
line=dict(color='#ff7f0e', dash='dash'),
visible=vis), row=1, col=2)
fig.add_trace(go.Scatter(x=phi_deg, y=np.sqrt(Fx**2 + Fy**2), mode='lines',
name='resultant',
line=dict(color='#2ca02c', dash='dot'),
visible=vis), row=1, col=2)
n_traces_g = Nt_c + 3 # four chip thickness curves, then Fx, Fy, resultant
steps = []
for gi, (name, ps, pe) in enumerate(GEOMS):
visv = [False] * (n_traces_g * len(GEOMS))
visv[n_traces_g * gi:n_traces_g * (gi + 1)] = [True] * n_traces_g
steps.append(dict(method='update', label=name, args=[{'visible': visv}]))
fig.update_layout(
sliders=[dict(active=START_G, currentvalue=dict(prefix='cut geometry: '), steps=steps)],
title='Fig. C.1 — Chip thickness and cutting force across cut geometries',
template='simple_white', height=430,
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=70, b=95))
fig.update_xaxes(title_text='Cutter Angle, φ (deg)', row=1, col=1)
fig.update_xaxes(title_text='Cutter Angle, φ (deg)', row=1, col=2)
fig.update_yaxes(title_text='Chip Thickness, h (mm)', row=1, col=1)
fig.update_yaxes(title_text='Force (N)', row=1, col=2)
fig.show()
Two observations from the explorer assist with our future analyses. Down milling mirrors up milling in time, with the force reducing as the cutter rotates toward the finished surface. Also, the x direction force changes sign between the two. For the slot, once the teeth are engaged continuously the summed force approaches a steadier profile; slotting with an even tooth count above two produces constant total force.
C.1.1 Tooth Passing Frequency¶
The previous forces are plotted against angle. We can use the spindle speed to convert angle to time and the force record becomes a periodic train of pulses, just like we discussed in Sect. B.7. Its fundamental frequency is the tooth passing frequency:
$$f_{tooth} = \frac{\Omega \, N_t}{60} \tag{4.14}$$
with $\Omega$ in rpm. The spectrum of the cutting force therefore contains $f_{tooth}$ and its integer multiples (harmonics), plus DC content from the nonzero mean value. The magnitudes of the individual harmonics depend on the sharpness of the force pulses and therefore the cutting conditions. A lower radial immersion results in more impulse-like forces and larger magnitudes for higher frequency harmonics. Because the forces excite the structure, a stable cut produces forced vibration at these frequencies. Changing the spindle speed moves the entire excitation family composed of the tooth passing frequency and its harmonics. The following cell lets you hear it.
# --- The tooth passing rhythm. Adjust and rerun. ---
omega_rpm = 7500 # spindle speed (rpm). Try changing these!
Nt = 4 # number of teeth
f_tooth = omega_rpm * Nt / 60
print(f'Eq. 4.14: f_tooth = {omega_rpm} x {Nt} / 60 = {f_tooth:.0f} Hz '
f'(spindle rotation at {omega_rpm / 60:.1f} Hz)')
dur = 1.5
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
phase = (t * f_tooth) % 1.0
pulses = np.where(phase < 0.12, np.sin(np.pi * phase / 0.12), 0.0)
display(Audio(fade(0.6 * (pulses - pulses.mean())), rate=AUDIO_RATE))
Eq. 4.14: f_tooth = 7500 x 4 / 60 = 500 Hz (spindle rotation at 125.0 Hz)
If the radial depth of cut is very small, there are long periods when no tooth is cutting at all and the force record is a series of very short impacts. In the frequency domain, this results in higher frequency harmonics with larger magnitudes. At larger radial depths the number of simultaneously engaged teeth may be constant or may alternate, depending on the geometry and tooth spacing. Interpreting the spectrum of a machining signal begins with Eq. 4.14: compute $f_{tooth}$ from the programmed spindle speed and tooth count and then every peak (or line) of the stable cut family should appear at an integer multiple of it.
▸ Key terms: milling, up milling, down milling, feed per tooth, radial immersion, chip thickness, tooth passing frequency
C.2 The Sound of a Stable Cut¶
Part B ended with a prediction: a repeating cutting force must produce harmonics at multiples of its fundamental frequency. A real machine adds two more ingredients. First, the idle machine has a baseline of broadband noise and amplifier and motor-related tones that exists without any cutting. Second, a real cutter never runs perfectly true: runout, the inequality in the radii of the individual teeth measured from the tool's axis of rotation, makes the teeth cut slightly unequal chips, which adds content at the once-per-revolution frequency $\Omega/60$ and its harmonics, spaced $N_t$ times more finely than the tooth passing family [1]. The following explorer describes a stable cut using one ingredient at a time, for the Example 4.11 conditions (7500 rpm, four teeth, so $f_{tooth} = 500$ Hz and the runout spacing is 125 Hz).
# --- Fig. C.2: a stable cut assembled ingredient by ingredient ---
dur = 2.0
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
rng = np.random.default_rng(26)
f_rev = 7500 / 60 # 125 Hz once-per-revolution
f_tp = 500.0 # tooth passing fundamental
baseline = 0.06 * rng.standard_normal(len(t))
for f_m in (60, 120, 240): # motor-related tones
baseline += 0.03 * np.sin(2 * np.pi * f_m * t)
runout = np.zeros_like(t)
for m in range(1, 8):
if m % 4: # skip multiples that overlap f_tooth
runout += (0.12 / m) * np.sin(2 * np.pi * m * f_rev * t)
tooth = np.zeros_like(t)
for m in range(1, 9):
tooth += (0.8 / m) * np.sin(2 * np.pi * m * f_tp * t)
STAGES = [('idle machine', baseline),
('idle + spindle runout', baseline + runout),
('idle + runout + tooth passing', baseline + runout + tooth)]
fig = go.Figure()
for si, (name, sig) in enumerate(STAGES):
f_sp, a_sp = compute_spectrum(sig, AUDIO_RATE, fmax=2600)
fig.add_trace(go.Scatter(x=f_sp, y=a_sp, mode='lines',
line=dict(color='#1f77b4'),
visible=(si == 2), showlegend=False))
steps = []
for si, (name, sig) in enumerate(STAGES):
vis = [False] * len(STAGES)
vis[si] = True
steps.append(dict(method='update', label=name, args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=2, currentvalue=dict(prefix='ingredients: '), steps=steps)],
title='Fig. C.2 — The spectrum of a stable cut, ingredient by ingredient',
template='simple_white', height=400,
xaxis_title='Frequency (Hz)', yaxis_title='Amplitude (normalized)',
yaxis_range=[0, 1.05],
margin=dict(l=60, r=20, t=60, b=95))
fig.show()
for name, sig in STAGES:
print(name + ':')
display(Audio(fade(0.6 * sig / np.max(np.abs(sig))), rate=AUDIO_RATE))
idle machine:
idle + spindle runout:
idle + runout + tooth passing:
Every line in the final spectrum has a name. The dense low family starting at 125 Hz is the spindle's once-per-revolution content from runout with its multiples, the tall family starting at 500 Hz is the tooth passing frequency with its multiples, and everything else is the idle baseline. This is the complete sound of a stable cut and the practical skill is interpreting its content. We also want to identify anything that does not match this picture. Any large magnitude lines that do not belong to the stable cutting family indicate something else is happening. Sect. C.7 provides the explanation.
▸ Key terms: runout, once-per-revolution frequency
C.3 Regenerative Chatter in Milling¶
Because endmills, holders, and spindles are not rigid, the cutting force deflects the tool, and, if the tool vibrates while removing material, the vibration is imprinted on the workpiece surface as a wavy profile. In milling, the time-delayed surface regeneration step occurs from tooth to tooth: the wavy surface left behind by tooth 1 is removed by tooth 2, one tooth period $\tau = 60/(\Omega N_t)$ later [1]. The instantaneous chip thickness therefore depends on both the current vibration and the surface left by the previous tooth, so vibration feeds back into force, which feeds back into vibration. This closed loop feedback behavior is called regeneration, and it makes chatter possible in milling.
The chip thickness is measured along the surface normal, whose direction rotates with the cutter angle. Projecting the x and y vibrations onto the normal,
$$n = x \sin(\phi) - y \cos(\phi) \tag{4.21}$$
with positive $n$ out of the cut, the instantaneous chip thickness becomes
$$h(t) = f_t \sin(\phi) + n(t - \tau) - n(t) \tag{4.22}$$
where the $n(t-\tau)$ term is the previous tooth's contribution along the current normal. Whether the feedback loop results in growing or shrinking vibrations depends on the phase between the wave left by one tooth and the vibration of the next tooth. When the two are in phase, the chip thickness variation is no worse than the circular tool path alone produces, and the cut produces forced vibration. In fact, matching the tooth passing frequency to the natural frequency encourages exactly this favorable phasing. When the passes are out of phase, however, the chip thickness varies with time, the force varies with it, and the vibration can grow tooth by tooth. This represents self-excited vibration, the third kind from Sect. A.2. The following explorer shows the two surfaces and the chip thickness produced between them as the phase $\varepsilon$ varies.
# --- Fig. C.3: two successive tooth passes and the chip between them ---
u = np.linspace(0, 4 * np.pi, 500)
amp = 1.0
EPS_STEPS = list(range(0, 361, 30))
START_E = 6 # slider starts at 180 deg
fig = go.Figure()
prev = amp * np.sin(u)
for ei, eps in enumerate(EPS_STEPS):
curr = amp * np.sin(u + np.deg2rad(eps)) - 3.0 # offset: one feed below
vis = (ei == START_E)
fig.add_trace(go.Scatter(x=u, y=prev, mode='lines',
name='surface left by previous tooth',
line=dict(color='#1f77b4'), visible=vis))
fig.add_trace(go.Scatter(x=u, y=curr, mode='lines',
name='path of current tooth',
line=dict(color='#ff7f0e', dash='dash'),
visible=vis))
fig.add_trace(go.Scatter(x=np.concatenate([u, u[::-1]]),
y=np.concatenate([prev, curr[::-1]]),
fill='toself', fillcolor='rgba(44,160,44,0.18)',
line=dict(width=0), name='chip', visible=vis,
hoverinfo='skip'))
steps = []
for ei, eps in enumerate(EPS_STEPS):
vis = [False] * (3 * len(EPS_STEPS))
vis[3 * ei:3 * ei + 3] = [True] * 3
steps.append(dict(method='update', label=f'{eps}', args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_E,
currentvalue=dict(prefix='phase between passes (deg) = '),
steps=steps)],
title='Fig. C.3 — The chip between successive tooth passes',
template='simple_white', height=400,
xaxis=dict(visible=False), yaxis=dict(visible=False),
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=20, r=20, t=60, b=95))
fig.show()
At 0 or 360 deg the two waves are aligned and the chip has nearly constant thickness. Even though the vibration exists, it does no regenerative harm. Near 180 deg the waves oppose and the chip thickness varies between nearly zero and nearly double, so the force varies with it and adds energy to the vibration. Whether that energy overcomes the system's damping depends on the axial depth $b$, which scales the force together with the chip thickness. The next section computes the boundary between stable conditions and regenerative chatter.
▸ Key terms: regeneration, tooth period, tooth-to-tooth phase, chatter
C.4 Stability Maps¶
The boundary between stable machining and chatter, drawn as the limiting axial depth $b_{lim}$ against spindle speed $\Omega$, is the stability map (also called a stability lobe diagram due to the shape of the individual curves, or lobes, that are combined to produce the final stability boundary). This section computes it twice, by two analytical solutions, each stated with its central assumption. Both linearize the same physics; they differ in how the rotating, time-varying cut geometry is reduced to a constant one. Each is summarized by its governing equations. The implementations appear in collapsed cells, one expand-click away, and both boundaries are drawn together on one figure.
C.4.1 Analytical Solution 1: The Average Tooth Angle Approach¶
The first approach evaluates the cut geometry at a single representative angle, the average angle of a tooth in the cut, $\phi_{ave} = (\phi_s + \phi_e)/2$. The x and y frequency response functions (FRFs) are combined into one oriented FRF through directional orientation factors obtained by projecting the cutting force onto each mode direction and that result onto the average surface normal [1]:
$$FRF_{orient} = \mu_x \, FRF_x + \mu_y \, FRF_y \tag{4.27}$$
For a slot, $\phi_{ave} = 90$ deg and the projections give $\mu_x = \cos\beta$ and $\mu_y = 0$; the general projections are computed in the code below. The stability boundary is then a mapping of the real part of the oriented FRF, valid where it is negative:
$$b_{lim} = \frac{-1}{2 \, K_s \, \mathrm{Re}\!\left[FRF_{orient}\right] N_t^*} \tag{4.23}$$
where $N_t^* = N_t(\phi_e - \phi_s)/360$ is the average number of teeth in the cut. Each valid chatter frequency $f_c$ maps to one spindle speed per lobe number $N$ through the tooth-to-tooth phase $\varepsilon$, with the best speeds at $\varepsilon = 2\pi$:
$$\Omega = \frac{60 \, f_c}{N_t \left(N + \varepsilon / 2\pi\right)} \tag{4.24}$$
# --- Analytical solution 1: average tooth angle (Tlusty) stability lobes ---
def frf_sdof(f, fn, k, zeta):
r = f / fn
return (1 / k) / ((1 - r**2) + 1j * 2 * zeta * r)
def tlusty_lobes(fn, k, zeta, Ks, beta_deg, Nt, phis_deg, phie_deg,
lobes=range(0, 6), n_f=4000):
'''Average tooth angle stability lobes per Sect. 4.3.1 and 4.3.2 of [1].
Assumption: the rotating cut geometry is evaluated at the single average
tooth angle. Returns a list of (omega_rpm, blim_mm) curves, one per lobe.'''
beta = np.deg2rad(beta_deg)
phi_ave = np.deg2rad(0.5 * (phis_deg + phie_deg))
theta_n = phi_ave - np.pi / 2 # surface normal from +x
mu_x = np.cos(beta + theta_n) * np.cos(theta_n)
mu_y = np.sin(beta + theta_n) * np.sin(theta_n)
Nt_star = Nt * (phie_deg - phis_deg) / 360
fc = np.linspace(fn * 1.0001, fn * 4, n_f)
FRFor = mu_x * frf_sdof(fc, fn, k, zeta) + mu_y * frf_sdof(fc, fn, k, zeta)
R, I = FRFor.real, FRFor.imag
valid = R < 0
fc, R, I = fc[valid], R[valid], I[valid]
blim = -1.0 / (2 * Ks * R * Nt_star) * 1e3 # m -> mm
eps = 2 * np.arctan2(-R, I) # tooth-to-tooth phase (rad)
curves = []
for N in lobes:
omega = 60 * fc / (Nt * (N + eps / (2 * np.pi))) # Eq. 4.24
curves.append((omega, blim, fc))
return curves
print('tlusty_lobes defined: the average tooth angle analytical solution.')
tlusty_lobes defined: the average tooth angle analytical solution.
C.4.2 Analytical Solution 2: The Fourier Series (Zero-Order) Approach¶
The second approach, presented by Altintas and Budak [2] and developed in Sect. 4.3.3 of [1], keeps the full time-varying directional coefficients of the milling equations and expands them in a Fourier series over the tooth period, but truncates the series at the average term, the zero-order approximation, to remove the time dependence. The averaged directional matrix is
$$[A_0] = \frac{N_t}{2\pi}\begin{bmatrix} \alpha_{xx} & \alpha_{xy} \\ \alpha_{yx} & \alpha_{yy} \end{bmatrix} \tag{4.44}$$
with, for example, the xx entry given by
$$\alpha_{xx} = \frac{1}{2}\left[\cos 2\phi - 2 K_r \phi + K_r \sin 2\phi\right]_{\phi_s}^{\phi_e} \tag{4.45}$$
and the companion $\alpha_{xy}$, $\alpha_{yx}$, and $\alpha_{yy}$ integrals written out in the collapsed cell. Stability of the resulting time-invariant system is an eigenvalue problem: the characteristic equation of the closed loop (Eq. 4.50 of [1]) is quadratic in the eigenvalue $\Lambda$,
$$a_0 \Lambda^2 + a_1 \Lambda + 1 = 0, \qquad a_0 = FRF_{xx}\,FRF_{yy}\left(\alpha_{xx}\alpha_{yy} - \alpha_{xy}\alpha_{yx}\right), \qquad a_1 = \alpha_{xx}\,FRF_{xx} + \alpha_{yy}\,FRF_{yy}$$
and each root $\Lambda = \Lambda_{Re} + i\Lambda_{Im}$ gives a candidate limiting depth and phase (Sect. 4.3.3 of [1]):
$$b_{lim} = -\frac{2\pi \Lambda_{Re}}{N_t k_t}\left(1 + \kappa^2\right), \qquad \kappa = \frac{\Lambda_{Im}}{\Lambda_{Re}}, \qquad \psi = \tan^{-1}\kappa, \qquad \varepsilon = \pi - 2\psi, \qquad \Omega = \frac{60 f_c}{N_t\left(\varepsilon/2\pi + j\right)}$$
for lobe numbers $j = 0, 1, 2, \ldots$ Only positive depths are physically meaningful and are retained. The composite boundary is the minimum value from the two eigenvalues at each spindle speed, exactly as Fig. 4.34 of [1] plots it. One detail matters here and Sect. 4.3.3 of [1] states it explicitly: the full chatter frequency range is scanned, including frequencies below the natural frequency, because the coupled two direction system can chatter there. Both approaches share the accuracy caveat stated in Chap. 4 of [1]: as radial immersion reduces and the force becomes more impulse-like, the accuracy obtained by averaging (by either approach) degrades. The following collapsed cells implement the algorithms; expand them whenever wanted.
# --- Analytical solution 2: Fourier series (zero-order) stability lobes ---
def zoa_lobes(fn, k, zeta, kt, kn, Nt, phis_deg, phie_deg,
lobes=range(0, 13), n_f=40000):
'''Fourier series (zero-order) stability lobes per Sect. 4.3.3 of [1],
following Altintas and Budak [2]. Only the average term of the Fourier
series of the time varying directional coefficients is retained (Eqs.
4.44 and 4.45). At each candidate chatter frequency the quadratic
characteristic equation a0*Lam^2 + a1*Lam + 1 = 0 gives two eigenvalues
(Eq. 4.53); each eigenvalue gives a candidate limiting depth (Eq. 4.58),
and only positive depths are physical and retained; the regenerative
phase then assigns one spindle speed per lobe number j (Eq. 4.59). Per
[1], the full chatter frequency range is scanned, including frequencies
below fn, where the coupled two direction system can also chatter.
Returns a list of (omega_rpm, blim_mm, fc_hz) branches.'''
Kr = kn / kt
def alpha_anti(p):
# antiderivatives of the Eq. 4.45 integrands, evaluated at angle p (rad)
return np.array([0.5 * ( np.cos(2*p) - 2*Kr*p + Kr*np.sin(2*p)), # xx
0.5 * (-np.sin(2*p) - 2*p + Kr*np.cos(2*p)), # xy
0.5 * (-np.sin(2*p) + 2*p + Kr*np.cos(2*p)), # yx
0.5 * (-np.cos(2*p) - 2*Kr*p - Kr*np.sin(2*p))]) # yy
axx, axy, ayx, ayy = (alpha_anti(np.deg2rad(phie_deg))
- alpha_anti(np.deg2rad(phis_deg)))
fc = np.linspace(30.0, fn * 8, n_f) # full frequency range
Phix = frf_sdof(fc, fn, k, zeta)
Phiy = frf_sdof(fc, fn, k, zeta)
a0 = Phix * Phiy * (axx * ayy - axy * ayx)
a1 = axx * Phix + ayy * Phiy
disc = np.sqrt(a1**2 - 4 * a0 + 0j)
curves = []
for Lam in ((-a1 + disc) / (2 * a0), (-a1 - disc) / (2 * a0)):
LR, LI = Lam.real, Lam.imag
with np.errstate(all='ignore'):
kappa = LI / LR
blim = -(2 * np.pi * LR / (Nt * kt)) * (1 + kappa**2) * 1e3
keep = np.isfinite(blim) & (blim > 0) # positive depths only
edges = np.where(np.diff(keep))[0] + 1 # contiguous valid runs
for run in np.split(np.arange(len(fc)), edges):
if len(run) < 2 or not keep[run[0]]:
continue
psi = np.arctan(kappa[run])
eps = np.pi - 2 * psi
for j in lobes:
om = 60 * fc[run] / (Nt * (eps / (2 * np.pi) + j))
curves.append((om, blim[run], fc[run]))
return curves
def lower_envelope(curves, om_grid):
'''Composite stability boundary: the minimum limiting depth over all lobe
branches at each spindle speed, as [1] plots in its Fig. 4.34. Each
branch is a continuous curve traced by the chatter frequency, so the
minimum is taken over linear interpolations along each branch, split
into monotone pieces wherever a branch folds back in spindle speed.
Interpolating along the curve matters on the branch tails, where a small
chatter frequency step sweeps a large speed range; collecting isolated
points into speed bins under-samples those tails and produces jagged,
grid-dependent artifacts at high spindle speed.'''
env = np.full(len(om_grid), np.inf)
for om, bl, fc in curves:
d = np.sign(np.diff(om))
brk = np.where(d[1:] != d[:-1])[0] + 1
for p in np.split(np.arange(len(om)), brk + 1):
if len(p) < 2:
continue
o, b = om[p], bl[p]
if o[0] > o[-1]:
o, b = o[::-1], b[::-1]
env = np.minimum(env, np.interp(om_grid, o, b,
left=np.inf, right=np.inf))
return env
print('zoa_lobes and lower_envelope defined: the Fourier series solution.')
zoa_lobes and lower_envelope defined: the Fourier series solution.
With both solutions in hand, we draw them together for the Example 4.11 system (fn = 500 Hz, k = 8 MN/m, ζ = 0.02 in both directions, four teeth, slotting, kt = 695 and kn = 281 N/mm², so Ks = 750 N/mm² and β = 68 deg), and mark the four operating points the simulation of Sect. C.5 will test.
# --- Fig. C.4: both analytical boundaries with the Example 4.11 cases ---
PAR = dict(fn=500.0, k=8e6, zeta=0.02, Nt=4, phis_deg=0.0, phie_deg=180.0)
KT, KN = 695e6, 281e6
KS = np.sqrt(KT**2 + KN**2)
BETA = np.rad2deg(np.arctan2(KT, KN))
CASES = [('1', 7500, 3.0, 'stable'), ('2', 7500, 5.0, 'unstable'),
('3', 5000, 0.1, 'stable'), ('4', 5000, 0.5, 'unstable')]
tl = tlusty_lobes(PAR['fn'], PAR['k'], PAR['zeta'], KS, BETA,
PAR['Nt'], PAR['phis_deg'], PAR['phie_deg'])
zo = zoa_lobes(PAR['fn'], PAR['k'], PAR['zeta'], KT, KN,
PAR['Nt'], PAR['phis_deg'], PAR['phie_deg'])
zo_om = np.linspace(2000, 16000, 561)
zo_bl = lower_envelope(zo, zo_om)
fig = go.Figure()
for i, (om, bl, fc) in enumerate(tl):
fig.add_trace(go.Scatter(x=om, y=bl, mode='lines',
line=dict(color='#1f77b4'),
name='average tooth angle', showlegend=(i == 0),
legendgroup='tl'))
fig.add_trace(go.Scatter(x=zo_om, y=zo_bl, mode='lines',
line=dict(color='#ff7f0e', dash='dash'),
name='Fourier series (zero-order)'))
for label, om, b, verdict in CASES:
color = '#2ca02c' if verdict == 'stable' else '#d62728'
fig.add_trace(go.Scatter(x=[om], y=[b], mode='markers+text',
marker=dict(color=color, size=11,
symbol='circle' if verdict == 'stable' else 'x'),
text=[label], textposition='top center',
showlegend=False))
fig.update_layout(title='Fig. C.4 — Stability lobes by two analytical solutions',
template='simple_white', height=470,
xaxis=dict(title='Spindle Speed (rpm)', range=[2000, 16000]),
yaxis=dict(title='Axial Depth of Cut, b (mm)', range=[0, 8]),
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=60, b=50))
fig.show()
def blim_at(curves, omega_rpm):
best = np.inf
for om, bl, fc in curves:
sel = (om > omega_rpm - 60) & (om < omega_rpm + 60)
if sel.any():
best = min(best, np.min(bl[sel]))
return best
print('Boundary near each case, and the verdict each method implies:')
for label, om, b, verdict in CASES:
bt = blim_at(tl, om)
bz = float(np.interp(om, zo_om, zo_bl))
vt = 'stable' if b < bt else 'unstable'
vz = 'stable' if b < bz else 'unstable'
print(f'case {label} ({om} rpm, {b} mm): '
f'average tooth angle {bt:5.2f} mm implies {vt}; '
f'zero-order {bz:5.2f} mm implies {vz}; Example 4.11 of [1] reports {verdict}')
Boundary near each case, and the verdict each method implies: case 1 (7500 rpm, 3.0 mm): average tooth angle 8.99 mm implies stable; zero-order 3.80 mm implies stable; Example 4.11 of [1] reports stable case 2 (7500 rpm, 5.0 mm): average tooth angle 8.99 mm implies stable; zero-order 3.80 mm implies unstable; Example 4.11 of [1] reports unstable case 3 (5000 rpm, 0.1 mm): average tooth angle 0.97 mm implies stable; zero-order 0.23 mm implies stable; Example 4.11 of [1] reports stable case 4 (5000 rpm, 0.5 mm): average tooth angle 0.97 mm implies stable; zero-order 0.23 mm implies unstable; Example 4.11 of [1] reports unstable
The two solutions disagree and the disagreement is the lesson of this section. For slotting, the average tooth angle is $\phi_{ave} = 90$ deg, which places the average surface normal parallel to the x direction. The projection of the y direction onto that normal is then zero, $\mu_y = 0$, and the average tooth angle approach discards the y direction dynamics entirely; its stability boundary reflects the x direction alone and reaches nearly 9 mm at the 7500 rpm pocket. The zero-order approach retains both directions through the off-diagonal terms of Eq. 4.44, and its boundary is lower, about 3.7 mm at 7500 rpm and about 0.23 mm at 5000 rpm, with wide low bands on both sides of the pocket. Much of the reduction comes from chatter frequencies below the natural frequency, which exist only because the two directions are coupled. Past the last pocket the boundary holds near 0.23 mm across the rest of this diagram, so raising the spindle speed alone does not raise the limiting depth for this system; the deeper stable regions are the pockets. The previous evaluation compares the four cases against each stability map: the average tooth angle boundary predicts that cases 2 and 4 will be stable, while the zero-order boundary reproduces all four Example 4.11 outcomes. Sect. 4.4 of [1] states the reason plainly: the disagreement for cases 2 and 4 is due to the orthogonality between the average surface normal and the y direction in slotting, which forces $\mu_y = 0$ and neglects the contribution of the y direction dynamics. Averaging the geometry at a single angle is a modeling assumption, and here it costs roughly a factor of two in the predicted limit for slotting. Sect. C.5 applies a time domain simulation, which makes no averaging assumption, to decide between them.
Both boundaries are useful in practice: near spindle speeds where the tooth passing frequency divides evenly into the natural frequency, the stable pockets reach several times the depth available elsewhere, and case 1 appears inside such a pocket. The minimum axial depth that provides the floor of the whole diagram depends on the damping. The following explorer computes both boundaries at three damping ratios; drag between them and watch the valleys rise and fall together.
# --- Fig. C.5: the stability boundary as damping varies ---
ZETAS = [0.01, 0.02, 0.05]
START_Z = 1
fig = go.Figure()
n_per = []
for zi, zval in enumerate(ZETAS):
tl_z = tlusty_lobes(PAR['fn'], PAR['k'], zval, KS, BETA,
PAR['Nt'], PAR['phis_deg'], PAR['phie_deg'])
zo_bl_z = lower_envelope(zoa_lobes(PAR['fn'], PAR['k'], zval, KT, KN,
PAR['Nt'], PAR['phis_deg'],
PAR['phie_deg']), zo_om)
vis = (zi == START_Z)
for i, (om, bl, fc) in enumerate(tl_z):
fig.add_trace(go.Scatter(x=om, y=bl, mode='lines',
line=dict(color='#1f77b4'),
name='average tooth angle',
showlegend=(i == 0), legendgroup='tl',
visible=vis))
fig.add_trace(go.Scatter(x=zo_om, y=zo_bl_z, mode='lines',
line=dict(color='#ff7f0e', dash='dash'),
name='Fourier series (zero-order)',
legendgroup='zo', visible=vis))
n_per.append(len(tl_z) + 1)
steps = []
total = sum(n_per)
for zi, zval in enumerate(ZETAS):
vis = [False] * total
ofs = sum(n_per[:zi])
vis[ofs:ofs + n_per[zi]] = [True] * n_per[zi]
steps.append(dict(method='update', label=f'{zval:.2f}', args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_Z, currentvalue=dict(prefix='Damping Ratio, ζ = '), steps=steps)],
title='Fig. C.5 — Damping sets the floor of the stability lobe diagram',
template='simple_white', height=440,
xaxis=dict(title='Spindle Speed (rpm)', range=[2000, 16000]),
yaxis=dict(title='Axial Depth of Cut, b (mm)', range=[0, 10]),
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=60, b=95))
fig.show()
▸ Key terms: stability map, stability lobe diagram, limiting axial depth, oriented FRF, directional orientation factors, average tooth angle approach, zero-order approximation
C.5 The Milling Time Domain Simulation¶
The analytical solutions take little time to compute, but are based on averaging assumptions. The time domain simulation of Sect. 4.4 of [1] makes neither approximation: it marches the cut forward in small time steps, keeping the full rotating geometry, the tooth-by-tooth surface regeneration of Eq. 4.22, and the nonlinearity that occurs when a tooth vibrates at a large enough amplitude that it leaves the cut entirely. This section provides a Python replication of p_4_10_1.m from [1], kept in a collapsed cell and expandable whenever the line-by-line detail is wanted.
The strategy divides one revolution into steps_rev equal angle increments $d\phi = 360/\texttt{steps\_rev}$ deg, with time step $dt = 60/(\texttt{steps\_rev} \cdot \Omega)$ s, and requires steps_rev to be a multiple of $N_t$ so the teeth land exactly on prespecified grid angles. A vector indexed by grid angle stores the surface left by the previous tooth at each orientation. Every time step performs three activities [1]: rotate the cutter by one increment; for each tooth inside the angular range $[\phi_s, \phi_e]$, calculate the chip thickness from Eq. 4.22 using the stored surface and the current normal-direction vibration (Eq. 4.21), compute the forces by Eqs. 4.7 to 4.10 with the stored surface updated to the current vibration, or, if the chip thickness is negative, set the forces to zero and add the feed contribution to the stored surface because the tooth has vibrated out of the cut; and numerically integrate the x and y equations of motion for that time step. The following function contains the complete simulation.
# --- The milling time domain simulation (Python equivalent of p_4_10_1.m) ---
def mill_tds(omega_rpm, b_mm, Nt=4, ft_mm=0.15, phis_deg=0.0, phie_deg=180.0,
fn_x=500.0, k_x=8e6, zeta_x=0.02,
fn_y=500.0, k_y=8e6, zeta_y=0.02,
kt=695e6, kn=281e6, steps_rev=652, revolutions=20):
"""Milling time domain simulation with straight teeth, per Sect. 4.4 of
[1]. Returns a dict with time (s), forces Fx, Fy (N), displacements
x, y (m), the number of steps per tooth period, and the sample rate (Hz).
steps_rev is rounded to the nearest multiple of Nt as in p_4_10_1.m."""
steps_rev = int(round(steps_rev / Nt) * Nt)
dphi = 2 * np.pi / steps_rev
dt = 60.0 / (steps_rev * omega_rpm)
steps = steps_rev * revolutions
b = b_mm * 1e-3
ft = ft_mm * 1e-3
ps, pe = np.deg2rad(phis_deg), np.deg2rad(phie_deg)
phi = np.arange(steps_rev) * dphi
sin_p, cos_p = np.sin(phi), np.cos(phi)
in_cut = (phi >= ps) & (phi <= pe)
teeth = (np.arange(Nt) * (steps_rev // Nt)).astype(int)
surf = np.zeros(steps_rev)
wnx, wny = 2 * np.pi * fn_x, 2 * np.pi * fn_y
mx, my = k_x / wnx**2, k_y / wny**2
cx, cy = 2 * zeta_x * np.sqrt(mx * k_x), 2 * zeta_y * np.sqrt(my * k_y)
Fx = np.zeros(steps); Fy = np.zeros(steps)
xh = np.zeros(steps); yh = np.zeros(steps)
x = y = vx = vy = 0.0
for cnt in range(steps):
teeth = (teeth + 1) % steps_rev
fx_s = fy_s = 0.0
for j in range(Nt):
idx = teeth[j]
if not in_cut[idx]:
continue
s, c = sin_p[idx], cos_p[idx]
n = x * s - y * c # Eq. 4.21
h = ft * s + surf[idx] - n # Eq. 4.22
if h < 0: # tooth vibrated out of the cut
surf[idx] += ft * s
continue
Ft = kt * b * h # Eq. 4.7
Fn = kn * b * h # Eq. 4.8
surf[idx] = n
fx_s += Ft * c + Fn * s # Eq. 4.9
fy_s += Ft * s - Fn * c # Eq. 4.10
Fx[cnt], Fy[cnt] = fx_s, fy_s
vx += dt * (fx_s - cx * vx - k_x * x) / mx # Euler integration, x
x += dt * vx
vy += dt * (fy_s - cy * vy - k_y * y) / my # Euler integration, y
y += dt * vy
xh[cnt], yh[cnt] = x, y
return dict(t=np.arange(steps) * dt, Fx=Fx, Fy=Fy, x=xh, y=yh,
steps_tooth=steps_rev // Nt, fs=1.0 / dt, steps_rev=steps_rev)
def chatter_metric(sim, threshold_um=1.0):
"""Periodic sampling stability metric per Sect. 4.4.5 of [1]. The x
displacement is sampled once per tooth passage over the second half of
the record, after the startup transient has decayed, and M is the
mean absolute difference between
successive samples, in micrometers. Stable cuts repeat at the tooth
period, so M sits near zero; chatter does not repeat and M reaches tens
to hundreds of micrometers ([1] demonstrates 0.06 versus 45.56 for its
Example 4.13). Cuts with M below one micrometer are classified stable.
Returns (M_um, verdict)."""
x = sim['x']
start = len(x) // 2 # after the startup transient
xs = x[np.arange(start, len(x), sim['steps_tooth'])] * 1e6
M = float(np.sum(np.abs(np.diff(xs))) / len(xs))
return M, ('stable' if M < threshold_um else 'chatter')
print('mill_tds and chatter_metric defined.')
mill_tds and chatter_metric defined.
The following parameter panel executes the simulation. Each input can be modified; the defaults are from Example 4.11 of [1], and the PRESET line selects one of the four Example 4.11 cases. Change any value and rerun. The output applies a stability metric to whatever parameters are set; the following figure runs all four cases and shows the x direction displacement for each with periodically sampled points added as circles. It applies the stability metric from Sect. 4.4.5 of [1]: the displacement is sampled once per tooth passage over the second half of the record, after the startup transient has decayed, and $M$ is the mean absolute difference between successive samples, in micrometers. For a stable cut, the samples repeat and $M$ is ideally zero; for chatter they do not repeat and $M$ reaches tens to hundreds of micrometers (Example 4.13 of [1] compares 0.06 for stable to 45.56 for unstable). A threshold value of one micrometer classifies every cut in this Part as stable (less than or equal to one) or chatter (more than one). Slotting with an even number of teeth makes the stable side especially quiet: Sect. C.1 showed the total force is constant there, so a stable cut carries almost no vibration and its samples are nearly identical.
# --- Run the simulation. Adjust the parameters and rerun. ---
PRESETS = {'1': dict(omega_rpm=7500, b_mm=3.0), # stable, per [1]
'2': dict(omega_rpm=7500, b_mm=5.0), # unstable
'3': dict(omega_rpm=5000, b_mm=0.1), # stable
'4': dict(omega_rpm=5000, b_mm=0.5)} # unstable
PRESET = '2' # pick a case, or set the parameters directly below
params = dict(
omega_rpm=7500, # spindle speed (rpm)
b_mm=5.0, # axial depth of cut (mm)
Nt=4, # number of teeth
ft_mm=0.15, # feed per tooth (mm)
phis_deg=0.0, # cut start angle (deg)
phie_deg=180.0, # cut exit angle (deg): slotting
fn_x=500.0, k_x=8e6, zeta_x=0.02, # x direction dynamics
fn_y=500.0, k_y=8e6, zeta_y=0.02, # y direction dynamics
kt=695e6, kn=281e6, # cutting force coefficients (N/m^2)
steps_rev=652, revolutions=20)
if PRESET in PRESETS:
params.update(PRESETS[PRESET])
sim = mill_tds(**params)
M, verdict = chatter_metric(sim)
f_tooth = params['omega_rpm'] * params['Nt'] / 60
print(f"Case: {params['omega_rpm']} rpm, {params['b_mm']} mm, "
f"f_tooth = {f_tooth:.0f} Hz.")
print(f"Periodic sampling metric M = {M:.2f} µm -> {verdict}.")
Case: 7500 rpm, 5.0 mm, f_tooth = 500 Hz. Periodic sampling metric M = 66.61 µm -> chatter.
# --- Fig. C.6: the four Example 4.11 cases ---
case_results = []
for lab in ('1', '2', '3', '4'):
p4 = dict(params)
p4.update(PRESETS[lab])
p4['revolutions'] = 30 # a little longer, so the sampled half is easy to see
s4 = mill_tds(**p4)
M4, v4 = chatter_metric(s4)
case_results.append((lab, p4, s4, M4, v4))
fig = make_subplots(rows=2, cols=2, horizontal_spacing=0.12,
vertical_spacing=0.17,
subplot_titles=[f"case {lab}: {p4['omega_rpm']:,} rpm, "
f"{p4['b_mm']} mm ({v4}, M = {M4:.2f} µm)"
for lab, p4, s4, M4, v4 in case_results])
for i_c, (lab, p4, s4, M4, v4) in enumerate(case_results):
r_, c_ = i_c // 2 + 1, i_c % 2 + 1
tt, xx = decimate_for_plot(s4['t'], s4['x'] * 1e6, max_points=3000)
fig.add_trace(go.Scatter(x=tt, y=xx, mode='lines',
line=dict(color='#1f77b4', width=0.8),
showlegend=False), row=r_, col=c_)
s0 = len(s4['t']) // 2
si = np.arange(s0, len(s4['t']), s4['steps_tooth'])
si = si[s4['t'][si] <= tt[-1]] # samples never extend past the line
fig.add_trace(go.Scatter(x=s4['t'][si], y=s4['x'][si] * 1e6,
mode='markers',
marker=dict(color='#d62728', size=4),
name='once-per-tooth samples',
showlegend=(i_c == 0)), row=r_, col=c_)
fig.update_xaxes(title_text='t (s)', row=r_, col=c_)
fig.update_yaxes(title_text='x (µm)', row=r_, col=c_)
fig.update_layout(title='Fig. C.6 — Time domain simulation of the four Example 4.11 cases',
template='simple_white', height=560,
legend=dict(x=1.0, y=1.12, xanchor='right', orientation='h'),
margin=dict(l=60, r=20, t=95, b=50))
fig.show()
The four panels replicate Example 4.11 of [1]: cases 1 and 3 are stable and 2 and 4 are chatter, where the once-per-tooth samples repeat for stable cutting and vary from one to the next for chatter. Compare each position against the boundaries of Fig. C.4. This matches the zero-order boundary and resolves the Sect. C.4 disagreement in its favor. Because this Part is about sound, the cell below provides the audible result for all four cases: the simulation time is selected so that each is long enough to hear, the displacement record is converted to audio, and the four clips share one loudness scale so they represent the relative vibration levels. The stable cuts are nearly silent, since the constant summed force of Sect. C.1 for this idealized slotting cut excited almost no vibration; a real stable cut still carries runout and ambient sound, as Sect. C.2 showed. The cuts that exhibit chatter are loud with different frequency content. The difference is clear.
# --- Hearing the simulation: the four Example 4.11 cases on one loudness scale ---
def sonify(omega_rpm, b_mm, seconds=1.2, **kw):
revs = int(np.ceil(seconds * omega_rpm / 60))
s = mill_tds(omega_rpm=omega_rpm, b_mm=b_mm, revolutions=revs, **kw)
a = resample_to_audio(s['x'] - s['x'].mean(), s['fs'])
M, verdict = chatter_metric(s)
return fade(a), verdict # unscaled; callers set loudness
clips = []
for lab in ('1', '2', '3', '4'):
p4 = PRESETS[lab]
a4, v4 = sonify(p4['omega_rpm'], p4['b_mm'])
clips.append((lab, p4, a4, v4))
peak = max(np.max(np.abs(a4)) for _, _, a4, _ in clips)
for lab, p4, a4, v4 in clips:
print(f"case {lab} ({p4['omega_rpm']:,} rpm, {p4['b_mm']} mm): "
f"simulation verdict {v4}")
display(Audio(0.8 * a4 / (peak + 1e-30), rate=AUDIO_RATE, normalize=False))
case 1 (7,500 rpm, 3.0 mm): simulation verdict stable
case 2 (7,500 rpm, 5.0 mm): simulation verdict chatter
case 3 (5,000 rpm, 0.1 mm): simulation verdict stable
case 4 (5,000 rpm, 0.5 mm): simulation verdict chatter
▸ Key terms: time domain simulation, once-per-tooth sampling
C.6 Mapping Stability with the Simulation¶
The analytical stability maps from Sect. C.4 and the simulation from Sect. C.5 answer the same question by different routes. The simulation can either study a single set of operating parameters or can be used to generate the entire map. To generate a stability map, the simulation is executed over a grid of spindle speeds and depths, applying the periodic sampling metric to each combination, and marking each grid point as stable or chatter using the one micrometer threshold. The following cell calculates a stability map. The grid density and number of revolutions per grid point are input parameters; the defaults trade resolution for a runtime of about a minute, with progress presented as it proceeds.
# --- Build a stability map from time domain simulations. Adjust and rerun. ---
import time as _time
n_omega = 20 # grid points in spindle speed
n_b = 13 # grid points in axial depth
revs_map = 25 # revolutions per simulation (transient plus judgement)
print(f'Grid: {n_omega} speeds by {n_b} depths, {revs_map} revolutions each.')
Grid: 20 speeds by 13 depths, 25 revolutions each.
# --- Stability map sweep ---
om_grid = np.linspace(3000, 15000, n_omega)
b_grid = np.linspace(0.4, 6.5, n_b)
Mmap = np.zeros((n_b, n_omega))
t0 = _time.perf_counter()
for i_o, om in enumerate(om_grid):
for i_b, bb in enumerate(b_grid):
s = mill_tds(omega_rpm=om, b_mm=bb, revolutions=revs_map)
M, verdict = chatter_metric(s)
Mmap[i_b, i_o] = M
print(f' column {i_o + 1:2d}/{n_omega} ({om:,.0f} rpm) done, '
f'{_time.perf_counter() - t0:5.1f} s elapsed')
print(f'{n_omega * n_b} simulations in {_time.perf_counter() - t0:.1f} s.')
column 1/20 (3,000 rpm) done, 1.6 s elapsed column 2/20 (3,632 rpm) done, 3.4 s elapsed column 3/20 (4,263 rpm) done, 5.0 s elapsed column 4/20 (4,895 rpm) done, 6.6 s elapsed column 5/20 (5,526 rpm) done, 8.2 s elapsed column 6/20 (6,158 rpm) done, 9.8 s elapsed column 7/20 (6,789 rpm) done, 11.2 s elapsed column 8/20 (7,421 rpm) done, 12.9 s elapsed column 9/20 (8,053 rpm) done, 14.4 s elapsed column 10/20 (8,684 rpm) done, 16.0 s elapsed column 11/20 (9,316 rpm) done, 17.6 s elapsed column 12/20 (9,947 rpm) done, 19.1 s elapsed column 13/20 (10,579 rpm) done, 20.7 s elapsed column 14/20 (11,211 rpm) done, 22.2 s elapsed column 15/20 (11,842 rpm) done, 23.7 s elapsed column 16/20 (12,474 rpm) done, 25.3 s elapsed column 17/20 (13,105 rpm) done, 26.8 s elapsed column 18/20 (13,737 rpm) done, 28.4 s elapsed column 19/20 (14,368 rpm) done, 29.8 s elapsed column 20/20 (15,000 rpm) done, 31.5 s elapsed 260 simulations in 31.5 s.
# --- Fig. C.7: the stability map ---
fig = go.Figure()
OMg, BBg = np.meshgrid(om_grid, b_grid)
stab = Mmap < 1.0
fig.add_trace(go.Scatter(x=OMg[stab], y=BBg[stab], mode='markers',
marker=dict(symbol='circle-open', size=8,
color='#2ca02c', line=dict(width=1.5)),
name='stable (simulation)', customdata=Mmap[stab],
hovertemplate='%{x:.0f} rpm, %{y:.2f} mm<br>'
'M = %{customdata:.2f} µm<extra></extra>'))
fig.add_trace(go.Scatter(x=OMg[~stab], y=BBg[~stab], mode='markers',
marker=dict(symbol='x', size=7, color='#d62728'),
name='chatter (simulation)', customdata=Mmap[~stab],
hovertemplate='%{x:.0f} rpm, %{y:.2f} mm<br>'
'M = %{customdata:.2f} µm<extra></extra>'))
for i, (om, bl, fc) in enumerate(tl):
fig.add_trace(go.Scatter(x=om, y=bl, mode='lines',
line=dict(color='#1f77b4'),
name='average tooth angle', showlegend=(i == 0)))
fig.add_trace(go.Scatter(x=zo_om, y=zo_bl, mode='lines',
line=dict(color='#ff7f0e', dash='dash'),
name='Fourier series (zero-order)'))
for label, om, b, verdict in CASES:
fig.add_trace(go.Scatter(x=[om], y=[b], mode='markers+text',
marker=dict(color='black', size=9, symbol='diamond'),
text=[label], textposition='top center', showlegend=False))
fig.update_layout(title='Fig. C.7 — Stability map from the time domain simulation',
template='simple_white', height=500,
xaxis=dict(title='Spindle Speed (rpm)', range=[3000, 15000]),
yaxis=dict(title='Axial Depth of Cut, b (mm)', range=[0, 6.5]),
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=60, b=50))
fig.show()
The crosses identify unstable results (this is the stability map's chatter region) and the open circles identify stable results (this is the stable region). The zero-order boundary agrees with the time domain simulation results, while the average tooth angle curve overestimates the stable region for slotting based on the reason Sect. C.4 identified: with $\mu_y = 0$ it cannot see the y direction dynamics that the simulation retains. The comparison helps us decide which approach to implement. The analytical solution time is milliseconds, but is based on averaging assumptions whose accuracies vary with the cut geometry. The time domain map solution time is about a minute, but does not depend on the same simplifying assumptions and retains everything, including the nonlinearity of teeth leaving the cut.
The periodic sampling metric that is used to automatically label each grid point can also be applied to real cuts in practice. We can sample the vibration signal once per tooth passage and quantify at what level the samples repeat. A periodic, stable cut repeats itself every tooth period, so the samples are nominally the same. Chatter, on the other hand, contains a new frequency that is not a multiple of the tooth passing frequency in general, so the samples vary from one tooth passage to the next. In the metric calculation, the samples are differenced and averaged. The one micrometer threshold used here reflects the difference demonstrated in Sect. 4.4.5 of [1]: a near zero value for stable cuts and tens of micrometers for chatter.
▸ Key terms: stability map
C.7 Hearing Chatter¶
Analyzing the frequency content for a cut is a great basis for identifying chatter. A stable cut's spectrum contains only the families of Sect. C.2: the tooth passing frequency and its harmonics, the runout family, and the idle baseline. Chatter adds a peak near a structural natural frequency that belongs to none of them in general [1]. The following cell analyzes the audible case 2 cut from Sect. C.5 in the three views, with the tooth passing family marked; the chatter peak is seen beside the expected peaks, close to them but separate.
# --- Fig. C.8: the chatter cut in three views, tooth passing family marked ---
audio_ch, _ = sonify(7500, 5.0)
audio_ch = 0.8 * audio_ch / (np.max(np.abs(audio_ch)) + 1e-30)
display(Audio(audio_ch, rate=AUDIO_RATE, normalize=False))
f_tooth = 500.0
f_sp, a_sp = compute_spectrum(audio_ch, AUDIO_RATE, fmax=2600)
t_g, f_g, S = compute_spectrogram(audio_ch, AUDIO_RATE, fmax=2600)
tt, xx = decimate_for_plot(np.arange(len(audio_ch)) / AUDIO_RATE, audio_ch,
max_points=6000)
fig = make_subplots(rows=1, cols=3, horizontal_spacing=0.10,
column_titles=('time domain', 'spectrum', 'spectrogram'))
fig.add_trace(go.Scatter(x=tt, y=xx, mode='lines',
line=dict(color='#1f77b4', width=0.7),
showlegend=False), row=1, col=1)
fig.add_trace(go.Scatter(x=f_sp, y=a_sp, mode='lines',
line=dict(color='#1f77b4'), showlegend=False), row=1, col=2)
for m in range(1, 6):
fig.add_shape(type='line', x0=m * f_tooth, x1=m * f_tooth, y0=0, y1=1.02,
xref='x2', yref='y2',
line=dict(color='black', width=1, dash='dot'), opacity=0.5)
fig.add_trace(go.Heatmap(x=t_g, y=f_g, z=S, colorscale='Inferno',
showscale=False), row=1, col=3)
fig.update_xaxes(title_text='t (s)', row=1, col=1)
fig.update_xaxes(title_text='Frequency (Hz)', row=1, col=2)
fig.update_xaxes(title_text='t (s)', row=1, col=3)
fig.update_yaxes(title_text='Amplitude', row=1, col=1)
fig.update_yaxes(title_text='Amplitude (normalized)', range=[0, 1.05], row=1, col=2)
fig.update_yaxes(title_text='Frequency (Hz)', range=[0, 2600], row=1, col=3)
fig.update_layout(title='Fig. C.8 — The chatter cut in three views',
template='simple_white', height=380,
margin=dict(l=60, r=20, t=70, b=50))
fig.show()
# identify the chatter frequency: the tallest peak away from the tooth family
mask = np.ones_like(f_sp, dtype=bool)
for m in range(1, 8):
mask &= np.abs(f_sp - m * f_tooth) > 10
mask &= f_sp > 100
fc_detected = f_sp[mask][np.argmax(a_sp[mask])]
print(f'Dotted lines: the tooth passing family at multiples of {f_tooth:.0f} Hz.')
print(f'Tallest peak outside the family: {fc_detected:.0f} Hz, the chatter frequency, '
f'near the {500:.0f} Hz natural frequency and on no family line.')
Dotted lines: the tooth passing family at multiples of 500 Hz. Tallest peak outside the family: 487 Hz, the chatter frequency, near the 500 Hz natural frequency and on no family line.
One caution before trusting this analysis on a real machine. Runout populates the spectrum at multiples of the once-per-revolution frequency, and a runout harmonic can land close to a plausible chatter frequency; Example 6.4 of [1] shows a case where the third runout harmonic and the chatter frequency coincide and cannot be separated by frequency alone. The once-per-tooth sampling metric is the robust companion test, since runout content is periodic with the rotation and does not affect the metric value, while chatter does.
This detection chain, spectrum, family identifications, and once-per-tooth sampling, runs on live audio as well: a microphone supplies the signal, the programmed spindle speed and tooth count identify the families, and any persistent peak outside them is flagged and its frequency is identified. That frequency is the input to the next section.
▸ Key terms: chatter frequency
C.8 Choosing a Better Spindle Speed¶
The chatter frequency represents key information: it identifies where the structure wants to vibrate, and Eq. 4.24 turns it into a prescription for a better spindle speed. Setting $\varepsilon = 2\pi$, the best speeds place the tooth passing frequency in step with the chatter frequency:
$$\Omega_{new} = \frac{60 \, f_c}{(N + 1) \, N_t}$$
for lobe numbers $N = 0, 1, 2, \ldots$, which is the spindle speed regulation procedure from Sect. 6.1 of [1]. The new speed drives the tooth-to-tooth phase toward the favorable nesting of Fig. C.3, which places the operating point at one of the peaks of the stable region in the stability map of Fig. C.4. Because the chatter frequency shifts slightly at the new speed, the procedure may require an iteration or two to converge, as Example 6.1 of [1] demonstrates. The following cell completes the full analysis on a cut that chatters: simulate, detect $f_c$ with the Sect. C.7 rule, compute the candidate speeds, pick the nearest, and repeat until the once-per-tooth metric reports a stable cut.
# --- Speed selection: starting point. Adjust and rerun. ---
omega_0 = 6000.0 # starting spindle speed (rpm)
b_run = 2.5 # axial depth (mm), held fixed throughout
Nt_run = 4
fn_ref = 500.0 # the structure's natural frequency (Hz), for the search band
print(f'Starting point: {omega_0:,.0f} rpm at b = {b_run} mm.')
Starting point: 6,000 rpm at b = 2.5 mm.
# --- The listen, identify, adjust loop ---
def detect_fc(omega_rpm, b_mm):
revs = int(np.ceil(0.9 * omega_rpm / 60)) # simulate about 0.9 s of cutting
s = mill_tds(omega_rpm=omega_rpm, b_mm=b_mm, revolutions=revs)
M, verdict = chatter_metric(s)
half = len(s['x']) // 2
f_sp, a_sp = compute_spectrum(s['x'][half:], s['fs'], fmax=2500)
f_tp = omega_rpm * Nt_run / 60
mask = f_sp > 0.95 * fn_ref # chatter sits near and above fn
for m in range(1, 9): # remove the tooth passing family
mask &= np.abs(f_sp - m * f_tp) > max(8.0, 0.015 * f_tp)
fc = f_sp[mask][np.argmax(a_sp[mask])]
return M, verdict, fc
omega = omega_0
history = []
for it in range(5):
M, verdict, fc = detect_fc(omega, b_run)
history.append((omega, verdict))
print(f'iteration {it}: {omega:7,.0f} rpm -> metric M = {M:8.2f} µm ({verdict})'
+ ('' if verdict == 'stable' else f', chatter frequency {fc:6.1f} Hz'))
if verdict == 'stable':
break
candidates = [60 * fc / ((N + 1) * Nt_run) for N in range(0, 4)]
omega = min(candidates, key=lambda c: abs(np.log(c / omega)))
print(f' best speeds {", ".join(f"{c:,.0f}" for c in candidates)} rpm; '
f'selecting {omega:,.0f} rpm')
print()
print(f'Converged: {history[-1][0]:,.0f} rpm at b = {b_run} mm, '
f'from a starting point of {omega_0:,.0f} rpm.')
print('Listen to the first and last cuts of the loop:')
a0, _ = sonify(history[0][0], b_run)
a1, _ = sonify(history[-1][0], b_run)
peak = max(np.max(np.abs(a0)), np.max(np.abs(a1)))
a0 = 0.8 * a0 / (peak + 1e-30)
a1 = 0.8 * a1 / (peak + 1e-30)
print(f'before, {history[0][0]:,.0f} rpm:')
display(Audio(a0, rate=AUDIO_RATE, normalize=False))
print(f'after, {history[-1][0]:,.0f} rpm:')
display(Audio(a1, rate=AUDIO_RATE, normalize=False))
iteration 0: 6,000 rpm -> metric M = 661.50 µm (chatter), chatter frequency 506.7 Hz
best speeds 7,600, 3,800, 2,533, 1,900 rpm; selecting 7,600 rpm
iteration 1: 7,600 rpm -> metric M = 0.00 µm (stable)
Converged: 7,600 rpm at b = 2.5 mm, from a starting point of 6,000 rpm.
Listen to the first and last cuts of the loop:
before, 6,000 rpm:
after, 7,600 rpm:
The result is a spindle speed just above 7500 rpm, beside the stable pocket that case 1 of Example 4.11 occupies, where the tooth passing frequency falls in step with the structure's preferred frequency. With a larger axial depth, the first adjusted speed can itself chatter at a new chatter frequency, and the loop simply repeats. Example 6.1 of [1] shows this behavior before converging. The audio pair summarizes the change: self-excited vibration at the first speed and forced vibration at the second, with nothing changed but the spindle speed. On the map of Fig. C.7, the operating point has moved horizontally from the chatter region into a stable pocket, and Fig. C.4 explains why the pocket is there.
▸ Key terms: best speeds, spindle speed regulation
Summary¶
- Milling removes material with a rotating, toothed tool; the chip thickness varies as $h = f_t \sin\phi$ (Eq. 4.2) over the angular range set by the radial immersion, and the cutting force follows it based on the coefficients $k_t$ and $k_n$ (Eqs. 4.7 to 4.10).
- The cutting force is periodic with its fundamental frequency $f_{tooth} = \Omega N_t / 60$ (Eq. 4.14). In practice, a stable cut's spectrum contains the tooth passing frequency family, the once-per-revolution runout frequency family, and the idle baseline. Every frequency can be identified and named.
- Regeneration closes a feedback loop: each tooth cuts the surface the previous tooth left behind (Eqs. 4.21 and 4.22), and the tooth-to-tooth phase determines whether vibration is harmlessly copied or pumped into chatter.
- The stability map gives the limiting depth as a function of spindle speed. Two analytical solutions are available to compute it: the average tooth angle approach (Eqs. 4.23, 4.24, 4.27) and the Fourier series zero-order approach (Eqs. 4.44, 4.45). For slotting, they disagree, because $\mu_y = 0$ removes the y direction from the average tooth angle result, and the zero-order boundary is the one the simulation confirms.
- The time domain simulation of Sect. 4.4 keeps the full geometry and the out-of-cut nonlinearity; its four Example 4.11 cases reproduce the Example 4.11 results of [1], and repeating it over a grid of points produces a stability map.
- Once-per-tooth sampling separates stable cuts from chatter: stable cuts have repeating samples, while cuts that chatter do not. The stability metric based on these samples is robust to runout content that can confuse frequency-based chatter identification.
- The chatter frequency can be used to select better spindle speeds through $\Omega_{new} = 60 f_c / ((N+1) N_t)$. Iterating the listen, identify, and adjust loop identifies a stable pocket, provided one is available at the selected axial depth of cut.
The series ends as intended: the stability map is computed by two analytical algorithms, checked by a time domain simulation, and navigated by listening. The same procedure applies at a machine: record the cut sound, identify and ignore the tooth passing family, find any remaining peak and identify it as the chatter frequency, and adjust the spindle speed to place a tooth passing harmonic on it.
Exercises¶
1. Using the parameter panel of Sect. C.5, change the feed per tooth from 0.15 to 0.05 mm/tooth for case 2. (a) Does the cut remain in chatter? (b) Explain, using Eqs. 4.7 and 4.22, why the feed per tooth scales the forces and vibrations but does not decide stability.
2. Run case 2 with steps_rev set to 200, 652, and 2000.
(a) How do the verdict and the displacement histories change?
(b) What does this suggest as a working rule for choosing the number of steps per revolution?
3. In Fig. C.5, the whole diagram scales with the damping ratio. (a) Using Eq. 4.23 and the Part A result that the FRF peak scales as $1/(2 k \zeta)$, show that the critical (minimum) limiting depth is proportional to $\zeta$ for small $\zeta$. (b) A tooling change doubles the damping ratio. What happens to the deepest stable slotting cut at the worst spindle speed?
4. Compute both analytical boundaries for a 25 percent radial immersion up milling cut ($\phi_s = 0$, $\phi_e = 60$ deg) using the Sect. C.4 cells. (a) How do the two solutions compare with each other at this low immersion, and with the slotting result? (b) Chap. 4 of [1] cautions that averaging degrades as the force becomes impulse-like. Which figure of Sect. C.1 shows why?
5. Rebuild the Sect. C.6 map with zeta_x = zeta_y = 0.04 (edit mill_tds's defaults or pass the values through).
(a) Compare the chatter region with Fig. C.7.
(b) Do the analytical boundaries recomputed at the new damping still trace its edge?
6. Retune the Sect. C.8 loop for a different structure, the system of Example 6.1 of [1] made symmetric: in the visible starting point cell set omega_0 = 10000, b_run = 1.5, and fn_ref = 800, and in the collapsed loop cell add fn_x=800, fn_y=800, k_x=5e6, k_y=5e6, zeta_x=0.01, zeta_y=0.01 to the mill_tds call inside detect_fc.
(a) Before running, compute the best speeds $60 f_n / ((N+1) N_t)$ for $N = 0$ and $N = 1$.
(b) Run the loop. Where does it land relative to your predictions, and how many iterations does it take?
Appendix: Utility Functions¶
The shared helper functions, defined once in the collapsed cell at the top of the notebook, are carried over from Parts A and B: fade, compute_spectrum, compute_spectrogram, and decimate_for_plot, plus one addition, resample_to_audio(x, fs_in, fs_out), which interpolates a simulated record onto the audio sample grid for playback. The stability solutions tlusty_lobes and zoa_lobes and the simulation mill_tds with its companion chatter_metric are defined in full view in Sects. C.4 and C.5, because they are the subject matter of this Part rather than plumbing.
References¶
[1] Schmitz, T. L., and Smith, K. S. (2019). Machining Dynamics: Frequency Response to Improved Productivity, 2nd ed. Springer. DOI: 10.1007/978-3-319-93707-6
[2] Altintas, Y., and Budak, E. (1995). Analytical prediction of stability lobes in milling. CIRP Annals, 44(1), 357 to 362.
[3] Taylor, F. W. (1907). On the art of cutting metals. Transactions of the ASME, 28, 31 to 350.
[4] Schmitz, T. L., and Gomez, M. F. Machining Dynamics: Jupyter Notebook Edition (in preparation).