Sage Journals: Discover world-class research

Abstract

Bi-rotary milling head is one of the core components of five-axis machining center, and its dynamic characteristics directly affect the machining stability and accuracy. During the sculptured surface machining, the bi-rotary milling head exhibits varying dynamics in various machining postures. To facilitate rapid evaluation of the dynamic behavior of the bi-rotary milling head within the whole workspace, this article presents a method for parametrically establishing dynamic equation at different postures. The rotating and swing shafts are treated as rigid bodies. The varying stiffness of the flexible joints (such as bearings and hirth coupling) affected by gravity and cutting force at different swing angles is analyzed and then a multi-rigid-body dynamic model of the bi-rotary milling head considering the pose-varying joint stiffness is established. The Lagrangian method is employed to deduce the parametric dynamic equation with posture parameters. The static stiffness, natural frequencies and frequency response functions at different postures are simulated using the parametric equation and verified by the impact testing experiments. The theoretical and experimental results show that the dynamics of the bi-rotary milling head vary with the machining postures, and the proposed method can be used for efficient and accurate evaluation of the pose-dependent dynamics at the design stage.

Keywords

Bi-rotary milling head joint stiffness dynamic model parametric modeling method pose-dependent dynamics

Introduction

Bi-rotary milling head with both rotation and swing degrees of freedom (DOFs) is one of the typical configurations of five-axis machine tool and is always used to achieve the machining of large parts with complex surfaces, especially in aerospace industry. In the machining process, many factors affect the machining accuracy such as kinematics errors, thermal errors, dynamic characteristic and tool wear. In the past decades, much research has been focused on the geometric error compensation of five-axis machine tools^1–4 and the thermal characteristics of milling head,⁵ aiming to improve the machining accuracy.

During the five-axis machining, however, the milling head would change the posture continuously. The different positions and postures may lead to the different static and dynamic characteristics for machine tools. Dhupia et al.⁶ performed modal tests on the arch-type reconfigurable machine tool and found the frequency response functions (FRFs) are different from one reconfiguration position (θ = 0°) to another (θ = 45°). Bonnemains et al.⁷ and Wu et al.⁸ studied the hybrid machine tool and concluded that the static stiffness behaves differently owing to the position variation. Law et al.⁹ also found that the natural frequency and dynamic stiffness of the three-axis milling machine behave position dependent. The position-varying dynamics have a great influence on machining stability^10,11 and control performance,¹² which limit the achievable productivity and machining trajectories in the whole working range of the machine.

In order to evaluate the varying dynamic behavior at the design stage, some modeling methods have been presented. The most common method is the finite element model (FEM) which is quite convenient to consider the structural flexibility.^13–15 But it needs to establish various dynamic models for different configurations and has to be remeshed and recomputed, which is very time-consuming. To improve the modeling efficiency, the co-simulation techniques (CSTs) were proposed to reduce the transition time of modeling from one state to another. One is to integrate FEM and multi-body simulation,^16,17 and the other is the integration of FEM and Kriging optimal interpolation.¹⁸ Substructure synthesis technique (SST) was also used to investigate the varying dynamics of machine tools.^19–21 In this method, those components are assembled to different configurations of machine tools by small DOFs. Consequently, its computational time and modeling efficiency are improved.

The component flexibility of machine tools can be treated conveniently in the three methods mentioned above, while the contact characteristics between joints or substructures were often ignored. One of the reasons for the varying dynamics of mechanical systems attributes to the variation of the system mass distribution;²² however, the other is stiffness variation of joints. Cao et al.²³ found that the speed-varying stiffness of bearing joints could lead to speed-varying dynamics of spindles. Yigit and Ulsoy²⁴ and Dhupia et al.²⁵ demonstrated that the weakly nonlinear stiffness of machine joints can cause significant variations of natural frequency and dynamic stiffness. The lumped parameter model (LPM), therefore, is proposed to describe the varying dynamics based on varying stiffness of joints. Zhang et al.²⁶ studied the effect of varying stiffness of kinematic joints with different feed rates on dynamics of the high-speed ball screw system. Wei et al.²⁷ discussed that the varying joint stiffness of the worm and worm wheel at different tilting angles caused the different dynamics of the tilting table. Wang et al.²⁸ analyzed the effect of axis coupling force on the stiffness changes of kinematic joints and the varying dynamic behavior of a three-axis gantry milling machine. Comparing to other approaches, LPM is more suitable for the parameterized analysis and reanalysis for the varying dynamics only due to the varying joint stiffness.

FRF at tool center point (TCP) reflects the dynamics of the entire machine tool, and its dominant mode determines the limit cutting depth during cutting process.²⁹ The machine frame properties can limit FRFs at TCP.^30,31 The milling head directly supports the spindle, and its dynamics inevitably affect the dynamic characteristics of TCP. Hung et al.³² confirmed the influence of the spindle orientation on tool point dynamics and machining stability. Law et al.³³ proposed a dynamic substructuring procedure to model the orientation-dependent dynamic behavior of machine tools with gimbal heads. Yan et al.³⁴ and Shi et al.³⁵ also found the variable dynamics of milling head at different postures by FEM. In the modeling process mentioned above,^32–35 the joint stiffness was ignored. However, the bi-rotary milling head system contains a large number of joints, such as rotating and swing shaft bearings, hirth coupling and locking joints. Those joints’ stiffness is often less than that of components and some of them vary with different postures. Therefore, to accurately evaluate the varying dynamics of the bi-rotary milling head system at different postures, the joint stiffness has to be considered. This research aims to propose a rapid modeling method to accurately evaluate the dynamic behavior of bi-rotary milling head in the whole working range. A multi-rigid-body dynamic model considering the joint stiffness is established, and the varying stiffness of flexible joints at different postures is analyzed. Finally, a parametric dynamic equation with posture parameters is deduced. Compared with other modeling methods, the pose-varying joint stiffness is considered in the multi-body dynamic model, and the explicit expression of dynamic equation is derived, which can easily and rapidly evaluate the dynamic behavior within the whole workspace.

This article is organized as follows: a dynamic equation with rotation and swing angles is presented in section “Dynamic equation of milling head with rotation and swing angles.” The nonlinear joint stiffness in different machining postures is calculated in section “Calculation of flexible joint stiffness.” The pose-dependent dynamics are analyzed in section “Pose-dependent dynamics.” In section “Experimental verification,” a direct-drive bi-rotary milling head is used to demonstrate the accuracy of the proposed method, followed by conclusions in section “Conclusion.”

Dynamic equation of milling head with rotation and swing angles

Equivalent dynamic model

There are mainly two kinds of bi-rotary milling head, the fork type (Figure 1(a)) and the bias type (Figure 1(b)). The only difference is that the spindle axis overlaps or offsets with the rotation axis.

Figure 1.

Bi-rotary milling head: (a) the fork type and (b) the bias type.

The fork-type direct-drive milling head considered here generally includes the rotating shaft C-axis, the swing shaft A-axis, A/C-axis drive motors, A/C-axis bearings and A/C-axis locking devices, as shown in Figure 2(a).

Figure 2.

Fork-type bi-rotary milling head: (a) geometric model, (b) dynamic model and (c) calculating model.

Generally, structural stiffness in machine tools is much stiffer than that of the joints and the dominant modes in low frequencies entirely dependent on those flexible joints.^28,36 Consequently, for the bi-rotary milling head, the shafts and housings are assumed to be a rigid body due to their larger stiffness at low frequencies, while the bearings are always treated as the flexible coupling. Positioning and locking of A/C-axis in different postures are usually accomplished by mechanical or servo control method. In this case, the stiffness of the locking joints is obviously smaller than that of the shafts, and therefore, they are also considered as the flexible coupling. Although the motorized spindle is a flexible part, it is assumed as a rigid rod and rigidly coupled to the spindle housing for investigating the dynamic characteristics of the bi-rotary milling head itself within the workspace.

Ignoring the motor inductance and all kinds of friction while A/C-axis rotating, we can obtain a multi-body dynamic model for the bi-rotary milling head system, including four rigid rods, eight flexible couplings and one rigid coupling, as shown in Figure 2(b). The four rigid rods are the C-axis (A), A-axis fork housing (B), A-axis (C) and motorized spindle (D), respectively. The eight flexible couplings are listed as the rear bearing (no. 1), front bearing (no. 2), hirth coupling (no. 3), left bearing (no. 4), right bearing (no. 5), C-axis locking joint (no. 6), A-axis left locking joint (no. 7) and A-axis right locking joint (no. 8). One rigid coupling is the connection of the motorized spindle and spindle housing.

The rods A and B are characterized as inertia 1 ( m ₁, J ₁) and inertia 2 ( m ₂, J ₂), respectively (see Figure 2(c)). The rigidly coupled rods C and D are equivalent to inertia 3 ( m ₃, J ₃). Under the assumption of small displacement, the flexible coupling is equivalent to spring and damping elements. The front bearing, rear bearing and C-axis locking joint are further treated as the equivalent joint A (EJ-A) due to their action on the C-axis. According to the series–parallel relationship, the equivalent stiffness of EJ-A is expressed as

{\begin{matrix} {k_{q 1}}_{x} = {k_{1}}_{x} + {k_{2}}_{x} + {k_{c 1}}_{x} \\ {k_{q 1}}_{y} = {k_{1}}_{y} + {k_{2}}_{y} + {k_{c 1}}_{y} \\ {k_{q 1}}_{z} = {k_{1}}_{z} + {k_{2}}_{z} + {k_{c 1}}_{z} \\ {k_{q 1}}_{θ x} = {k_{1}}_{θ x} + {k_{12}}_{θ x} + {k_{c 1}}_{θ x} + k_{1 y} l_{1 s 1}^{2} + k_{2 y} l_{1 s 2}^{2} + k_{c 1 y} l_{1 c}^{2} \\ {k_{q 1}}_{θ y} = {k_{1}}_{θ V} + {k_{2}}_{θ x} + {k_{c 1}}_{θ y} + k_{1 x} l_{1 s 1}^{2} + k_{2 x} l_{1 s 2}^{2} + k_{c 1 x} l_{1 c}^{2} \\ {k_{q 1}}_{θ z} = {k_{c 1}}_{θ z} \end{matrix}

(1)

where k_1x, k_1y, k_1z, k_1θx and k_1θy are the stiffness of the rear bearing, the subscripts x, y, z, θ_x and θ_y represent x–y–z translational direction and rotational direction about x- and y-axes (see section “Calculation of flexible joint stiffness”) and similarly hereinafter; k_2x, k_2y, k_2z, k_2θx and k_2θy are the stiffness of the front bearing; and k_c_1x, k_c_1y, k_c_1z, k_c_1θx, k_c_1θy and k_c_1θz are the stiffness of the C-axis locking joint.

The joint of hirth coupling is defined as EJ-B. Similar to EJ-A, the left bearing, right bearing and left and right locking joints of A-axis are equivalent to EJ-C, with the stiffness written by

{\begin{matrix} {k_{q 3}}_{x} = k_{4 x} + k_{5 x} + k_{c 2 z} + k_{c 3 z} \\ {k_{q 3}}_{y} = k_{4 y} + k_{5 y} + k_{c 2 x} + k_{c 3 x} \\ {k_{q 3}}_{z} = k_{4 z} + k_{5 z} + k_{c 2 y} + k_{c 3 y} \\ {k_{q 3}}_{θ x} = {k_{c 2}}_{θ z} + {k_{c 3}}_{θ z} \\ {k_{q 3}}_{θ y} = {k_{4}}_{θ v} + {k_{5}}_{θ v} + {k_{c 2}}_{θ x} + {k_{c 3}}_{θ x} + k_{4 z} l_{3 s 1}^{2} + k_{5 z} l_{3 s 2}^{2} + k_{c 2 y} l_{3 c 1}^{2} + k_{c 3 y} l_{3 c 2}^{2} \\ {k_{q 3}}_{θ z} = {k_{4}}_{θ z} + {k_{5}}_{θ z} + {k_{c 2}}_{θ y} + {k_{c 3}}_{θ y} + k_{4 y} l_{3 s 1}^{2} + k_{5 y} l_{3 s 2}^{2} + k_{c 2 x} l_{3 c 1}^{2} + k_{c 3 x} l_{3 c 2}^{2} \end{matrix}

(2)

where k_4x, k_4y, k_4z, k_4θv and k_4θz are the stiffness of the left bearing; k_5x, k_5y, k_5z, k_5θv and k_5θz are the stiffness of the right bearing; k_c_2x, k_c_2y, k_c_2z, k_c_2θx, k_c_2θy and k_c_2θz are the stiffness of the A-axis left locking joint; and k_c_3x, k_c_3y, k_c_3z, k_c_3θx, k_c_3θy and k_c_3θz are the stiffness of the A-axis right locking joint.

As for the damping of joints, the equivalent process is similar to that of stiffness described above.

Description of the rotation and swing postures

The postures of the rotating and swing shafts can be expressed by an inertial coordinate system fixed on the C-axis housing. Using Cartesian coordinate system, the local coordinate of each rod is shown in Figure 2(c), and its origin is located in the static equilibrium position of gravity and cutting forces. In this figure, θ_A and θ_C represent the swing angle and the rotation angle, respectively. The position and posture of each point of the model can be calculated by the coordinate transformation.³⁷ Set the coordinates of any point P as P(x_i, y_i, z_i) in the moving coordinate system {O_i}, then the coordinates in the relative coordinate system {O_i₋₁} are calculated by

X_{i - 1} = {}_{i}^{i - 1}T X_{i}

(3)

where $X_{i - 1} = [\begin{matrix} x_{i - 1} & y_{i - 1} & z_{i - 1} & 1 \end{matrix}]^{T}$ , $X_{i} = [\begin{matrix} x_{i} & y_{i} & z_{i} & 1 \end{matrix}]^{T}$ and ${}_{i}^{i - 1}T$ is the coordinate transformation matrix from {O_i} to {O_i₋₁} and is expressed as equation (4) under the large displacement

{}_{i}^{i - 1}T = T' = [\begin{matrix} \cos α \cos β & \cos α \sin β \sin γ - \sin α \cos γ & \cos α \sin β \cos γ + \sin α \sin γ & x \\ \sin α \cos β & \sin α \sin β \sin γ + \cos α \cos γ & \sin α \sin β \cos γ - \cos α \sin γ & y \\ - \sin β & \cos β \sin γ & \cos β \cos γ & z \\ 0 & 0 & 0 & 1 \end{matrix}]

(4)

where α, β and γ are the z–y–x Euler angles of the origin of the moving coordinate system in the relative coordinate system; x, y and z represent the position coordinates of the origin of the moving coordinate system in the relative coordinate system.

As the vibration occurs, according to the Maclaurin formula, sin(δα) ≈ δα, cos(δα) ≈ 1 and so on, and ${}_{i}^{i - 1}T$ is given by

{}_{i}^{i - 1}T = T ″ = [\begin{matrix} 1 & - δ α & δ β & δ x \\ δ α & 1 & - δ γ & δ y \\ - δ β & δ γ & 1 & δ z \\ 0 & 0 & 0 & 1 \end{matrix}]

(5)

where δα, δβ and δγ are the micro angles of the origin of the moving coordinate system in the relative coordinate system and δx, δy and δz represent the origin’s micro displacements.

By multiplying the transformation matrix of each point, a 4 × 4 transformation matrix of P relative to fixed coordinate system {O₀} is

{}_{i}^{0}T = {}_{1}^{0}T {}_{i - 1}^{1}T {}_{1}^{i - 1}T = [\begin{matrix} {}_{i}^{0}R & {}_{i}^{0}P \\ 0 & 1 \end{matrix}]

(6)

where ${}_{i}^{0}R$ and ${}_{i}^{0}P$ are the posture matrix and position vector, respectively.

As described above, the absolute position and posture of any point of the bi-rotary milling head in the inertial coordinate system can be obtained by the rotation angle, swing angle and geometric parameters. These transformation matrices are included in Appendix 2.

Parametric dynamic equation

Based on the assumption of an ideal constrained system, the system kinetic energy T consisting of the translational kinetic energy and rotational kinetic energy is written as

T = \frac{1}{2} \sum_{i = 1}^{3} v_{i}^{T} m_{i} v_{i} + \frac{1}{2} \sum_{i = 1}^{3} ω_{i}^{T} J_{i} ω_{i}

(7)

where i is the number of the inertia; the mass matrix of the inertia i is given as

m_{i} = [\begin{matrix} m_{i} & 0 \\ m_{i} \\ 0 & m_{i} \end{matrix}]

The translational velocity is determined by

v_{i} = \frac{d}{dt} ({}_{i}^{0}P) = {[\begin{matrix} \overset{\cdot}{x} \\ \overset{\cdot}{y} \\ \overset{\cdot}{z} \end{matrix}]}_{i}

where $\overset{\cdot}{x}$ , $\overset{\cdot}{y}$ and $\overset{\cdot}{z}$ are the translational velocity in the x–y–z direction.

The moment of inertia matrix is determined by

J_{i} = {}_{i}^{0}R [\begin{matrix} I_{xx}^{i} & 0 \\ I_{yy}^{i} \\ 0 & I_{zz}^{i} \end{matrix}] {{}_{i}^{0}R}^{T}

where $I_{xx}^{i}$ , $I_{yy}^{i}$ and $I_{zz}^{i}$ are principal moments of inertia of the inertia i.

The rotational velocity is determined by

ω_{i} = ω_{i - 1} + {}_{i}^{i - 1}R ω'_{i} = {[\begin{matrix} {\overset{\cdot}{θ}}_{x} \\ {\overset{\cdot}{θ}}_{y} \\ {\overset{\cdot}{θ}}_{z} \end{matrix}]}_{i - 1} + {}_{i}^{i - 1}R {[\begin{matrix} {\overset{\cdot}{θ}}_{x} \\ {\overset{\cdot}{θ}}_{y} \\ {\overset{\cdot}{θ}}_{z} \end{matrix}]}^{'}_{i}

where $ω_{i - 1}$ is the rotational velocity vector of the origin of the moving coordinate system in the relative coordinate system and $ω'_{i}$ is the rotational velocity vector in the moving coordinate system.

Suppose that the system vibrates around the gravity balance position, and consequently only elastic potential energy is left, the potential energy V is expressed as

V = \frac{1}{2} \sum_{i = 1}^{3} x_{i}^{T} k_{i} x_{i} + \frac{1}{2} \sum_{i = 1}^{3} φ_{i}^{T} k_{θ i} φ_{i}

(8)

where $x_{i} = [\begin{matrix} x & y & z \end{matrix}]_{i}^{T} - [\begin{matrix} x & y & z \end{matrix}]_{i - 1}^{T}$ is the relative displacement vector of the inertia i, $φ_{i} = [\begin{matrix} θ_{x} & θ_{y} & θ_{z} \end{matrix}]_{i}^{T} - [\begin{matrix} θ_{x} & θ_{y} & θ_{z} \end{matrix}]_{i - 1}^{T}$ is the relative rotational angle vector of the inertia i, $k_{i} = {}_{i}^{0}R [\begin{matrix} k_{qix} & 0 \\ k_{qiy} \\ 0 & k_{qiz} \end{matrix}] {{}_{i}^{0}R}^{T}$ is the translational stiffness matrix of the ith equivalent joint and $k_{θ i} = {}_{i}^{0}R [\begin{matrix} k_{qi θ x} & 0 \\ k_{qi θ y} \\ 0 & k_{qi θ z} \end{matrix}] {{}_{i}^{0}R}^{T}$ is the rotational stiffness matrix of the ith equivalent joint.

The Rayleigh loss function is defined as

Ψ = \frac{1}{2} \sum_{i = 1}^{3} {\overset{\cdot}{x}}_{i}^{T} c_{i} {\overset{\cdot}{x}}_{i} + \frac{1}{2} \sum_{i = 1}^{3} {\overset{\cdot}{φ}}_{i}^{T} c_{θ i} {\overset{\cdot}{φ}}_{i}

(9)

where ${\overset{\cdot}{x}}_{i} = {[\begin{matrix} \overset{\cdot}{x} & \overset{\cdot}{y} & \overset{\cdot}{z} \end{matrix}]}_{i}^{T} - {[\begin{matrix} \overset{\cdot}{x} & \overset{\cdot}{y} & \overset{\cdot}{z} \end{matrix}]}_{i - 1}^{T}$ is the relative translational velocity vector of the inertia i, ${\overset{\cdot}{φ}}_{i} = {[\begin{matrix} {\overset{\cdot}{θ}}_{x} & {\overset{\cdot}{θ}}_{y} & {\overset{\cdot}{θ}}_{z} \end{matrix}]}_{i}^{T} - {[\begin{matrix} {\overset{\cdot}{θ}}_{x} & {\overset{\cdot}{θ}}_{y} & {\overset{\cdot}{θ}}_{z} \end{matrix}]}_{i - 1}^{T}$ is the relative rotational velocity vector of the inertia i, $c_{i} = {}_{i}^{0}R [\begin{matrix} c_{qix} & 0 \\ c_{qiy} \\ 0 & c_{qiz} \end{matrix}] {{}_{i}^{0}R}^{T}$ is the translational damping matrix of the ith equivalent joint and $c_{θ i} = {}_{i}^{0}R [\begin{matrix} c_{qi θ x} & 0 \\ c_{qi θ y} \\ 0 & c_{qi θ z} \end{matrix}] {{}_{i}^{0}R}^{T}$ is the rotational damping matrix of the ith equivalent joint.

Substituting equations (7)–(9) into the Lagrange equation

\frac{d}{dt} (\frac{\partial T}{\partial {\overset{\cdot}{x}}_{j}}) - \frac{\partial T}{\partial x_{j}} + \frac{\partial V}{\partial x_{j}} + \frac{\partial Ψ}{\partial {\overset{\cdot}{x}}_{j}} = F_{j}, j = 1, 2, \dots, n

(10)

where F_j is the external force in the jth DOF and n is the number of DOF.

We can obtain the dynamic equations of bi-rotary milling head. The centrifugal force and the Coriolis force generated by the rigid-body motion are ignored due to the low rotating and swing speeds. Therefore, linear vibration differential equations are derived in matrix form as

M (θ_{A}, θ_{C}) \overset{\cdot\cdot}{x} + C (θ_{C}) \overset{\cdot}{x} + K (θ_{C}) x = Q

(11)

where M is the mass matrix; C is the damping matrix; K is the stiffness matrix; Q is the external force vector; and x is the micro displacement vector, x = [x₁, y₁, z₁, θ_x₁, θ_y₁, θ_z₁, x₂, y₂, z₂, θ_x₂, θ_y₂, θ_z₂, x₃, y₃, z₃, θ_x₃, θ_y₃, θ_z₃]^T.

As seen from equation (11), the coefficients M , C and K are not constant. M is a function of θ_A and θ_C, which describes the distribution of the system mass at different postures, and K changes with θ_C. However, it can be known that K is also associated with θ_A in section “Calculation of flexible joint stiffness.” This is because the matrix elements representing contact stiffness of joints vary with different swing angles. The damping is simplified using the Rayleigh viscous equivalent damping,³⁸ given by

C = α M + β K

(12)

where α and β are proportional coefficients.

Finally, substituting equation (12) into equation (11), the pose-dependent dynamic equation of the bi-rotary milling head is expressed as

\begin{matrix} M (θ_{A}, θ_{C}) \overset{\cdot\cdot}{x} + (α M (θ_{A}, θ_{C}) \\ + β K (θ_{A}, θ_{C})) \overset{\cdot}{x} + K (θ_{A}, θ_{C}) x = Q \end{matrix}

(13)

where the mass and stiffness matrices are listed in Appendix 3.

Equation (13) is a second-order linear equation with 18 DOFs. It is convenient to calculate the system pose-dependent dynamics by equation (13) if the posture parameters and physical parameters are known.

Calculation of flexible joint stiffness

Contact stiffness of the rotating shaft bearing

The rotating shaft is usually supported by two bearings: the front one using three-row cylindrical roller slewing bearing is the main load bearing and the rear one using cylindrical roller bearing which is not a load bearing and only radial positioning. Considering the influence of gravity and cutting force, the front bearing loads are written as F ₂ = {F_2x, F_2y, F_2z, M_2x, M_2y}, as shown in Figure 3.

Figure 3.

Force diagram of the three-row roller slewing bearing.

The outer ring of the bearing is considered fixed in space as the load is applied to the bearing, and the inner ring could be displaced. The relationship between roller contact deformation and contact force meets Hertz contact deformation. Under the external force vector F ₂, force and moment balance equations of the bearing can be expressed as^39,40

{\begin{matrix} F_{2 x} = \sum_{i = 1}^{Z_{3}} Q_{{ψ_{3,}}_{i}} \cos ψ_{3, i} \\ F_{2 y} = \sum_{i = 1}^{Z_{3}} Q_{ψ_{3, i}} \sin ψ_{3, i} \\ F_{2 z} = - \sum_{i = 1}^{Z_{1}} Q_{{ψ_{1,}}_{i}} + \sum_{i = 1}^{Z_{2}} Q_{{ψ_{2,}}_{i}} \\ M_{2 x} = - \frac{D_{m 1}}{2} \sum_{i = 1}^{Z_{1}} Q_{{ψ_{1,}}_{i}} \sin ψ_{1, i} + \frac{D_{m 2}}{2} \sum_{i = 1}^{Z_{2}} Q_{{ψ_{2,}}_{i}} \sin ψ_{2, i} \\ M_{2 y} = \frac{D_{m 1}}{2} \sum_{i = 1}^{Z_{1}} Q_{{ψ_{1,}}_{i}} \cos ψ_{1, i} - \frac{D_{m 2}}{2} \sum_{i = 1}^{Z_{2}} Q_{{ψ_{2,}}_{i}} \cos ψ_{2, i} \end{matrix}

(14)

where $Q_{ψ_{1, i}}$ , $Q_{ψ_{2, i}}$ and $Q_{ψ_{3, i}}$ are the normal contact loads of the ith roller in the first, second and third rows, respectively; Z₁, Z₂ and Z₃ are the cylindrical roller numbers in the first, second and third rows, respectively; $ψ_{1, i}$ , $ψ_{2, i}$ and $ψ_{3, i}$ are the position angles in the first, second and third rows, respectively; and D_m₁ and D_m₂ are the diameters of the circle of cylindrical roller center in the first and second rows, respectively.

The external force vector F ₂ is definitely expressed as

{\begin{matrix} F_{2 x} = F_{t} \\ F_{2 y} = F_{a} \sin θ_{A} + F_{r} \cos θ_{A} \\ F_{2 z} = - (m_{1} + m_{2} + m_{3}) g + F_{r} \sin θ_{A} - F_{a} \cos θ_{A} \\ M_{2 x} = - m_{3} g l_{3 m} \sin θ_{A} + (F_{r} \sin θ_{A} - F_{a} \cos θ_{A}) l_{3} \sin θ_{A} + (F_{a} \sin θ_{A} + F_{r} \cos θ_{A}) \times (l_{1} - l_{1 m} - l_{1 s 2} + l_{2} + l_{3} \cos θ_{A}) \\ M_{2 y} = - F_{t} (l_{1} - l_{1 m} - l_{1 s 2} + l_{2} + l_{3} \cos θ_{A}) \end{matrix}

(15)

where F_t, F_a and F_r are the three components of average static cutting force, respectively, and g = 9.8 N/kg.

Substituting equation (15) into equation (14), we obtain the relationship of the force vector F ₂ and its deformation Δ₂ = {δ_2x, δ_2y, δ_2z, θ_2x, θ_2y}, then the contact stiffness is calculated by

{k_{2}}_{j} = \frac{F_{2 j}}{Δ_{2 j}}

(16)

where the subscript j represents the x–y–z–θx–θy deformation direction in turn.

The stiffness of rear bearing which is a general cylindrical roller bearing is calculated directly by the literature.⁴⁰

Contact stiffness of the coupling joint for replaceable milling head

For replaceable milling head, the hirth coupling is always used as the coupling between rotating shaft and swing shaft housing, and its loads are noted as F ₃ = {F_3x, F_3y, F_3z, M_3x, M_3y, M_3z}, as shown in Figure 4.

Figure 4.

Force diagram of the hirth coupling.

The following assumptions are made: (a) the upper and lower tooth disks of hirth coupling are undeformed except for the contact deformation of the mating surfaces between the upper and lower teeth. (b) The contact between tooth surfaces is uniform. The deformation of centroid of the tooth surface is defined as its overall deformation due to the small area. (c) The radial loads and moments act on the center plane of the hirth coupling. Based on these assumptions, Figure 5 gives the surface contact forces of single tooth as well as the definition of coordinate system.

Figure 5.

Force diagram of tooth surfaces of the hirth coupling.

The position angle of each pair of mating teeth surfaces is determined by

ψ_{i} = \frac{360 \circ}{4 Z} + \frac{360 \circ}{2 Z} (i - 1) (i = 1, 2, \dots, 2 Z)

(17)

where Z is the number of teeth.

Supposing the upper tooth disk is constrained, the lower tooth disk would produce a displacement Δ₃ described as Δ₃ = {δ_3x, δ_3y, δ_3z, θ_3x, θ_3y, θ_3z} under the external load F ₃. Accordingly, the elastic deformation of each pair of mating surfaces can be determined by

δ_{i} = {\begin{matrix} δ_{ai} \sin θ_{d} + δ_{3 x} \sin ψ_{i} - δ_{3 y} \cos ψ_{i} - R_{m} θ_{3 z} (i = 1, 3, \dots, 2 Z - 1) \\ δ_{ai} \sin θ_{d} - δ_{3 x} \sin ψ_{i} + δ_{3 y} \cos ψ_{i} + R_{m} θ_{3 z} (i = 2, 4, \dots, 2 Z) \end{matrix}

(18)

where θ_d is the pressure angle of the mating surface, R_m is the circle radius of the mating surface centroid and δ_ai is given by

δ_{ai} = δ_{3 z} + R_{m} θ_{3 x} \sin ψ_{i} - R_{m} θ_{3 y} \cos ψ_{i} (i = 1, 2, \dots, 2 Z)

(19)

For the single tooth surface, the relationship between elastic deformation δ and the average surface pressure P_n can be expressed as^41,42

δ = a {P_{n}}^{b}

(20)

where a and b are constants determined by machining methods, materials and so on.

Substituting equations (18) and (19) into equation (20), we obtain the normal contact load of each pair of the mating tooth surfaces

Q_{ψ i} = A {(\frac{δ_{i}}{a})}^{1 / b}

(21)

where A is the contact area.

The force and moment balance equations of hirth coupling are established as

{\begin{matrix} F_{3 x} = \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \cos θ_{d} \sin ψ_{2 i - 1} + f_{2 ψ_{2 i - 1}} \sin θ_{d} | \sin ψ_{2 i - 1} | + f_{1 ψ_{2 i - 1}} | \cos ψ_{2 i - 1} | - Q_{ψ_{2 i}} \cos θ_{d} \sin ψ_{2 i} + f_{2 ψ_{2 i}} \sin θ_{d} | \sin ψ_{2 i} | + f_{1 ψ_{2 i}} | \cos ψ_{2 i - 1} |) \\ F_{3 y} = - \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \cos θ_{d} \cos ψ_{2 i - 1} - f_{2 ψ_{2 i - 1}} \sin θ_{d} | \cos ψ_{2 i - 1} | - f_{1 ψ_{2 i - 1}} | \sin ψ_{2 i - 1} | - Q_{ψ_{2 i}} \cos θ_{d} \cos ψ_{2 i} - f_{2 ψ_{2 i}} \sin θ_{d} | \cos ψ_{2 i} | - f_{1 ψ_{2 i}} | \sin ψ_{2 i - 1} |) \\ F_{3 z} = \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \sin θ_{d} + f_{2 ψ_{2 i - 1}} \cos θ_{d} + Q_{ψ_{2 i}} \sin θ_{d} + f_{2 ψ_{2 i}} \cos θ_{d}) \\ M_{3 x} = R_{m} \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \sin θ_{d} \sin ψ_{2 i - 1} + f_{2 ψ_{2 i - 1}} \cos θ_{d} | \sin ψ_{2 i - 1} | + Q_{ψ_{2 i}} \sin θ_{d} \sin ψ_{2 i} + f_{2 ψ_{2 i}} \cos θ_{d} | \sin ψ_{2 i} |) \\ M_{3 y} = - R_{m} \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \sin θ_{d} \cos ψ_{2 i - 1} - f_{2 ψ_{2 i - 1}} \cos θ_{d} | \cos ψ_{2 i - 1} | + Q_{ψ_{2 i}} \sin θ_{d} \cos ψ_{2 i} - f_{2 ψ_{2 i}} \cos θ_{d} | \cos ψ_{2 i} |) \\ M_{3 z} = - R_{m} \sum_{i = 1}^{Z} (Q_{ψ_{2 i - 1}} \cos θ_{d} - f_{2 ψ_{2 i - 1}} \sin θ_{d} - Q_{ψ_{2 i}} \cos θ_{d} + f_{2 ψ_{2 i}} \sin θ_{d}) \end{matrix}

(22)

where $f_{1 ψ_{i}} = f_{2 ψ_{i}} = μ Q_{ψ_{i}}$ and µ = 0.15 is the static friction coefficient between contact surfaces.

The external force vector is

{\begin{matrix} F_{3 x} = F_{t} \\ F_{3 y} = F_{a} \sin θ_{A} + F_{r} \cos θ_{A} \\ F_{3 z} = - (m_{2} + m_{3}) g + F_{r} \sin θ_{A} - F_{a} \cos θ_{A} + F_{pre} \\ M_{3 x} = - m_{3} g l_{3 m} \sin θ_{A} + (F_{r} \sin θ_{A} - F_{a} \cos θ_{A}) l_{3} \sin θ_{A} + (F_{a} \sin θ_{A} + F_{r} \cos θ_{A}) \times (l_{2} + l_{3} \cos θ_{A}) \\ M_{3 y} = - F_{t} (l_{2} + l_{3} \cos θ_{A}) \\ M_{3 z} = F_{t} l_{3} \sin θ_{A} \end{matrix}

(23)

where F_pre is the initial preload.

The least square fitting method or Newton–Raphson method can be employed to calculate the deformation Δ₃ through equations (17)–(23) and then the contact stiffness is given by

k_{3 j} = \frac{F_{3 j}}{Δ_{3 j}}

(24)

where the subscript j represents the x–y–z–θx–θy–θz deformation direction in turn.

Contact stiffness of the swing shaft bearing

Cross-roller bearing is always used to support the swing shaft. Affected by the gravity, cutting force and swing posture, the bearing is also loaded by a 5-DOF external force F ₄ = {F_4x, F_4y, F_4z, M_4y, M_4z}, as shown in Figure 6.

Figure 6.

Force diagram of the cross-roller bearing.

Assuming that the outer ring of the bearing is fixed, then the inner ring produces a displacement as the load is applied to the bearing. The relationship between roller contact deformation and contact force meets Hertz contact deformation. The force and moment balance equations of bearing are as follows⁴³

{\begin{matrix} F_{4 x} = - \sum_{i = 1}^{Z / 2} (Q_{ψ_{2 i - 1}} \cos a_{1} - Q_{ψ_{2 i}} \cos a_{2}) \\ F_{4 y} = - \sum_{i = 1}^{Z / 2} (Q_{ψ_{2 i - 1}} \sin a_{1} \cos ψ_{2 i - 1} + Q_{ψ_{2 i}} \sin a_{2} \cos ψ_{2 i}) \\ F_{4 z} = - \sum_{i = 1}^{Z / 2} (Q_{ψ_{2 i - 1}} \sin a_{1} \sin ψ_{2 i - 1} + Q_{ψ_{2 i}} \sin a_{2} \sin ψ_{2 i}) \\ M_{4 y} = - \frac{D_{m}}{2} \sum_{i = 1}^{Z / 2} (Q_{ψ_{2 i - 1}} \cos a_{1} \sin ψ_{2 i - 1} - Q_{ψ_{2 i}} \cos a_{2} \sin ψ_{2 i}) \\ M_{4 z} = - \frac{D_{m}}{2} \sum_{i = 1}^{Z / 2} (Q_{ψ_{2 i - 1}} \cos a_{1} \cos ψ_{2 i - 1} - Q_{ψ_{2 i}} \cos a_{2} \cos ψ_{2 i}) \end{matrix}

(25)

where $Q_{ψ_{2 i - 1}}$ and $Q_{ψ_{2 i}}$ are the normal contact loads of roller in the first and second columns, respectively; a₁ and a₂ are the roller contact angles in the first and second columns, respectively; $ψ_{2 i - 1}$ and $ψ_{2 i}$ are the roller position angles in the first and second columns, respectively; Z is the number of cylindrical rollers; and D_m is the diameter of the circle of cylindrical roller center. And the external force vector F ₄ is derived as

{\begin{matrix} F_{4 x} = F_{t} \\ F_{4 y} = \frac{m_{3} g l_{3 s 2}}{(l_{3 s 1} + l_{3 s 2})} + F_{a} \cos θ_{A} - F_{r} \sin θ_{A} \\ F_{4 z} = F_{a} \sin θ_{A} + F_{r} \cos θ_{A} \\ M_{4 y} = (F_{a} \sin θ_{A} + F_{r} \cos θ_{A}) l_{3 s 1} + F_{t} l_{3} \sin θ_{A} \\ M_{4 z} = - (\frac{m_{3} g l_{3 s 2}}{(l_{3 s 1} + l_{3 s 2})} + F_{a} \cos θ_{A} - F_{r} \sin θ_{A}) l_{3 s 1} + F_{t} l_{3} \cos θ_{A} \end{matrix}

(26)

Substituting equation (26) into equation (25), we can establish the relationship of the force vector F ₄ and its deformation Δ₄ = {δ_4x, δ_4y, δ_4z, θ_4y, θ_4z}, then contact stiffness is calculated by

{k_{4}}_{j} = \frac{F_{4 j}}{Δ_{4 j}}

(27)

where the subscript j represents the x–y–z–θy–θz deformation direction in turn.

Contact stiffness of the locking joint

As mentioned above, two methods are commonly used for the positioning and locking of the rotating and swing shafts. Considering the fixed angle machining, the mechanical method is discussed here, and it is achieved by the friction of contact surface between the friction ring and the shaft, as shown in Figure 7(a). The contact surface (S_abcd) is cylindrical, as shown in Figure 7(b).

Figure 7.

Schematic diagram of cylindrical locking joint: (a) operating diagram, (b) structure diagram and (c) schematic diagram of xoy plane.

It is assumed that the contact surface of the friction ring is uniformly pressured. For cylindrical contact surface per unit area (ΔS: r × Δθ × Δz), three orthogonal contact stiffness are considered including one normal contact stiffness a_n and two tangential stiffness a_t, respectively. From equation (20), taking the derivative of P_n with respect to δ we can yield

a_{n} = \frac{d P_{n}}{d δ} = \frac{1}{ab} P_{n}^{1 - b}

(28)

The relationship between tangential contact stiffness and normal contact stiffness is⁴⁴

a_{t} = α_{nt} a_{n}

(29)

where α_nt is constant determined by materials, surface roughness and so on.

As shown in Figure 7(c), the normal and tangential contact stiffness can be transferred to the values in the Cartesian coordinates {O_xyz}, expressed as

{\begin{matrix} a_{eqx} = a_{n} \cos^{2} θ + a_{t} \sin^{2} θ \\ a_{eqy} = a_{n} \sin^{2} θ + a_{t} \cos^{2} θ \\ a_{eqz} = a_{t} \end{matrix}

(30)

where θ is the position angle in the circumferential direction.

Finally, integrating in the whole contact surface, we can obtain the total contact stiffness of the cylindrical joint as follows

{\begin{matrix} k_{cx} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} a_{eqx} rd θ dz = π rh a_{n} (1 + α_{nt}) \\ k_{cy} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} a_{eqy} rd θ dz = π rh a_{n} (1 + α_{nt}) \\ k_{cz} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} a_{eqz} rd θ dz = 2 π rh α_{nt} a_{n} \\ k_{c θ_{x}} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} (a_{eqz} r^{2} \sin^{2} θ + a_{eqy} z^{2}) d θ dz = π h a_{n} (α_{nt} r^{2} + \frac{(α_{nt} + 1) h^{2}}{12}) \\ k_{c θ_{y}} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} (a_{eqz} r^{2} \cos^{2} θ + a_{eqx} z^{2}) d θ dz = π h a_{n} (α_{nt} r^{2} + \frac{(α_{nt} + 1) h^{2}}{12}) \\ k_{c θ_{z}} = \int_{- h / 2}^{h / 2} \int_{0}^{2 π} a_{t} r^{3} d θ dz = 2 π r^{3} h α_{nt} a_{n} \end{matrix}

(31)

where h and r are the height and radius of the cylindrical joint, respectively.

Substituting the parameters of A/C-axis locking joint into equation (31), we can obtain the stiffness of the C-axis locking joint k _c1, the stiffness of the A-axis left locking joint k _c2 and the stiffness of the A-axis right locking joint k _c3.

Pose-dependent dynamics

Static stiffness at different postures

The dynamics of a bi-rotary milling head shown in Figure 2 within its workspace were calculated by the parametric equation (13). The simulation parameters are included in Appendix 4. The static stiffness maps of the milling head at the spindle nose is shown in Figure 8, where it can be seen that the translational and rotational stiffness are symmetric about θ_A = 0°. The values of K_x, K_y, K_θx and K_θz reach the minimum at θ_A = 0°, while their maximum are at θ_A = 90°. But the minimum and maximum values of K_z and K_θy are opposite.

Figure 8.

Static stiffness maps at the spindle nose.

Natural frequencies at different postures

Figure 9 shows the variation of natural frequencies during the whole working postures. As it can be seen, the rotation angle has little effect on the system natural frequencies, and the first two natural frequencies at positive and negative swing angles are almost the same due to the structural symmetry. However, it is a remarkable fact that the two natural frequencies both increase with θ_A.

Figure 9.

Pose-dependent natural frequencies of the bi-rotary milling head: (a) first natural frequency and (b) second natural frequency.

FRFs at different swing angles

The FRFs at the spindle nose within swing angles from 0° to 90° are shown in Figure 10. The first and second natural frequencies increase with the swing angle which is consistent with the discussion in the previous section. The amplitude at the first natural frequency slowly increases with θ_A; however, the amplitude at the second natural frequency decreases significantly with the increase in θ_A.

Figure 10.

FRFs of simulation at the spindle nose at different swing angles: (a) YY FRFs and (b) XX FRFs.

Reasons for pose-dependent dynamics

From Figures 8 –10, the dynamics of the bi-rotary milling head are pose dependent and sensitive to the swing angle. The variation of the swing angle can change the system moments of inertia according to equation (13), thereby changing its mass matrix, as shown in Figure 11. Moreover, the stiffness matrix is also changed owing to the nonlinear change of flexible joints at different swing angles.

Figure 11.

Influence of the swing angle on the mass matrix.

The eccentric mass m₃ would produce a moment M_2x which can be calculated by equation (15) when the swing shaft swings without other loads, as plotted in Figure 12. It can be clearly seen that the gravity moment varies with the swing angle. Subjected to the varying moment, the equivalent contact stiffness of the front bearing at different swing angles is calculated by equations (14)–(16), as shown in Figure 13(a). It indicates that k_2θx varies obviously with the swing angle. For the joint of the hirth coupling, as shown in Figure 13(b), the contact stiffness calculated by equations (17)–(24) is also varying.

Figure 12.

Variation of gravity moment with the swing angle.

Figure 13.

Influence of the swing angle on the contact stiffness: (a) the joint of three-row cylindrical roller bearings and (b) the joint of hirth coupling.

Figure 14 describes the influence of the varying moment of inertia and stiffness on natural frequencies with the swing angle. For the first natural frequency, plotted in Figure 14(a), the change of mass matrix Δ M and stiffness matrix Δ K both cause an increase in natural frequencies with the swing angle except for the range of 70°–90°, which is the result of a sharp decline of the contact stiffness k_2θx (see Figure 13(a)). Moreover, Δ K is more important than Δ M when the swing angle is <90°. While for the second natural frequency, plotted in Figure 14(b), Δ K leads to a decrease with the increase in the swing angle instead. However, when combining the influence of Δ M and Δ K , the predicted values agree better with measurements for either of the two frequencies.

Figure 14.

Influence of the swing angle on the natural frequency: (a) first natural frequency and (b) second natural frequency.

Experimental verification

Experimental setup

A fork-type direct-drive bi-rotary milling head with A/C-axis (θ_A ∈ −110° to 110°; θ_C ∈ 0°–360°) was employed as the test object, as shown in Figure 15, where the front bearing adopts three-row roller slewing bearing (model: YRT 325), the rear one is cylindrical roller bearing (model: N 1922) and both the left and the right bearings use cross-roller bearings (model: SX 011824) with a supporting span of 368 mm. The circumferential hydraulic brake with a pressure of 4 MPa is used for the locking of A/C-axis. The maximum torque and the highest speed of C-axis motor are 2394 N m and 60 r/min, while for A-axis motor, the values are 1450 N m and 60 r/min, respectively.

Figure 15.

Impact testing experiment.

Modal testing was performed by LMS Test.Lab with a hammer (model: PCB 086D05, sensitivity: 0.21 mV/N) exciting on the spindle nose and three-direction acceleration sensors (model: PCB 352A16, sensitivity: 100 mV/g) used to pick acceleration response of the measuring points distributed on the milling head. The number of measurement points is 40, the frequency bandwidth is 512 Hz and the number of spectral line is 2048. And then the modal frequencies, mode shapes and FRFs were calculated by model analysis module of the test system. The tests were repeated every 20° within swing angles from 0° to 100°.

Comparison between experimental and simulated results

Test results of natural frequencies at different swing angles are plotted in Figure 14. The first two natural frequencies both increase with the swing angle, and the variation is consistent with the evaluated results. Moreover, the evaluated value at each swing angle is very close to the experimental value. The mode shape of the first frequency is a first-order swing mode around the x-axis (Figure 16(a)), and the second mode shape is a first-order swing mode around the y-axis (Figure 16(b)). Figure 16 shows that the low-frequency modes of the milling head are the relative vibration between structural components, so the rigid-body assumption is tenable.

Figure 16.

Mode shape comparison between experimental and simulated results: (a) first mode (126 Hz) and (b) second mode (131 Hz).

The measured FRFs at different swing angles are shown in Figure 17. The change trend of the first two natural frequencies and corresponding amplitudes are also close to the evaluated values (Figure 10) although there exists a little discrepancy at the swing angles of 80° and 90°. The different amplitudes at the resonance peaks between the simulation and the experiment are mainly affected by the damping. The Rayleigh viscous equivalent damping is used in the model, while the damping of the actual system is a nonviscous damping. In addition, the residual modes and the neglect of the cross stiffness of the joint are also the factors causing the error. Compared with the tested natural frequencies and FRFs, the simulated results can match well with the experimental results, which validate the modeling method proposed above.

Figure 17.

FRFs of experiment at the spindle nose at different swing angles: (a) YY FRFs and (b) XX FRFs.

Conclusion

In this article, a parametric modeling method is presented to evaluate the pose-dependent dynamics of bi-rotary milling heads. The explicit expressions of the system mass, damping and stiffness matrices with rotation and swing angles are obtained in a straightforward way, which is conducive to elegant and efficient simulation and analysis of the dynamic behavior within the whole workspace. Moreover, the pose-varying stiffness of joints is considered. Using this method, a fork-type direct-drive bi-rotary milling head is analyzed, and the main conclusions are as follows:

The predicted and measured results indicate that the system dynamics are pose dependent in terms of static stiffness maps, natural frequencies and FRFs at different postures. The reasons are not only the variation of the mass distribution of the system but also the varying stiffness of three-row cylindrical roller slewing bearing and hirth coupling affected by the gravity at different swing angles.

The pose-dependent dynamics are sensitive to the swing angle while nonsensitive to the rotation angle. It is because that the rotating angle does not change the relative position between the rotating shaft and swing shaft, and the joint stiffness is a constant at a certain swing angle.

The parametric modeling method can also be used for other milling heads with different structures and transmission forms to facilitate rapid assessment of design alternatives at an early stage of development. Furthermore, the dynamic model of tool–spindle–milling head system can be established, which helps to explore the effect of pose-dependent dynamics of the milling head on cutting stability to improve machining accuracy and efficiency of five-axis machine tools.

Footnotes

Appendix 1

Appendix 2

Table 1.

Transformation matrix of each point on the bi-rotary milling head.

Parameters	Type	α	β	γ	x	y	z
${}_{1}^{0}T$	$T ″$	θ_z ₁	θ_y ₁	θ_x1	x ₁	y ₁	z ₁
${}_{m_{1}}^{1}T$	$T'$	θ_C	0	0	0	0	0
${}_{1_{1}}^{1}T$	$T'$	θ_C	0	0	0	0	l _1m
${}_{1_{2}}^{1}T$	$T'$	θ_C	0	0	0	0	−l₁+l_1m
${}_{2}^{0}T$	$T ″$	θ_z2	θ_y ₂	θ_x ₂	x ₂	y ₂	z ₂
${}_{2_{1}}^{2}T$	$T'$	θ_C	0	0	0	0	0
${}_{m_{2}}^{2}T$	$T'$	θ_C	0	0	0	0	−l_2m
${}_{2_{2}}^{2}T$	$T'$	θ_C	0	0	0	0	−l₂
${}_{3}^{0}T$	$T ″$	θ_z ₃	θ_y ₃	θ_x ₃	x ₃	y ₃	z ₃
${}_{3_{1}}^{3}T$	$T'$	θ_C	0	−π/2	0	0	0
${}_{m_{3}}^{3_{1}}T$	$T'$	0	0	θ_A	0	l_3mcos θ_A	l_3msin θ_A
${}_{3_{1}}^{3_{2}}T$	$T'$	0	0	θ_A	0	l₃cos θ_A	l₃sin θ_A

Appendix 3

The system mass matrix M is

(32)

M = [\begin{matrix} m_{1} \\ 0 & m_{1} \\ 0 & 0 & m_{1} \\ 0 & 0 & 0 & a_{1} \\ 0 & 0 & 0 & a_{2} & a_{3} & S \\ 0 & 0 & 0 & 0 & 0 & I_{z z}^{1} & Y \\ 0 & 0 & 0 & 0 & 0 & 0 & m_{2} & W \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{2} l_{2 m} & 0 & a_{4} \\ 0 & 0 & 0 & 0 & 0 & 0 & - m_{2} l_{2 m} & 0 & 0 & a_{5} & a_{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & I_{z z}^{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m_{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{9} & a_{11} & a_{13} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{7} & 0 & a_{12} & a_{14} & a_{16} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{8} & a_{10} & 0 & a_{15} & a_{17} & a_{18} \end{matrix}]

where

\begin{matrix} a_{1} = I_{xx}^{1} \cos^{2} θ_{C} + I_{yy}^{1} \sin^{2} θ_{C} \\ a_{2} = I_{xx}^{1} \cos θ_{C} \sin θ_{C} - I_{yy}^{1} \cos θ_{C} \sin θ_{C} \\ a_{3} = I_{xx}^{1} \sin^{2} θ_{C} + I_{yy}^{1} \cos^{2} θ_{C} \\ a_{4} = I_{xx}^{2} \cos^{2} θ_{C} + I_{yy}^{2} \sin^{2} θ_{C} + m_{2} l_{2 m}^{2} \\ a_{5} = I_{xx}^{2} \cos θ_{C} \sin θ_{C} - I_{yy}^{2} \cos θ_{C} \sin θ_{C} \\ a_{6} = I_{xx}^{2} \sin^{2} θ_{C} + I_{yy}^{2} \cos^{2} θ_{C} + m_{2} l_{2 m}^{2} \\ a_{7} = - m_{3} l_{3 m} \cos θ_{A} \\ a_{8} = - m_{3} l_{3 m} \sin θ_{A} \cos θ_{C} \\ a_{9} = m_{3} l_{3 m} \cos θ_{A} \\ a_{10} = - m_{3} l_{3 m} \sin θ_{A} \sin θ_{C} \\ a_{11} = m_{3} l_{3 m} \sin θ_{A} \cos θ_{C} \\ a_{12} = m_{3} l_{3 m} \sin θ_{A} \sin θ_{C} \\ a_{13} = I_{xx}^{3} \cos^{2} θ_{C} + I_{yy}^{3} \sin^{2} θ_{A} \sin^{2} θ_{C} + I_{zz}^{3} \cos^{2} θ_{A} \sin^{2} θ_{C} + m_{3} l_{3 m}^{2} \cos^{2} θ_{C} \sin^{2} θ_{A} + m_{3} l_{3 m}^{2} \cos^{2} θ_{A} \\ a_{14} = I_{xx}^{3} \cos θ_{C} \sin θ_{C} - I_{yy}^{3} \sin θ_{C} \cos θ_{C} \sin^{2} θ_{A} - I_{zz}^{3} \cos^{2} θ_{A} \sin θ_{C} \cos θ_{C} + m_{3} l_{3 m}^{2} \sin θ_{C} \cos θ_{C} \sin^{2} θ_{A} \\ a_{15} = I_{yy}^{3} \sin θ_{A} \sin θ_{C} \cos θ_{A} - I_{zz}^{3} \sin θ_{A} \sin θ_{C} \cos θ_{A} - m_{3} l_{3 m}^{2} \cos θ_{A} \sin θ_{A} \sin θ_{C} \\ a_{16} = I_{xx}^{3} \sin^{2} θ_{C} + I_{yy}^{3} \sin^{2} θ_{A} \cos^{2} θ_{C} + I_{zz}^{3} \cos^{2} θ_{A} \cos^{2} θ_{C} + m_{3} l_{3 m}^{2} \cos^{2} θ_{C} \cos^{2} θ_{A} + m_{3} l_{3 m}^{2} \sin^{2} θ_{C} \\ a_{17} = - I_{yy}^{3} \sin θ_{A} \cos θ_{C} \cos θ_{A} + I_{zz}^{3} \sin θ_{A} \cos θ_{C} \cos θ_{A} + m_{3} l_{3 m}^{2} \sin θ_{A} \cos θ_{A} \cos θ_{C} \\ a_{18} = I_{yy}^{3} \cos^{2} θ_{A} + I_{zz}^{3} \sin^{2} θ_{A} + m_{3} l_{3 m}^{2} \sin^{2} θ_{A} \end{matrix}

The system stiffness matrix K is

(33)

K = [\begin{matrix} b_{1} \\ b_{2} & b_{5} \\ 0 & 0 & b_{8} \\ 0 & 0 & 0 & b_{9} \\ 0 & 0 & 0 & b_{10} & b_{13} & S \\ 0 & 0 & 0 & 0 & 0 & b_{16} & Y \\ b_{3} & b_{6} & 0 & 0 & 0 & 0 & b_{17} & M \\ b_{4} & b_{7} & 0 & 0 & 0 & 0 & b_{18} & b_{21} \\ 0 & 0 & - k_{3 z} & 0 & 0 & 0 & 0 & 0 & b_{24} \\ 0 & 0 & 0 & b_{11} & b_{14} & 0 & 0 & 0 & 0 & b_{25} \\ 0 & 0 & 0 & b_{12} & b_{15} & 0 & 0 & 0 & 0 & b_{26} & b_{29} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b_{32} \\ 0 & 0 & 0 & 0 & 0 & 0 & b_{19} & b_{22} & 0 & 0 & 0 & 0 & b_{33} \\ 0 & 0 & 0 & 0 & 0 & 0 & b_{20} & b_{23} & 0 & 0 & 0 & 0 & b_{34} & b_{35} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - k_{q 3 y} & 0 & 0 & 0 & 0 & 0 & k_{q 3 y} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b_{27} & b_{30} & 0 & 0 & 0 & 0 & b_{36} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b_{28} & b_{31} & 0 & 0 & 0 & 0 & b_{37} & b_{38} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - k_{q 3 θ y} & 0 & 0 & 0 & 0 & 0 & k_{q 3 θ y} \end{matrix}]

where

\begin{matrix} b_{1} = {k_{3}}_{x} \cos^{2} θ_{C} + {k_{3}}_{y} \sin^{2} θ_{C} + k_{q 1 v} \sin^{2} θ_{C} + k_{q 1 x} \cos^{2} θ_{C} \\ b_{2} = - \cos θ_{C} \sin θ_{C} (k_{q 1 y} - k_{q 1 x} - {k_{3}}_{x} + {k_{3}}_{y}) \\ b_{3} = - {k_{3}}_{x} \cos^{2} θ_{C} - {k_{3}}_{y} \sin^{2} θ_{C} \\ b_{4} = - \cos θ_{C} \sin θ_{C} ({k_{3}}_{x} - {k_{3}}_{y}) \\ b_{5} = {k_{3}}_{x} \sin^{2} θ_{C} + {k_{3}}_{y} \cos^{2} θ_{C} + k_{q 1 x} \sin^{2} θ_{C} + k_{q 1 y} \cos^{2} θ_{C} \\ b_{6} = - \cos θ_{C} \sin θ_{C} ({k_{3}}_{x} - {k_{3}}_{y}) \\ b_{7} = - {k_{3}}_{x} \sin^{2} θ_{C} - {k_{3}}_{y} \cos^{2} θ_{C} \\ b_{8} = {k_{3}}_{z} + k_{q 1 z} \\ b_{9} = k_{3 θ x} \cos^{2} θ_{C} + k_{3 θ y} \sin^{2} θ_{C} + k_{q 1 θ x} \cos^{2} θ_{C} + k_{q 1 θ y} \sin^{2} θ_{C} \\ b_{10} = \cos θ_{C} \sin θ_{C} (k_{3 θ x} - k_{3 θ y} + k_{q 1 θ x} - k_{q 1 θ y}) \\ b_{11} = - k_{3 θ y} \sin^{2} θ_{C} - k_{3 θ x} \cos^{2} θ_{C} \\ b_{12} = - \cos θ_{C} \sin θ_{C} (- k_{3 θ y} + k_{3 θ x}) \\ b_{13} = k_{q 1 θ x} \sin^{2} θ_{C} + k_{3 θ x} \sin^{2} θ_{C} + k_{3 θ y} \cos^{2} θ_{C} + k_{q 1 θ y} \cos^{2} θ_{C} \\ b_{14} = - \cos θ_{C} \sin θ_{C} (k_{3 θ x} - k_{3 θ y}) \\ b_{15} = - k_{3 θ y} \cos^{2} θ_{C} - k_{3 θ x} + k_{3 θ x} \cos^{2} θ_{C} \\ b_{16} = k_{3 θ z} + k_{q 1 θ z} \\ b_{17} = {k_{3}}_{y} \sin^{2} θ_{C} + k_{q 3 x} \cos^{2} θ_{C} + k_{q 3 z} \sin^{2} θ_{C} + k_{3 x} \cos^{2} θ_{C} \\ b_{18} = - \cos θ_{C} \sin θ_{C} ({k_{3}}_{y} + k_{q 3 z} - k_{q 3 x} - {k_{3}}_{x}) \\ b_{19} = - k_{q 3 x} \cos^{2} θ_{C} - k_{q 3 z} \sin^{2} θ_{C} \\ b_{20} = \cos θ_{C} \sin θ_{C} (- k_{q 3 x} + k_{q 3 z}) \\ b_{21} = {k_{3}}_{x} \sin^{2} θ_{C} + {k_{3}}_{y} \cos^{2} θ_{C} + k_{q 3 x} \sin^{2} θ_{C} + k_{q 3 z} \cos^{2} θ_{C} \\ b_{22} = \cos θ_{C} \sin θ_{C} (- k_{q 3 x} + k_{q 3 z}) \\ b_{23} = - k_{q 3 x} \sin^{2} θ_{C} - k_{q 3 z} \cos^{2} θ_{C} \\ b_{24} = {k_{3}}_{z} + k_{q 3 y} \\ b_{25} = k_{3 θ x} \cos^{2} θ_{C} + k_{3 θ y} \sin^{2} θ_{C} + {k_{q 3 θ}}_{x} \cos^{2} θ_{C} + k_{q 3 θ z} \sin^{2} θ_{C} \\ b_{26} = \cos θ_{C} \sin θ_{C} (k_{q 3 θ x} + k_{3 θ x} - k_{3 θ y} - k_{q 3 θ z}) \\ b_{27} = - k_{q 3 θ x} \cos^{2} θ_{C} + k_{q 3 θ z} \cos^{2} θ_{C} - k_{q 3 θ z} \\ b_{28} = - \cos θ_{C} \sin θ_{C} (k_{q 3 θ x} - k_{q 3 θ z}) \\ b_{29} = k_{3 θ x} \sin^{2} θ_{C} + k_{3 θ y} \cos^{2} θ_{C} + k_{q 3 θ x} \sin^{2} θ_{C} + k_{q 3 θ z} \cos^{2} θ_{C} \\ b_{30} = - \cos θ_{C} \sin θ_{C} (k_{q 3 θ x} - k_{q 3 θ z}) \\ b_{31} = k_{q 3 θ x} \cos^{2} θ_{C} - k_{q 3 θ z} \cos^{2} θ_{C} - k_{q 3 θ x} \\ b_{32} = k_{3 θ z} + k_{q 3 θ y} \\ b_{33} = k_{q 3 z} \sin^{2} θ_{C} + k_{q 3 x} \cos^{2} θ_{C} \\ b_{34} = - \cos θ_{C} \sin θ_{C} (- k_{q 3 x} + k_{q 3 z}) \\ b_{35} = k_{q 3 x} \sin^{2} θ_{C} + k_{q 3 z} \cos^{2} θ_{C} \\ b_{36} = k_{q 3 θ x} \cos^{2} θ_{C} + k_{q 3 θ z} \sin^{2} θ_{C} \\ b_{37} = - \cos θ_{C} \sin θ_{C} (- k_{q 3 θ x} + k_{q 3 θ z}) \\ b_{38} = k_{q 3 θ z} \cos^{2} θ_{C} + k_{q 3 θ x} \sin^{2} θ_{C} \end{matrix}

Appendix 4

Parameters for the dynamic model of the bi-rotary milling head are presented in the following.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the key program of the National Natural Science Foundation (No. 51235009), the National Science and Technology Major Project of China (No. 2011zx04016-031) and the Fundamental Research Funds for the Central Universities and the Natural Science Basic Research Plan in Shaanxi Province of China (No. 2015JM5220).

References

Shaowei

Guofu

Shengfeng

. Integrated geometric error modeling, identification and compensation of CNC machine tools. Int J Mach Tool Manu 2012; 52: 24–29.

Lin

M-T

S-K

. Modeling and improvement of dynamic contour errors for five-axis machine tools under synchronous measuring paths. Int J Mach Tool Manu 2013; 72: 58–72.

Ming-Tzong

Shih-Kai

. Modeling and analysis of servo dynamics errors on measuring paths of five-axis machine tools. Int J Mach Tool Manu 2013; 66: 1–14.

Yang

Altintas

. Generalized kinematics of five-axis serial machines with non-singular tool path generation. Int J Mach Tool Manu 2013; 75: 119–132.

Yang

Peng

. Direct-drive bi-rotary milling head variable load thermal characteristics analysis. AASRI Procedia 2012; 3: 270–276.

Dhupia

Powalka

Katz

. Dynamics of the arch-type reconfigurable machine tool. Int J Mach Tool Manu 2007; 47: 326–334.

Bonnemains

Chanal

Bouzgarrou

B-C

. Stiffness computation and identification of parallel kinematic machine tools. J Manuf Sci Eng 2009; 131: 041013.

Wang

. Study on the stiffness of a 5-DOF hybrid machine tool with actuation redundancy. Mech Mach Theory 2009; 44: 289–305.

Law

Phani

Altintas

. Position-dependent multibody dynamic modeling of machine tools based on improved reduced order models. J Manuf Sci Eng 2013; 135: 021008–021008; 1-11.

10.

Albertelli

Cau

Bianchi

. The effects of dynamic interaction between machine tool subsystems on cutting process stability. Int J Adv Manuf Tech 2012; 58: 923–932.

11.

Law

Ihlenfeldt

Wabner

. Position-dependent dynamics and stability of serial-parallel kinematic machines. CIRP Ann: Manuf Techn 2013; 62: 375–378.

12.

Maj

Modica

Bianchi

. Machine tools mechatronic analysis. Proc IMechE, Part B: J Engineering Manufacture 2006; 220: 345–353.

13.

Hungsun

Hae-Jin

Hyung Wook

. Design and dynamic analysis of an arch-type desktop reconfigurable machine. Int J Mach Tool Manu 2010; 50: 575–584.

14.

Rizk

Fauroux

Mumteanu

. A comparative stiffness analysis of a reconfigurable parallel machine with three or four degrees of mobility. J Mach Eng 2006; 6: 45–55.

15.

Henninger

Eberhard

. Analysis of dynamic stability for milling processes with varying workpiece dynamics. Proc Appl Math Mech 2008; 8: 10367–10368.

16.

Zaeh

Siedl

. A new method for simulation of machining performance by integrating finite element and multi-body simulation for machine tools. CIRP Ann: Manuf Techn 2007; 56: 383–386.

17.

Fleischer

Munzinger

Tröndle

. Simulation and optimization of complete mechanical behaviour of machine tools. Prod Engineer 2008; 2: 85–90.

18.

Gao

X-S

Zhang

Y-D

Gao

H-D

. Dynamic characteristic analysis of whole machine tools based on Kriging model. J Cent South Univ 2013; 20: 3094–3102.

19.

Xiangsheng

Zhang

. Effects of machine tool configuration on its dynamics based on orthogonal experiment method. Chinese J Aeronaut 2012; 25: 285–291.

20.

Law

Altintas

Phani

. Rapid evaluation and optimization of machine tools with position-dependent stability. Int J Mach Tool Manu 2013; 68: 81–90.

21.

Luo

H-W

Wang

Zhang

. Rapid evaluation for position-dependent dynamics of a 3-DOF PKM module. Adv Mech Eng 2014; 6: 238928.

22.

Luo

Pan

Cai

. A method to predict position-dependent structural natural frequencies of machine tool. Int J Mach Tool Manu 2015; 92: 72–84.

23.

Cao

. Chatter stability of milling with speed-varying dynamics of spindles. Int J Mach Tool Manu 2012; 52: 50–58.

24.

Yigit

Ulsoy

. Dynamic stiffness evaluation for reconfigurable machine tools including weakly non-linear joint characteristics. Proc IMechE, Part B: J Engineering Manufacture 2002; 216: 87–101.

25.

Dhupia

Powalka

Ulsoy

. Effect of a nonlinear joint on the dynamic performance of a machine tool. J Manuf Sci Eng 2007; 129: 943–950.

26.

Zhang

Liu

. Dynamic modeling and analysis of the high-speed ball screw feed system. Proc IMechE, Part B: J Engineering Manufacture 2015; 229: 870–877.

27.

Wei

Zhang

. Effect of tilting angle on the dynamics of tilting table driven by worm and worm wheel. Proc IMechE, Part C: J Mechanical Engineering Science 2015; 229: 1782–1791.

28.

Wang

Liu

Yang

. The effect of axis coupling on machine tool dynamics determined by tool deviation. Int J Mach Tool Manu 2015; 88: 71–81.

29.

Altintaş

Budak

. Analytical prediction of stability lobes in milling. CIRP Ann: Manuf Techn 1995; 44: 357–362.

30.

Sulitka

Kolar

Janota

. Properties of the spindle and machine frame coupled system. Jpn Soc Mech Eng 2009; 5: 745–750.

31.

Kolar

Sulitka

Janota

. Simulation of dynamic properties of a spindle and tool system coupled with a machine tool frame. Int J Mach Tool Manu 2011; 54: 11–20.

32.

Hung

Chen

Luo

. Effect of tool orientation on the machining stability of a milling machine with swinging head. World Acad Sci Eng Tech 2013; 77: 958–964.

33.

Law

Grossi

Scippa

. Modeling the orientation-dependent dynamics of machine tools with gimbal heads. In: Proceedings of the 3rd international Chemnitz manufacturing colloquium (ICMC 2014), Chemnitz, 8–9 April 2014, pp.1–16. Fraunhofer IWU.

34.

Yan

Tang

Zheng

. Modal analysis about direct-drive bi-rotary milling head. In: Xu

(eds) Advanced materials research. Zürich: Trans Tech Publications, 2011, pp.299–303.

35.

Shi

Huang

Chen

. Finite element static stiffness analysis of the A/C axes bi-rotary milling head. In: Xu

(eds) Advanced materials research. Zürich: Trans Tech Publications, 2013, pp.1195–1199.

36.

Gomand

Bearee

Kestelyn

. Physical dynamic modelling and systematic control structure design of a double linear drive moving gantry stage industrial robot. In: 2007 European conference on power electronics and applications, Aalborg, 2–5 September 2007, pp.1–9. New York: IEEE.

37.

Craig

. Introduction to robotics: mechanics and control. Upper Saddle River, NJ: Pearson Prentice Hall, 2005.

38.

Rao

. Mechanical vibrations. Boston, MA: Addison-Wesley Reading, 2004.

39.

Aguirrebeitia

Abasolo

Avilés

. Theoretical calculation of general static load-carrying capacity for the design and selection of three row roller slewing bearings. Mech Mach Theory 2012; 48: 52–61.

40.

Harris

Kotzalas

. Essential concepts of bearing technology. Boca Raton, FL: CRC Press, 2006.

41.

Huang

Zhang

. Experimental investigation of dynamic normal characteristics of machined joint surfaces. J Vib Acoust 2000; 122: 393.

42.

Wang

Liu

Zhang

. Analysis and modeling for flexible joint interfaces under micro and macro scale. Precis Eng 2013; 37: 817–824.

43.

Harris

Kotzalas

. Advanced concepts of bearing technology: rolling bearing analysis. Boca Raton, FL: CRC Press, 2006.

44.

Gonzalez-Valadez

Baltazar

Dwyer-Joyce

. Study of interfacial stiffness ratio of a rough surface in contact using a spring model. Wear 2010; 268: 373–379.

A parametric modeling method for the pose-dependent dynamics of bi-rotary milling head

Abstract

Keywords

Introduction

Dynamic equation of milling head with rotation and swing angles

Equivalent dynamic model

Description of the rotation and swing postures

Parametric dynamic equation

Calculation of flexible joint stiffness

Contact stiffness of the rotating shaft bearing

Contact stiffness of the coupling joint for replaceable milling head

Contact stiffness of the swing shaft bearing

Contact stiffness of the locking joint

Pose-dependent dynamics

Static stiffness at different postures

Natural frequencies at different postures

FRFs at different swing angles

Reasons for pose-dependent dynamics

Experimental verification

Experimental setup

Comparison between experimental and simulated results

Conclusion

Footnotes

Appendix 1

Appendix 2

Appendix 3

Appendix 4

Declaration of conflicting interests

Funding

References