Research on the dynamic response of magnetic levitation milling electric spindle considering the variable stiffness characteristics of the spindle-holder-tool joint system

Abstract

The time-varying contact stiffness at the spindle-holder-tool joint in a magnetic levitation milling spindle significantly affects system dynamics during high-speed machining. To address this issue, a dynamic model of the rotor system is developed by integrating Timoshenko beam theory, rotor dynamics, and joint contact theory. The joint interface is modeled using distributed spring-damping elements, and the influence of centrifugal expansion on contact stiffness degradation is analyzed. The model incorporates the effects of magnetic bearing support stiffness, gyroscopic moments, and milling force excitation. The dynamic response is obtained using time-domain numerical integration. Results indicate that increasing rotational speed reduces joint contact stiffness, alters mode shapes, and amplifies tool tip vibration amplitudes, thereby diminishing the system’s vibration resistance and increasing the risk of instability. This study provides a theoretical foundation for improving the stability and dynamic performance of high-speed magnetic levitation spindles.

Keywords

magnetic levitation electric spindle spindle-holder-tool joint time-varying stiffness dynamic response centrifugal expansion

1. Introduction

Magnetic levitation milling spindles are increasingly employed in ultra-precision machining, such as aerospace manufacturing, due to their contact-free support, lubrication-free, and high rotational speed capabilities. The dynamic behavior of the spindle system directly influences machining accuracy and operational stability. Among various factors, the connection performance of the spindle-holder-tool joint plays a critical role in determining the overall dynamic characteristics (Ertürk et al., 2007a, 2007b; Salahshoor and Ahmadian, 2009). During high-speed machining, the joint interface is susceptible to stiffness degradation and dynamic instability (Gagnol et al., 2011). It has been reported that 25-50% of the total tool tip deformation originates from joint interfaces (Shamine and Shin, 2000). Therefore, accurate modeling of the joint’s dynamic characteristics and investigation of its vibration response are essential for improving machining stability and precision.

Considerable research has been devoted to the dynamic modeling and parameter identification of spindle-holder-tool systems. Experimental techniques, particularly impact hammer testing, are widely used to identify dynamic parameters due to their simplicity and efficiency (Iglesias et al., 2022). In milling dynamics modeling, this approach is frequently applied to determine modal parameters such as stiffness, damping, and modal mass at the tool tip (Li et al., 2020). For modeling contact stiffness at the spindle-holder interface, Namazi et al. (2007) proposed a Timoshenko beam-based method using uniformly distributed translational and rotational springs to quantify radial and rotational stiffness and identified spring parameters by minimizing errors in spindle assembly experiments. Xu et al. (2014) incorporated clamping and centrifugal forces into their model and identified axial and radial joint stiffness experimentally. Similarly, considering flexible deformation effects, Zhao et al. (2024) investigated the contact characteristics of mechanical transmission systems, highlighting the significant influence of structural compliance on interface stiffness. In addition, frequency response function (FRF) identification methods have been widely used to analyze the dynamic behavior of spindle interfaces. For instance, Özşahin et al. (2009) used tool tip FRF to evaluate the influence of spindle-holder interface stiffness on vibration amplitude. Du et al. (2015) proposed an ultrasonic method to examine the effects of nominal contact area, contact pressure distribution, and interface stiffness on tool tip vibration.

In addition to experimental approaches, analytical models have been developed to predict joint dynamics (Sakamoto et al., 2012). Ertürk et al. (2006) presented an analytical algorithm based on Timoshenko beam theory, calculating point interface contact stiffness through electro-nano coupling and structural modification. Zhao et al. (2013) developed a dynamic model of the spindle rotor system considering contact at tool-jacket, jacket-holder, and holder-spindle interfaces. Using Euler-Bernoulli beam theory, Luo et al. (2022) segmented the spindle-holder system into three substructures and applied the transfer matrix method together with multi-point response coupling to predict the end FRF. Kiran (2022) treated the holder-tool joint as a flexible component and proposed a frequency response coupling analysis accounting for all translational and rotational mechanical characteristics. Liu et al. (2019) and Cui and Su (2020) developed analytical models of contact stiffness at the holder-tool interface using fractal and multi-scale contact theories and integrated them with finite element models to predict chatter stability. Miao et al. (2023) employed a distributed nonlinear spring-damper model to describe the contact behavior at spindle-holder-tool joints, deriving a stiffness matrix dependent on vibration displacement via fractal theory and the slicing method. Furthermore, nonlinear dynamic behavior and vibration characteristics of mechanical systems under complex operating conditions have been extensively studied, demonstrating that contact nonlinearity plays a crucial role in system dynamic response (Zhao et al., 2025a). In addition, related studies on contact modeling and system reliability provide further insights into the coupling between interface stiffness and overall dynamic performance (Zhao et al., 2023, 2025b).

Although existing studies have provided valuable insights into spindle dynamics, several critical issues remain insufficiently addressed, particularly for magnetic levitation milling spindles.

First, most previous research has focused on conventional spindle systems supported by rolling bearings, whose boundary conditions and support characteristics differ significantly from those of magnetic levitation systems. The absence of mechanical contact in magnetic bearings leads to fundamentally different dynamic behaviors, which are not adequately captured by existing models.

Second, the evolution mechanism of joint contact stiffness under high-speed operating conditions has not been clearly quantified. In particular, limited theoretical work has been conducted on the stiffness variation of the spindle-holder-tool joint under the combined effects of centrifugal forces and drawbar preload variation. As a result, the mechanism of stiffness degradation at high rotational speeds remains unclear.

Third, the relationship between joint stiffness variation and system dynamic response, especially under varying rotational speeds, has not been systematically investigated. This limitation restricts the accurate prediction of vibration characteristics and the effective optimization of system performance.

These issues make it difficult to accurately predict the dynamic behavior of magnetic levitation spindle systems under high-speed operating conditions. To address these issues, this study develops a comprehensive dynamic model of a magnetic levitation milling spindle system that incorporates the variable stiffness characteristics of the spindle-holder-tool joint. The main contributions of this work are summarized as follows:

(1) A coupled dynamic model of a magnetic levitation spindle system considering the spindle-holder-tool joint characteristics is established.

(2) The stiffness variation of the joint induced by centrifugal effects is systematically analyzed based on elastic theory.

(3) The influence of joint stiffness variation on the dynamic response of the system is quantitatively evaluated under different rotational speeds.

The results of this study provide a theoretical basis for the design optimization and stability improvement of high-speed magnetic levitation milling spindles.

2. Dynamic model of magnetic levitation milling electric spindle

Figure 1(a) illustrates the structure of the magnetic levitation milling spindle. The main components include the spindle, front and rear radial magnetic bearings, motor rotor, holder, collet (jacket), lock nut, and cutter. The inner conical surface of the spindle and the outer conical surface of the holder are held together by the drawbar pulling force. The holder and collet are fixed onto the holder via the lock nut. The motor supplies driving torque to the spindle, while the radial magnetic bearings provide electromagnetic support to achieve stable suspension.

Figure 1.

Spindle-holder-tool system.

To analyze the contact characteristics of the joint, the finite element method is used to discretize the spindle-holder-tool system into interconnected “mass-spring-damping” elements, as shown in Figure 1(b).

In Figure 1(b), in order to quantify the elastic deformation of the joint under milling loads, the tapered joint is divided into stepped shaft segments, and the straight-shank end mill is simplified as a uniform circular cross-section beam element. Axial displacements caused by radial deformations are neglected. The system is ultimately represented by 2-node, 8-degree-of-freedom Timoshenko beam elements, each node having translational degrees of freedom along the x- and y-axes and rotational degrees of freedom about the x- and y-axes. The spindle-holder and holder-tool joints are modeled using distributed nonlinear spring-damper elements, while the magnetic bearing support is introduced with equivalent electromagnetic stiffness. The spindle-magnetic bearing interface is modeled using a 4-DOF nonlinear stiffness formulation.

2.1. Dynamic model of the spindle-holder taper joint

When the spindle and the holder are combined, they are subject to the pulling force of the drawbar. Therefore, when it rotates, there are mechanical characteristics such as contact surface positive pressure, torsional moment, and conical surface friction force between its joint surfaces. When modeling, based on the finite element method, each conical unit is first simplified into a Timoshenko beam with an equivalent diameter of its median diameter, and then based on the force characteristics of each beam element node, it is simplified into a spring-damping model, as shown in Figure 2.

Figure 2.

Distributed spring-damping model of the spindle-holder joint.

According to Figure 2, when analyzing the stress characteristics of each Timoshenko beam unit node on the conical surface of the spindle-tool holder joint, if the nodes on the spindle are represented by and the nodes on the holder are represented by, then the radial force, axial force, torque, friction, etc., experienced by each node are represented by springs, torsion springs, and damping in turn. The equivalent physical model is shown in Figure 3.

Figure 3.

Distributed spring-damping node coupling model.

However, during the operation of the magnetic levitation electric spindle, the spindle-holder joint will produce corresponding displacements under the action of force and torque, which is reflected in Figure 3 as corresponding elastic deformation, which in turn triggers the mutual coupling of forces and moments between each node.

In order to describe the coupling relationship between the force and torque of each joint surface during rotation, when the $i, i + 1, \dots, i + n$ node of the spindle and the $j, j + 1, \dots, j + n$ node of the tool holder are in relative motion, the coupling force $Δ F_{p} (p = i, i + 1, \dots, i + n)$ and coupling moment $Δ M_{p} (p = i, i + 1, \dots, i + n)$ acting on the spindle node, and the coupling force $Δ F_{q} (q = j, j + 1, \dots, j + n)$ and coupling moment $Δ M_{q} (q = i, i + 1, \dots, i + n)$ on the holder node, their expressions are:

{\begin{cases} Δ F_{p} = k_{c p} (x_{i} - x_{j}) + c_{c p} (\dot{x_{i}} - \dot{x_{j}}) \\ Δ F_{q} = k_{c p} (x_{j} - x_{i}) + c_{c p} (\dot{x_{j}} - \dot{x_{i}}) \\ Δ M_{p} = k_{c p} (θ_{i} - θ_{j}) + c_{c p} (\dot{θ_{i}} - \dot{θ_{j}}) \\ Δ M_{q} = k_{c p} (θ_{j} - θ_{i}) + c_{c p} (\dot{θ_{j}} - \dot{θ_{i}}) \end{cases}

(1)

In the formula,

x_{i}, x_{j}

are the displacement vectors at the spindle-tool holder joint;

θ_{i}, θ_{j}

are the angular displacement vectors at the spindle-holder joint; and

k_{c p}, c_{c p}

are the coupling stiffness and damping of each node between the spindle and the holder joint surface, their expressions are:

\begin{array}{l} k_{c p} = (\begin{array}{l} k_{i v} \\ k_{i w} \\ k_{i θ v} \\ k_{i θ w} \end{array}) \end{array}

(2)

c_{c p} = (\begin{array}{l} c_{i v} \\ c_{i w} \\ c_{i θ v} \\ c_{i θ w} \end{array})

(3)

in the formula,

k_{i v}

k_{i θ v}

k_{i w}

k_{i θ w}

, respectively, represent the translational stiffness and rotational stiffness of the

i^{t h}

unit in the

v

w

directions;

c_{i v}

c_{i θ v}

c_{i w}

c_{i θ w}

, respectively, represent the translational damping and rotational damping of the

i^{t h}

unit in the

v

w

directions.

The overall stiffness matrix and overall damping matrix of the spindle-holder joint can be expressed as:

\begin{array}{l} K_{c} = (\begin{array}{l} k_{c p} & 0 & 0 & - k_{c p} \\ 0 & 0 \\ 0 & 0 \\ - k_{c p} & 0 & 0 & k_{c p} \end{array}) \end{array}

(4)

C_{c} = (\begin{array}{l} c_{c p} & 0 & 0 & - c_{c p} \\ 0 & 0 \\ 0 & 0 \\ - c_{c p} & 0 & 0 & c_{c p} \end{array})

(5)

2.1.1. Static contact stiffness model of spindle-holder joint

As shown in Figure 4, under static conditions, the pretightening force $f_{1}$ of the tool holder causes the surface of the inner taper hole in the spindle-tool holder joint to form stable contact.

Figure 4.

Stress analysis diagram of the spindle-tool holder joint.

Assuming that the conical contact pressure between the spindle and the holder is uniformly distributed along the axial direction, then according to the conical contact mechanics theory, the contact pressure per unit area $P_{n 1}$ can be obtained. The expressions between each force and the contact surface pressure are as follows:

{\begin{cases} f_{1} = f_{1} \sin θ_{1} + f_{μ 1} \cos θ_{1} \\ A_{1} = π D_{m} \frac{L_{1}}{\cos θ_{1}} \\ P_{n 1} = \frac{f_{n_{1}}}{A_{1}} = \frac{f_{1}}{π D_{m} L_{1} (μ_{1} + \tan θ_{1})} \\ D_{m} = \frac{(D_{1} + d_{1})}{2} \end{cases}

(6)

in the formula,

f_{μ 1}

is the friction force of the outer conical surface of the tool holder;

f_{n 1}

is the positive pressure of the outer conical surface;

μ_{1}

is the friction coefficient of the outer conical surface; and

θ_{1}

is the conical angle of the spindle-tool holder joint surface, which is 8.297°.

A_{1}

is the spindle-tool holder contact taper surface area;

D_{m}

is the tool holder-spindle taper surface matching diameter; and

L_{1}

is the conical surface matching distance.

According to the relationship between the contact pressure and contact stiffness of the joint in the literature (Zhao et al., 2013), and using the finite element method to discretize the joint surface into $N$ segments, the median diameter of each segment is $D_{i m}$ , the normal stiffness $k_{i n}$ , tangential stiffness $k_{i τ}$ , and angular stiffness $k_{i θ}$ of each contact unit of the joint are obtained. The calculation formula is:

{\begin{cases} k_{i n} = α {p_{i}}^{β} A_{i} = \frac{α {p_{i}}^{β} π D_{i m} L_{i}}{N \cos θ} \\ k_{i τ} = 0.5 k_{i n} \\ k_{i θ} = \frac{α {p_{i}}^{β} π L_{i} (\sin θ^{2} + 0.5 \cos θ^{2})}{4 N \cos θ} \\ \times {[d + (2 i - 1) (\frac{L_{i}}{N} \tan θ)]}^{3} \end{cases}

(7)

in the formula,

p_{i}

is the average contact pressure of the

i^{t h}

unit;

A_{i}

is the contact area of the

i^{t h}

unit;

L_{i}

is the fitting distance of the

i^{t h}

conical surface distance; and

α, β

are the contact surface characteristic parameters.

2.1.2. Dynamic stiffness model of spindle-holder joint

When the spindle is stationary, the spindle-holder joint surface produces an interference fit under the tension of the broach rod, forming an initial interference amount $δ_{0}$ . According to the elastic deformation theory, the expression is:

δ_{0} = 2 (u_{s} - u_{h}) = \frac{4 a b^{2} P_{0}}{E (b^{2} - a^{2})}

(8)

{\begin{cases} u_{s} = \frac{a^{2} P}{E (b^{2} - a^{2})} [(1 - υ) a + \frac{(1 + υ^{2}) b^{2}}{a}] \\ u_{h} = - \frac{P}{E} [(1 - υ) a] \end{cases}

(9)

In the formula:

u_{s}

a n d u_{h}

are the deformation of the spindle and the tool holder at any contact position, respectively;

P_{0}

is the static pressure at any position on the conical surface of the spindle-tool holder joint;

a

is the radius of any contact position of the conical part of the tool holder; and

b

is the outer diameter of the spindle taper hole part.

E

is the elastic modulus and

υ

is Poisson’s ratio.

When the spindle rotates at high speed, due to the hollow conical structure of the spindle joint and the solid conical structure of the tool holder, the deformation difference between the spindle and the tool holder at any contact position will cause the axial distribution of the interference at the joint to change. This results in a centrifugal expansion effect as shown in Figure 5.

Figure 5.

Simplified diagram of the centrifugal expansion of the spindle-holder joint.

According to the elastic deformation theory (Xu and Liu, 1995), the deformation expression of the spindle-holder joint at any contact point on the tapered surface shown in Figure 5 is:

{\begin{cases} u_{s^{'}} = \frac{(3 - 2 υ) (1 + υ)}{4 E} ρ ω^{2} a [b^{2} + \frac{1 - 2 υ}{3 - 2 υ} a^{2}] \\ u_{h^{'}} = \frac{(1 - 2 υ) (1 + υ)}{4 E} ρ ω^{2} a^{3} \end{cases}

(10)

u_{s^{'}}, u_{h^{'}}

are, respectively, the deformation of the spindle and tool holder at any contact position under the centrifugal expansion effect.

Furthermore, according to formula (10), the expression for the radial gap $δ_{1}$ is:

δ_{1} = u_{s^{'}} - u_{h^{'}} = \frac{(3 - 2 υ) (1 + υ)}{4 E} ρ ω^{2} a b^{2}

(11)

The dynamic contact pressure $P$ at any position of the spindle-holder joint can be expressed as the result of the interaction between the static interference amount $δ_{0}$ and the radial gap $δ_{1}$ caused by the centrifugal expansion effect. According to equations (8) and (11), the expression of the contact pressure $P$ at any position of the spindle-tool holder joint contacting the tapered surface under dynamic conditions is:

P = \frac{E (δ_{0} - δ_{1}) (b^{2} - a^{2})}{4 a b^{2}}

(12)

From this, the expression of the contact pressure $P$ at any position of the joint at any rotation speed can be established as:

P = P_{0} - \frac{ρ (b^{2} - a^{2}) (3 - 2 υ) (1 + υ)}{16} ω^{2}

(13)

From equation (13), the variable stiffness model of the spindle-holder joint can be obtained by integrating the pressure of each shaft segment unit and bringing it into equation (7).

2.2. Construction of dynamic model of holder-tool joint

The holder-tool joint is structurally composed of a tool holder, a collet, a lock nut, and a tool. Torque is applied through the lock nut, and the tool is assembled to the tool holder through an elastic jacket. In order to describe the contact relationship between the tool holder and the tool, a distributed spring-damping model is used to replace the jacket to simulate the contact characteristics between the tool holder and the tool, as shown in Figure 6.

Figure 6.

Distributed spring-damping model of tool holder-tool joint.

Assuming a uniform distribution of pressure at the contact interface among the tool, tool holder, jacket, and preload nut joint, where the tool diameter is denoted as $d_{t o o l}$ and the contact length between the tool holder and the tool is represented as $L_{2}$ , the concentrated force at each interface can be determined as illustrated in Figure 7.

Figure 7.

Stress analysis diagram of holder-tool joint.

The expressions of each concentration force in Figure 7 are as follows:

{\begin{cases} M = F_{e} k_{0} d \\ F_{n_{2}} \sin θ_{3} + F_{μ_{2}} \cos θ_{3} = F_{e} \\ F_{n 2} \sin θ_{3} + F_{μ 2} \cos θ_{3} = F_{n 1} \sin θ_{2} + F_{μ 1} \cos θ_{2} + F_{μ 3} \\ F_{n 2} \cos θ_{3} - F_{μ 2} \sin θ_{3} + F_{n 1} \cos θ_{2} - F_{μ 1} \sin θ_{2} = F_{n 3} \end{cases}

(14)

In the formula,

M

is the locking torque,

F_{e}

is the equivalent preload force, and

k_{0}

is the preload coefficient;

F_{n_{3}}

is the holder-tool joint surface positive pressure,

F_{μ 3}

is the holder-tool joint surface friction force, its magnitude is

μ_{3} F_{n 3}

, and

μ_{2}

is the jacket surface friction factor;

F_{n 1}

is the jacket-holder conical surface positive pressure,

F_{μ 1}

is the jacket-holder conical surface friction force, its value is

{μ_{1} F}_{n 1}

, and

μ_{1}

is the friction factor of the collet-tool holder taper surface;

F_{n 2}

is the positive pressure of the collet-nut contact surface,

F_{μ 2}

is the friction force of the collet-nut contact surface, and its value is

{μ_{2} F}_{n 2}

, and

μ_{2}

is the friction factor between the collet and the nut;

d

is the major diameter of the round nut thread,

θ_{2}, θ_{3}

are the collet-tool holder taper angle, and the collet-nut taper angle.

By solving equation (14), the positive pressure $F_{n 3}$ of the collet-tool interface can be obtained, which is equivalent to the positive pressure on the tool holder-tool interface, and the jacket-tool contact area is equivalent to the holder-tool contact area. Therefore, based on the previous assumption that the pressure on the contact interface is uniformly distributed, the average contact pressure on the holder-tool interface can be obtained as $P_{n 3}$ :

p_{n 3} = \frac{F_{n 3}}{A_{3}}

(15)

In the formula,

A_{3}

is the contact area of the tool holder-tool joint surface. The contact stiffness of each segmented unit can be obtained by discretizing the holder-tool joint and bringing it into equation (7).

2.3. Electromagnetic stiffness model of active magnetic bearing

This paper uses the equivalent magnetic circuit analysis method to model the electromagnetic force of the active magnetic bearing, and divides the radial magnetic bearing into four parts: stator yoke, magnetic pole, air gap, and rotor yoke, as shown in Figure 8:

Figure 8.

Structural diagram of active magnetic bearing.

Based on the static suspension position of the rotor, during the high-speed rotation of the spindle, due to uneven mass distribution and external unbalanced excitation, the rotor’s spatial position shifts, as shown in Figure 9.

Figure 9.

Schematic diagram of rotor offset of magnetic bearing.

The offset of the rotor from the balance reference position is expressed by the eccentric angle $φ$ and the displacement $e$ . When the rotor is in the balanced position, the normal air gap length between the rotor and the stator is $g_{0}$ , and $g (θ)$ is the air gap length at the circumferential position of the air gap $θ$ . According to geometric relationships, when the rotor is eccentric, the expression of $g (θ)$ can be obtained:

g (θ) = g_{0} - e \cos (θ - φ)

(16)

When the coupled and saturated magnetic flux density is not considered, combined with reference (Jiang et al., 2020), the magnetic flux density

F_{i}

under the initial air gap of each node can be obtained based on the magnetomotive force

g_{(θ_{i})}

at each magnetic pole and the air gap length

B_{g, θ}

{\begin{cases} \begin{array}{l} B_{g, θ} = \frac{μ_{0} F_{i}}{g (θ_{i})} λ_{i}, i \in (1, 8) \\ F_{i} = N I_{0} + N I_{i} \end{array} \\ θ_{i} = \frac{π}{8} + \frac{π}{4} (i - 1) \end{cases}

(17)

In the formula,

μ_{0}

is the vacuum permeability,

F_{i}

is the magnetomotive force at the circumferential node position,

I_{i}

is the control current,

I_{0}

is the bias current,

θ_{i}

is the air gap deflection angle at the center of the magnetic pole, and

λ_{i}

is the leakage compensation coefficient.

In order to obtain the total electromagnetic carrying capacity $F (I, e, φ)$ of the active magnetic bearing, the air gap is divided into n position points at equal intervals along the circumferential direction of the magnetic bearing. The air gap magnetic flux density at each position is calculated. The force vector at each position point is superimposed to obtain the total electromagnetic force. At the same time, the total electromagnetic force is decomposed into two components in the x and y directions:

{\begin{cases} F (I, e, φ) = \sqrt{F_{x}^{2} + F_{y}^{2}} \\ F_{x} = \sum_{θ = 0}^{2 π} \frac{{B_{g, θ}}^{2} A_{θ}}{2 μ_{0}} \cos θ \\ F_{y} = \sum_{θ = 0}^{2 π} \frac{{B_{g, θ}}^{2} A_{θ}}{2 μ_{0}} \sin θ \end{cases}

(18)

In the formula,

A_{θ}

is the equivalent magnetic pole area at the offset position. According to the electromagnetic force, the active magnetic bearing current stiffness and displacement stiffness are calculated:

{\begin{cases} k_{I} = \frac{\partial F (I, d, φ) {​ |}_{d = 0, φ = 0}}{\partial I} \\ k_{S} = \frac{\partial F (I_{i}, d, φ) {​ |}_{I = 0}}{\partial d} \end{cases}

(19)

For different currents and rotor positions, the current stiffness $k_{I}$ and displacement stiffness $k_{S}$ of the magnetic bearing at different positions in the axial and radial directions can be obtained:

{\begin{cases} K_{s} \approx \frac{μ_{0} N^{2} A_{α} I_{0}^{2}}{g_{0}^{3}} \\ K_{I} \approx \frac{μ_{0} N^{2} A_{α} I_{0}}{g_{0}^{2}} \end{cases}

(20)

2.4. Milling force model

The magnetic levitation milling electric spindle will be affected by the milling force load during the cutting process. The milling force formula in the literature (Tang, 2014) is quoted here to decompose the milling force on the micro-element cutting edge into tangential milling force ${d f}_{t}$ and radial milling force. ${d f}_{r}$ The expression is:

{\begin{cases} d f_{t} (i, j, k) = K_{t c} h (i, j, k) d z + K_{t p} d s \\ d f_{r} (i, j, k) = K_{r c} h (i, j, k) d z + K_{r p} d s \end{cases}

(21)

In the formula,

K_{t c} 、 K_{r c}

represent the tangential and radial shear force coefficients, respectively,

K_{t p} 、 K_{r p}

represent the tangential and radial plowing force coefficients, respectively,

h (i, j, k)

is the instantaneous cutting thickness of a single micro-element cutting edge, and

d s = d z / \cos β

is the length of the micro-element cutting edge. Decompose the cutting force on the micro-element cutting edge with the spatial position angle

φ (i, j, k)

into the local coordinate system

O x y z

of the tool, and the following expression can be obtained:

[\begin{array}{l} d f_{x} (i, j, k) \\ d f_{y} (i, j, k) \end{array}] = T [\begin{array}{l} d f_{t} \\ d f_{r} \end{array}]

(22)

T = [\begin{array}{l} - \sin φ (i, j, k) & \cos φ (i, j, k) \\ - \cos φ (i, j, k) & - \sin φ (i, j, k) \end{array}]

(23)

The total milling force component is obtained by summing the components of the milling force in the x and y directions of all micro-elements participating in the cutting:

{\begin{cases} f_{x} (j) = \sum_{i = 1}^{N_{i}} \sum_{k = 1}^{N_{k}} d f_{x} (i, j, k) \\ f_{y} (j) = \sum_{i = 1}^{N_{i}} \sum_{k = 1}^{N_{k}} d f_{y} (i, j, k) \end{cases}

(24)

3. Dynamic characteristic analysis

To investigate the influence of the spindle-holder and holder-tool joints on the dynamic characteristics of the magnetic levitation milling spindle, simulations are performed using the dynamic model described previously. The operating speed is assumed to be 10,000 rpm, with feed per tooth 0.06 mm, radial depth of cut 2 mm, and axial depth of cut 2 mm. The equations of motion are solved using the ode45 solver in MATLAB. System structural parameters are listed in Tables 1 and 2.

Table 1.

Structural parameters and contact stiffness of spindle-holder joint.

$D_{1} / m m$	$d_{1} / m m$	$L_{1} / m m$	$θ_{1} / (°)$
44.45	28.4	55	8.297
$f_{1} / K N$	$K_{n 1} (N / m)$	$K_{τ 1} (N / m)$	$K_{θ 1} (N \cdot m / r a d)$
13	$1.69 \times 10^{10}$	$8.45 \times 10^{9}$	$3.02 \times 10^{6}$

Table 2.

Structural parameters and contact stiffness of the holder-tool joint.

$d_{t o o l} / m m$	$L_{2} / m m$	$M (N \cdot m)$	$K_{n 2} (N / m)$
20	40	60	$1.16 \times 10^{10}$
$K_{τ 2} (N / m)$	$K_{θ 2} (N \cdot m / r a d)$
$5.78 \times 10^{9}$	$4.5 \times 10^{5}$

The data in Tables 1 and 2 are based on the structural parameters of the BT40 holder, which reference the international standard ISO 7388-2:2007(E). In the table: $D_{1}$ is the diameter of the large end of the tool holder; $d_{1}$ is the diameter of the small end of the tool holder; $K_{n 1}$ , $K_{τ 1}$ , and $K_{θ 1}$ are the total normal, tangential, and angular stiffness of the spindle-holder joint, respectively; $d_{t o o l}$ is the diameter of the tool; $L_{2}$ is the contact length between the tool and the tool holder; and $K_{n 2}$ , $K_{τ 2}$ , and $K_{θ 2}$ are the total normal, tangential, and angular stiffness of the holder-tool joint, respectively.

Based on formulas (7) and (13), three sets of initial tool pulling forces are designated as the initial tool pulling forces. Simulation analysis reveals that with an increase in spindle speed, the contact stiffness and angular stiffness at the spindle-tool holder joint exhibit a notable decreasing trend, as illustrated in Figure 10.

Figure 10.

Stiffness variation diagram of the spindle-holder joint.

3.1. Modal analysis

Utilizing the dynamic model of the magnetic levitation milling electric spindle, which incorporates the spindle-holder-tool junction established previously, the first three natural modes of the system and their associated mode shapes were determined using MATLAB software, as illustrated in Figures 11 –13.

Figure 11.

First-order modal shape diagram.

Figure 12.

Modal shape diagram (8000 r/min).

Figure 13.

Modal shape diagram (10,000 r/min).

For the first-order mode, the system exhibits dominant rigid-body behavior, with a relatively uniform mode shape distribution along the structure. In this case, the deformation energy is minimal and not concentrated in the spindle-holder-tool joint region. As a result, the influence of joint stiffness variation on this mode is negligible. This is confirmed by the simulation results, where the first-order natural frequency shows only a slight change from $106.2047 Hz$ to $106.184 Hz$ (a reduction of $0.019 %$ ) as the rotational speed increases from $8000 rpm$ to $10, 000 rpm$ .

In contrast, the second-order mode is characterized by significant bending deformation, indicating the flexible dynamic behavior of the rotor system. Due to the abrupt change in cross-sectional geometry and moment of inertia between the spindle and the tool holder-tool assembly, a considerable portion of the modal strain energy is concentrated in the joint region. Consequently, the joint becomes a critical compliance location that strongly influences the modal characteristics. Under the effect of centrifugal expansion, the degradation of joint stiffness leads to a pronounced reduction in the second-order natural frequency, decreasing from $287.925 Hz$ to $242.1125 Hz$ (a reduction of $15.91 %$ ).

For the third-order mode, the deformation pattern becomes more complex, with multiple nodes appearing near the tool tip and the joint interface. In this higher-frequency range, the effects of shear deformation and rotary inertia, as described by Timoshenko beam theory, become more significant. However, the modal strain energy is distributed over a wider region of the structure, reducing the dominance of the joint region. As a result, the influence of local stiffness degradation is less pronounced, and the third-order natural frequency decreases from $1451.1382 Hz$ to $1392.6003 Hz$ (a reduction of $4.03 %$ ).

These modal variations are closely related to the degradation of joint stiffness induced by centrifugal effects. Specifically, as the rotational speed increases from $8000 rpm$ to $10, 000 rpm$ , the normal contact stiffness of the spindle-holder-tool joint decreases from ${1.14309 \times 10}^{10} N / m$ to ${8.78228 \times 10}^{9} N / m$ (a reduction of $23.17 %$ ), while the angular stiffness decreases from ${2.09644 \times 10}^{6} N \cdot m / r a d$ to ${1.6521 \times 10}^{6} N \cdot m / r a d$ (a reduction of $21.19 %$ ). According to structural dynamics principles, the effective stiffness of the system is largely governed by its most compliant region. Therefore, the reduction in joint stiffness directly leads to the observed decrease in natural frequencies, particularly for higher-order modes with strong strain energy localization.

3.2. Dynamic response analysis

To elucidate the impact of variations in contact stiffness at the spindle-holder-tool junction on the dynamic behavior of the magnetic levitation milling electric spindle, the milling force was employed as the external excitation to determine the time-domain vibration response of the tool tip node across different rotational speed conditions. The vibration responses at the tool tip point are presented in Figures 14 and 15.

Figure 14.

Vibration response diagram of the tool tip (8000 r/min).

Figure 15.

Vibration response diagram of the tool tip (10,000 r/min).

When the spindle speed is $8000 rpm$ , the displacement time-domain responses of the tool tip node in the x and y directions are shown in Figure 14(a) and (b). After the electric spindle starts and reaches a stable rotational speed, the milling cutter enters the cutting state. Under the action of periodic milling force, the vibration amplitude at the tool tip node increases rapidly, and the system’s equilibrium position deviates to a certain extent from the initial center axis. Subsequently, the system enters a stable cutting state.

When the spindle speed is increased to $10, 000 rpm$ , the vibration displacement of the tool tip node is shown in Figure 15(a) and (b). Compared with the $8000 rpm$ working condition, the displacement amplitudes in both the x and y directions have significantly increased. This is because as the spindle speed increases, the centrifugal softening effect causes the contact pressure at the joint to decrease. According to equation (13), the decrease in contact pressure will directly lead to the synchronous attenuation of the normal stiffness, tangential stiffness, and angular stiffness of the joint surface. The reduction of the equivalent stiffness at the junction weakens the overall bending resistance and vibration resistance of the spindle-holder-tool system, causing the tool tip to generate a greater dynamic displacement response under the milling force excitation.

4. Conclusion

In this paper, a dynamic model of a magnetic levitation milling spindle system considering the variable stiffness characteristics of the spindle-holder-tool joint was established. The effects of centrifugal expansion and rotational speed on joint stiffness and system dynamic response were systematically investigated. The results show that, as the rotational speed increases from $8000 rpm$ to $10, 000 rpm$ , the normal contact stiffness of the joint decreases by $23.17 %$ , while the angular stiffness decreases by $21.19 %$ due to centrifugal effects. This stiffness degradation leads to a reduction in system natural frequencies, where the second-order natural frequency decreases by $15.91 %$ and the third-order natural frequency decreases by $4.03 %$ , while the first-order frequency remains nearly unchanged. Modal analysis reveals that the sensitivity of the system to joint stiffness variation is strongly mode-dependent. The first-order mode, dominated by rigid-body motion, shows minimal sensitivity to stiffness changes. In contrast, higher-order bending modes exhibit significant sensitivity due to the concentration of modal strain energy in the joint region. As a result, stiffness degradation at the joint has a pronounced influence on system dynamics.

Furthermore, time-domain analysis indicates that the reduction in joint stiffness leads to a noticeable increase in tool tip vibration amplitude at higher rotational speeds. This demonstrates that joint stiffness plays a critical role in determining the dynamic stability of the spindle system.

The findings of this study provide important theoretical support for the design and optimization of high-speed magnetic levitation milling spindles. In particular, controlling joint stiffness and mitigating its degradation under centrifugal effects are essential for improving vibration performance and ensuring machining stability.

Footnotes

ORCID iDs

Rui Yang

Huimin Kang

Quan Zhou

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China grant number 51875198, the Education Department of Hunan Province grant number 23A0360, and the Xiangtan Science and Technology Bureau grant number GX-YB20221004.

Declaration of conflicting interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

(25)

{\begin{cases} \frac{\partial σ_{r}}{\partial r} + \frac{1}{r} \frac{\partial τ_{r φ}}{\partial_{r}} + \frac{σ_{r} - σ_{φ}}{r} + F_{r} = 0 \\ \frac{1}{r} \frac{\partial σ_{φ}}{\partial φ} + \frac{\partial τ_{r φ}}{\partial_{r}} + \frac{2 τ_{r φ}}{r} + F_{φ} = 0 \end{cases}

in the formula:

σ_{r}

is the radial normal stress;

σ_{φ}

is the circumferential normal stress;

τ_{r φ}

is the shear stress, and its value is 0;

F_{r}

is the radial force component; and

F_{φ}

is the component of the circumferential force, and its value is 0.

Simplifying equation (25), we can obtain: (26)

\frac{\partial σ_{r}}{\partial r} + \frac{σ_{r} - σ_{φ}}{r} + ρ r ω^{2} = 0

The geometric equation is: (27)

ε_{r} = \frac{\partial u_{r}}{\partial r}, ε_{φ} = \frac{u_{r}}{r} + \frac{1}{r} \frac{\partial u_{r}}{\partial r}

The physical equation is: (28)

{\begin{cases} ε_{r} = \frac{1 - υ^{2}}{E} (σ_{r} - \frac{υ}{1 - υ} σ_{φ}) \\ ε_{φ} = \frac{1 - υ^{2}}{E} (σ_{φ} - \frac{υ}{1 - υ} σ_{r}) \end{cases}

According to Hooke’s Law, it can be obtained that: (29)

\frac{d σ_{r}}{d r} - υ \frac{d σ_{r}}{d r} = \frac{1 + υ}{r} (σ_{r} - σ_{φ})

in the formula:

ε_{r}

is the radial shear strain;

ε_{φ}

is the circumferential shear strain;

u_{r}

is the radial displacement;

E

is the elastic modulus; and

υ

is poisson’s ratio.

The boundary conditions can be obtained from Figure 16. (30)

{\begin{cases} h o l d e r : σ_{r | r = d} = 0, σ_{r | r = c} = 0 \\ s p i n d l e : σ_{r | r = b} = 0, σ_{r | r = d} = 0 \end{cases}

In the formula: $c$ is the inner diameter of the cone of the holder; $d$ is the outer diameter of the holder cone (the inner diameter of the spindle taper hole); and $b$ is the outer diameter of the main taper hole. Figure 16.

Schematic diagram of the taper hole fit between the spindle and the holder.

References

Cui

(2020) Chatter stability prediction method of the spindle-tool holder-tool system with interface contact characteristics. Mathematical Problems in Engineering 2020: 15. https://doi.org/10.1155/2020/7121328

Zhang

, et al. (2015) Ultrasonic measurement of contact stiffness and pressure distribution on spindle–holder taper interfaces. International Journal of Machine Tools and Manufacture 97: 18–28. https://doi.org/10.1016/j.ijmachtools.2015.06.007

Ertürk

Ozgüven

Budak

(2006) Analytical modeling of spindle–tool dynamics on machine tools using Timoshenko beam model and receptance coupling for the prediction of tool point FRF. International Journal of Machine Tools and Manufacture 46(15): 1901–1912. https://doi.org/10.1016/j.ijmachtools.2006.01.032

Ertürk

Ozgüven

Budak

(2007a) Effect analysis of bearing and interface dynamics on tool point FRF for chatter stability in machine tools by using a new analytical model for spindle–tool assemblies. International Journal of Machine Tools and Manufacture 47(1): 23–32. https://doi.org/10.1016/j.ijmachtools.2006.03.001

Ertürk

Budak

Ozgüven

(2007b) Selection of design and operational parameters in spindle-holder-tool assemblies for maximum chatter stability by using a new analytical model. International Journal of Machine Tools and Manufacture 47(9): 1401–1409. https://doi.org/10.1016/j.ijmachtools.2006.08.016

Gagnol

Ray

(2011) Modal identification of spindle-tool unit in high-speed machining. Mechanical Systems and Signal Processing 25(7): 2388–2398. https://doi.org/10.1016/j.ymssp.2011.02.019

Iglesias

Taner

Zsahin

, et al. (2022) Alternative experimental methods for machine tool dynamics identification: a review. Mechanical Systems and Signal Processing 170: 108837. https://doi.org/10.1016/j.ymssp.2022.108837

Jiang

Wang

(2020) Analytical calculation of active magnetic bearing based on distributed magnetic circuit method. IEEE Transactions on Energy Conversion 36(3): 1. https://doi.org/10.1109/tec.2020.3040975

Kiran

(2022) A receptance coupling procedure considering frequency-dependent behavior of holder-tool contact dynamics. Journal of Manufacturing Processes 80: 624–641. https://doi.org/10.1016/j.jmapro.2022.06.007

10.

Kilic

Altintas

(2020) General cutting dynamics model for five-axis ball-end milling operations. Journal of Manufacturing Science and Engineering 142(12): 121003. https://doi.org/10.1115/1.4047625

11.

Liu

Wang

, et al. (2019) Contact stiffness of spindle-tool holder based on fractal theory and multi-scale contact mechanics model. Mechanical Systems and Signal Processing 119: 363–379. https://doi.org/10.1016/j.ymssp.2018.09.037

12.

Luo

Chen

, et al. (2022) FRF prediction modeling and joint parameter identification of motorized spindle-tool handle. Journal of the Brazilian Society of Mechanical Sciences and Engineering 44(7): 303. https://doi.org/10.1007/s40430-022-03552-5

13.

Miao

, et al. (2023) A fully analytical nonlinear dynamic model of spindle-holder-tool system considering contact characteristics of joint interfaces. Mechanical system and signal processing 202: 110693. https://doi.org/10.1016/j.ymssp.2023.110693

14.

Namazi

Altintas

Abe

, et al. (2007) Modeling and identification of tool holder-spindle interface dynamics. International Journal of Machine Tools and Manufacture 47(9): 1333–1341. https://doi.org/10.1016/j.ijmachtools.2006.08.003

15.

Özşahin

Ertürk

Özgüven

, et al. (2009) A closed form approach for identification of dynamical contact parameters in spindle-holder-tool assemblies. International Journal of Machine Tools and Manufacture 49(1): 25–35. https://doi.org/10.1016/j.ijmachtools.2008.08.007

16.

Sakamoto

Matsuda

Shimizu

(2012) Effect of clamped holders on dynamic characteristics of spindle system of machining center. International Journal of Automation Technology 6(2): 168–174. https://doi.org/10.20965/ijat.2012.p0168

17.

Salahshoor

Ahmadian

(2009) Continuous model for analytical prediction of chatter in milling. International Journal of Machine Tools and Manufacture 49(14): 1136–1143. https://doi.org/10.1016/j.ijmachtools.2009.07.016

18.

Shamine

Shin

(2000) Analysis of no. 50 taper joint stiffness Under axial and radial loading. Journal of Manufacturing Processes 2(3): 167–173. https://doi.org/10.1016/s1526-6125(00)70118-6

19.

Tang

(2014) Three-Dimensional Stability and Parameter Optimization in Milling of Thin-Walled Parts. Beijing: National Defense Industry Press, 31–35.

20.

Liu

(1995) Applied Elastic and Plastic Mechanics. Beijing: Tsinghua University Press, 237–244.

21.

Zhang

Feng

, et al. (2014) Characteristics of stiffness and contact stress distribution of a spindle–holder taper joint under clamping and centrifugal forces. International Journal of Machine Tools and Manufacture 82-83: 21–28. https://doi.org/10.1016/j.ijmachtools.2014.03.006

22.

Zhao

Zhang

(2013) Analytical modeling method of dynamics for the spindle rotor system. Journal of Mechanical Engineering 49(6): 44. https://doi.org/10.3901/jme.2013.06.044

23.

Zhao

Jiang

, et al. (2023) Static model and reliability analysis of double nut ball screw considering screw bending. Tribology International 187: 108733. https://doi.org/10.1016/j.triboint.2023.108733

24.

Zhao

Gong

Zhao

, et al. (2024) Investigations on contact characteristics of ball screw considering flexible deformation of screw and nut. Nonlinear Dynamics 112(20): 17745–17779. https://doi.org/10.1007/s11071-024-09949-2

25.

Zhao

Gong

, et al. (2025a) Nonlinear vibration characteristics and reliability analysis of dynamic model of linear motion platform supported by double rolling linear guide rails. Mechanical Systems and Signal Processing 229: 112507. https://doi.org/10.1016/j.ymssp.2025.112507

26.

Zhao

Gong

Yao

, et al. (2025b) Reliability-based design optimization of key components in the ball screw feed system. Mechanism and Machine Theory 214: 106151. https://doi.org/10.1016/j.mechmachtheory.2025.106151