Design of eddy current dampers for vibration suppression in robotic milling

Abstract

The milling robot normally has a low stiffness which may easily cause chatter during machining. This article presents a novel eddy current damper design for chatter suppression in the robotic milling process. The designed eddy current dampers are installed on the milling spindle to damp the tool tip vibrations. The structural design and the working principle of the eddy current dampers are explained. The magnetic flux density distribution and the magnetic force generation of the designed eddy current damper are analyzed with the finite element method. The tool tip dynamics without and with eddy current dampers are modeled, and the damping performance of the proposed eddy current dampers in the robotic milling process is verified through both simulations and experiments. The results show that the peaks of the tool tip frequency response function caused by the milling tool modes are damped significantly, and the stable depth of cut is improved greatly with eddy current dampers.

Introduction

Industrial robot is a promising technology which has been implemented in the grinding, wielding, and spray-painting processes in industry. Comparing to multi-axis machine tools, industrial robot has much lower price and better flexibility, and it is able to work in parallel. These advantages make robotic milling a new choice for machining of large complex aerospace structure parts. However, due to its serial structure, the industrial robot has its natural limitations: (1) low accuracy (about 0.3 mm) ¹ and (2) low stiffness (normally less than 1 N/μm). ² Researchers have proposed different methods^3–5 to compensate for the positioning error of the milling robots under large machining load and to improve the accuracy of robotic milling significantly. But for the chatter vibration caused by the low stiffness of the milling robot, it is still a tough problem that needs to be solved.

In machining applications with machine tools, researchers have proposed different methods for chatter suppression. Altintas^6,7 proposed a spindle speed selection method in which the chatter-free axial depth of cut and spindle speed are properly selected based on the chatter stability lobes that are calculated with the tool tip frequency response function (FRF). Tarng et al. ⁸ used tuned vibration absorber to suppress the chatter vibration in the turning operation. The natural frequency of the vibration absorber was tuned to be equal to the natural frequency of the cutting tool. Yang et al. ⁹ presented an optimal tuning method for multiple tuned mass dampers to increase chatter stability. Munoa et al. ¹⁰ installed a two-axis inertial actuator on a ram-type traveling column milling machine for chatter suppression. Different control laws were implemented for chatter suppression. Van Dijk ¹¹ developed an active chatter control methodology for high-speed milling process with magnetic bearings, which can guarantee chatter-free cutting operations in an a priori defined range of process parameters such as spindle speed and depth of cut.

Although the above methods may suppress chatter in machine tool applications, they may not work well in the robotic milling case, since the tool tip dynamics of the milling robot changes at different robot postures. The stability lobes can hardly be calculated in advance to select the chatter-free cutting parameters. The tuned mass dampers may not work as well since they require accurate tuning of their natural frequencies. The performance of the model-based active chatter control methods will also be depressed due to the changing dynamics of the milling robot.

In order to improve the chatter stability in the robotic milling process, Pan et al.^2,12 studied the chatter mechanism in the milling process with an industrial robot and suggested that proper tool path and feed direction should be carefully planned in advance, so that the cutting force direction can be properly adjusted to avoid the mode coupling chatter during the robotic milling process. Zaghbani et al. ¹³ proposed a spindle speed selection method for robotic milling to achieve better vibration suppression performance. The average value of the root mean square (RMS_directional) of the vibration signals at the spindle or the workpiece location in the x-, y-, and z-directions is used as an indicator for defining the degree of stability of the machining process. The lower RMS_directional value represents better machining stability. Before the machining process, cutting tests at different spindle speeds are carried out, and the spindle speed which has the lowest RMS_directional value will be implemented for robotic milling. Guo et al. ¹⁴ presented an optimization method for robot postures in machining applications, which can be used to find the robot posture that has the largest stiffness at a certain tool center location. By optimizing the robot postures at different tool center locations along the machining path, the robot stiffness and the chatter stability can be improved. The methods proposed in previous works^2,12–14 avoid chatter by properly selecting the feed direction, spindle speed, and robot postures. However, these methods require large amounts of experiments or calculations prior to each machining process, and therefore they are time consuming.

Vibration suppression by improving the damping properties of the structures may be a better way, and this can be realized by attaching eddy current dampers (ECDs) onto the structure that need to be damped. Bae and colleagues^15–20 studied the theory of eddy current damping and verified its effectiveness experimentally on vibration suppression of flexible beams. Yang et al.^21,22 designed an ECD for vibration suppression of a thin-walled workpiece and achieved promising vibration reduction performance. From the previous studies, it can be found that ECDs can increase the damping property of the structures they are attached to in a broader frequency range over the classic tuned mass dampers. Also, they are not quite sensitive to the change of structural modal frequencies, thus having good robustness. Moreover, ECDs are passive dampers and do not require complex control laws, and therefore they are easy to be implemented. These advantages make ECDs a good choice for vibration attenuation in the robotic milling process.

In this work, a novel ECD is designed for vibration suppression in the robotic milling process. The ECDs are installed onto a milling spindle and tuned to damp the peaks of the tool tip FRFs which are caused by the milling tool modes. The damping performance of the presented ECDs in the robotic milling process is verified through both simulations and experiments. The rest of the article is organized as follows: the proposed vibration suppression method for robotic milling and the structural design of the proposed ECD are introduced in section “Proposed methodology for vibration suppression in robotic milling.” The damping force of the ECD is modeled in section “Damping force modeling of the designed ECD.” Finite element analysis is carried out to analyze the magnetic flux density and damping coefficient of the ECD in section “Finite element analysis.” The tool tip dynamics without and with ECD are modeled in section “Modeling of the tool tip dynamics.” Both the tool tip FRFs and time-domain chatter stability are analyzed through simulations in section “Simulations.” Experiments are carried out to further validate the performance of the proposed method in section “Experiments.” The article is concluded in section “Conclusion.”

Proposed methodology for vibration suppression in robotic milling

The configuration of the proposed vibration suppression method for the robotic milling process is shown in Figure 1. A high-speed spindle is installed to the end effector of a Comau NJ220 robot through fixture 1. Four ECDs are clamped to the spindle using fixture 2. The two ECDs installed facing the x-direction are working in parallel to damp the spindle vibrations in the z-direction, while the other two parallel installed ECDs facing the z-direction are working together to suppress the spindle vibrations in the x-direction.

Figure 1.

Configuration of the proposed vibration suppression method for robotic milling.

The cross-section of the designed ECD which is installed facing the z-direction (Figure 1) is shown in Figure 2. Two sets of six permanent magnets (material: NdFeB N35) with alternate magnetization directions are aligned and glued to the top and the bottom sides of the middle steel plate, respectively. The middle steel plate is installed between two copper springs and will vibrate in the x-direction during the milling process. A pair of conductive copper plates is installed on the upper and lower sides of the permanent magnets, respectively, with a small air gap in between, and a steel plate is installed next to each copper plate to strengthen the magnetic flux density inside the copper plates. ²³ The upper steel plate, lower steel plate, and the middle steel plate are all made of 1020 steel, which is a type of high-permeability magnetic steel. The design parameters of the proposed ECD are given in Table 1. The weight of each ECD is about 1 kg. The design parameters are selected based on the following rules:

1. The diameter of the spindle is 100 mm. Therefore, the depth and width of each permanent magnet are selected so that the total length and width of each ECD are smaller than 80 mm for easy installation.

2. A smaller air gap gives a larger magnetic flux density inside the coppers. But if the air gap is too small, it adds the difficulty for fabrication. Therefore, an air gap between 0.5 and 1 mm is a reasonable value, considering both magnetic flux density and feasibility of fabrication. It is selected to be 0.5 mm in this study.

3. The thicknesses of the steel plates and the permanent magnets are selected to make sure that the magnetic flux densities inside them are not saturated. This is verified through finite element analysis.

4. If the thicknesses of the copper plates are too small, the generated amount of eddy current is small. However, if the thicknesses of the copper plates are very high, the ECD will become too heavy. Therefore, the thickness of each copper plate is selected to be twice the thickness of the permanent magnets. This is to allow sufficient eddy current generation, without increasing the total weight of each ECD too much.

Figure 2.

Cross-section of the designed ECD.

Table 1.

Design parameters of the proposed ECD.

Design parameters	Values (mm)
Height of permanent magnet	2.5
Width of permanent magnet	10.0
Depth of permanent magnet	50.0
Air gap	0.5
Thickness of copper plates	5.0
Thickness of upper and lower steel plates	3.0
Thickness of middle steel plate	5.0

ECD: eddy current damper.

During the robotic milling process, the spindle vibration transfers to the ECDs and causes the middle steel plate to vibrate. The relative movement between the permanent magnets and the copper plates induces eddy currents in both the upper and the lower copper plates. These eddy currents circulate in such a way that they induce a magnetic field with an opposite polarity as the applied field, causing an electromagnetic force which damps the relative velocity between the permanent magnets and the copper plates, ²⁴ therefore leading to an improved spindle damping property. Once the spindle damping is increased, the spindle vibration and the tool tip vibration will attenuate accordingly. ²⁵

Damping force modeling of the designed ECD

For clarity of coordinates, the ECD shown in Figure 2 is taken as an example in this section for the damping force modeling. In order to determine the damping force, the magnetic flux density generated by the permanent magnets need to be calculated first. It can be obtained using the equivalent current model. Since the magnetic flux density in the z-direction is the dominant one that generates the eddy current damping force for our designed ECD structure, only the magnetic flux density in the z-direction is analyzed in this section. The effect of the magnetic flux density in the y-direction is ignored due to small magnitude. Figure 3 shows the schematic of a magnet array with alternating magnetization directions. The z-direction magnetic flux density at point P (x, y, z) generated by the jth permanent magnet can be obtained as ²⁶

$$B_z^j(x,y,z)={(-1)}^j\frac{1}{4\pi }{B}_r\sum _{i=1}^2{(-1)}^i{b}_z(x-x_{\mathrm{ji}},y-y_{\mathrm{ji}},z-z_{\mathrm{ji}})$$ (1a)

$$\begin{array}{c}\multicolumn{1}{c}{b_z(x,y,z)=\mathrm{acrtan}\left(\frac{\mathrm{yz}}{x\sqrt{x^2+y^2+z^2}}\right)}\\ \multicolumn{1}{c}{+\arctan \left(\frac{\mathrm{xz}}{y\sqrt{x^2+y^2+z^2}}\right)}\end{array}$$ (1b)

where B_r is the remanence of the permanent magnet. (x_j1, y_j1, z_j1) and (x_j2, y_j2, z_j2) are the coordinates of the left-down corner and the right-up corner of the jth permanent magnet, respectively. The total magnetic flux density B_zn at point P (x, y, z) generated by the whole magnet array can be obtained by superposition as

$$B_{\mathrm{zn}}(x,y,z)=\sum _{j=1}^N{B}_z^j(x,y,z)$$ (2)

Figure 3.

Modeling of the magnetic flux density of the permanent magnets at point P (x, y, z) in the space.

Since the steel plates are installed next to the copper plates and between the two sets of magnet arrays (Figure 2), the magnetic flux density in the copper plates is strengthened. The image method is adopted to analyze the influence of the steel plates. For double-sided ferromagnetic boundaries, the infinite number of image groups need to be considered. Using coordinate transformation, the final magnetic flux density is equal to the sum of the contributions of each imaginary magnetic field and the original magnetic field which can be expressed as ²⁷

$$\begin{array}{c}\multicolumn{1}{c}{B_z(x,y,z)=B_{\mathrm{zn}}(x,y,z)+\sum _{k=1}^{\infty }}\\ \multicolumn{1}{c}{(B_{\mathrm{zn}}(x,y,z+\mathrm{kH})+B_{\mathrm{zn}}(x,y,z-\mathrm{kH}))}\end{array}$$ (3)

where H is the distance between the bottom surface of the lower steel plate and the upper surface of the top steel plate.

Figure 4 shows the electric field distribution of the upper copper plate in an alternating magnetic field generated by the permanent magnets. It is in a bottom view of the upper copper plate. The length and width of the projection area of each magnetic pole in the upper copper plate are 2a and 2b, respectively. The relative velocity between the copper plate and the permanent magnets is v, and the width of the copper plate is L. Under the magnetic pole areas in which the magnetic flux directions are pointing inside the copper plate (red color in Figure 4), the negative charges move downward and are stacked at the lower ends of the pole projections. Under the magnetic pole areas in which the magnetic flux directions are pointing outside the copper plate (blue color in Figure 4), the negative charges move upward and are stacked at the upper ends of the pole projections.

Figure 4.

Electric field distribution of the upper copper plate in an alternating magnetic field.

The magnetic flux is assumed to be uniform inside the pole areas and zero outside. Therefore, an average magnetic flux density is used here for analysis. The average magnetic flux density $B_{\mathrm{za}}^i$ for the ith magnetic pole can be calculated by integrating the magnetic flux density $B_z$ in the copper plate volume V under the ith magnetic pole as

$$B_{\mathrm{za}}^i=\underset{V}{\int }{B}_z(x,y,z)\mathrm{dV}$$ (4)

It can be seen from Figure 4 that the upper copper plate has six magnetic pole projection areas inside it. The electrostatic fields generated by Coulomb charges in the six pole areas interact with all other electrostatic fields in the copper plate. Therefore, the net electrostatic field E_y can be obtained by superposition ²⁸

$$E_y(x,y)=\sum _{i=1}^6{E}_y^i(x,y)$$ (5a)

$$\begin{array}{c}\multicolumn{1}{c}{E_y^i(x,y)=-\frac{{\mathrm{vB}}_{\mathrm{za}}^i}{2\pi }[\arctan \frac{x-b-2a(i-1)}{y-a}-\arctan \frac{x+b-2a(i-1)}{y-a}}\\ \multicolumn{1}{c}{+\arctan \frac{x+b-2a(i-1)}{y+a}-\arctan \frac{x-b-2a(i-1)}{y+a}](i=1,2,\dots ,6)}\end{array}$$ (5b)

where $B_{\mathrm{za}}^i$ is the average magnetic flux density in the z-direction in the copper plate for the ith magnetic pole. $E_y^i(x,y)$ is the electrostatic field generated by the ith magnetic pole. The current density $J_y$ of the eddy currents in the upper copper plates can be obtained as

$$J_y(x,y)=\sigma (E_y(x,y)+v{B}_{\mathrm{za}}(x,y))$$ (6)

where $\sigma$ is the electrical conductivity of the copper plate, v is the relative velocity between the permanent magnets and the copper plates, and $B_{\mathrm{za}}(x,y)$ is the average magnetic flux density of the magnetic pole area in which (x, y) is located.

The current density is calculated using equation (6) under the assumption that the copper plate has an infinite width. However, the width of the copper plate is finite in reality. This finite width can be taken into account by imposing the boundary condition of zero eddy current on the lateral edges of the copper plate (x′ = ±L/2). The method of images can be used to satisfy this boundary condition. First, the expression of the current density $J_y$ is transformed from the x–y coordinate system to the x′–y′ coordinate system (Figure 4) to make it symmetric

$$J_y(x\prime ,y\prime )=J_y(x+5b,y)$$ (7)

Then, an infinite number of imaginary eddy currents at the right and left sides of the copper plate are considered in the current density calculation ²⁶

$$\begin{array}{c}\multicolumn{1}{c}{J_{\mathrm{fy}}(x\prime ,y\prime )=J_y(x\prime ,y\prime )+\sum _{n=1}^{\infty }{(-1)}^n}\\ \multicolumn{1}{c}{[J_y(x\prime -\mathrm{nL},y\prime )+J_y(x\prime +\mathrm{nL},y\prime )]}\end{array}$$ (8)

where $J_{\mathrm{fy}}$ is the current density of the copper plate with a finite width.

Thus, the damping force $(F_u)$ acting on the upper copper plate can be calculated based on the Lorentz force law as follows

$$F_u=\underset{V_t}{\int }{J}_{\mathrm{fy}}{B}_{\mathrm{za}}\mathrm{dV}$$ (9)

where V_t is the volume of the upper copper plate under all the magnetic pole areas. Since the eddy current damping force is generated on both the upper and lower copper plates of the designed ECD, the total eddy current damping force (F) of the ECD is

$$F=2F_u$$ (10)

Finite element analysis

Magnetic field analysis

In this section, finite element analysis is carried out to analyze the magnetic flux density of the designed ECD and calculate the induced eddy current damping force. The commercial software COMSOL Multiphysics is used for the analysis.

Figure 5 shows the magnetic flux density distribution of the designed ECD in the zx plane. The color grade represents for the level of the magnetic flux density norm, while the red arrows indicate the distribution of the magnetic flux density in the x- and z-axes. Figure 6 shows the three-dimensional (3D) magnetic flux density distribution of the designed ECD. The color grade in Figure 6 represents for the z component of the magnetic flux density. In Figure 6, L1 to L6 represent the selected lines that are used for the analysis of the magnetic flux density distribution. L1 is parallel to the y-axis, while the others are all parallel to the x-axis. L1, L2, L3, and L4 are on the lower surface of the upper copper plate. L5 is on the upper surface of the upper copper plate, and L6 is in the middle of the upper steel plate. L1 is located above one of the south poles of the permanent magnets. L3, L5, and L6 share the same y coordinate.

Figure 5.

2D magnetic flux density distribution of the designed ECD in the zx plane.

Figure 6.

3D magnetic flux density distribution of the designed ECD.

The distribution of the magnetic flux density components B_z and B_y along L1 is shown in Figure 7, and it can be seen that the magnitude of B_z is small around the two ends (y = 0–2 mm and y = 48–50 mm) and is relatively large and even in the rest of the range. However, it is a different case for the magnitude of B_y, which is getting larger when it is closer to the ends (y = 0 or y = 50 mm) and getting close to zero in the middle range (y = 10–40 mm). Figure 8 shows the magnetic flux density distribution inside the upper copper plate along L2, L3, and L4. It can be seen that the directions of the magnetic flux densities B_z and B_y alternate along the x-axis, and the magnitudes of B_z and B_y are different along L2, L3, and L4. The way how these magnitudes change matches the curves shown in Figure 7. The magnetic flux density distribution along L3, L5, and L6 is shown in Figure 9. The magnitudes of B_z decrease along the z-axis in the direction far away from the permanent magnets, while the magnitudes of B_y are almost zero along L3, L5, and L6. The analysis results further verify the assumption in section “Damping force modeling of the designed ECD,” in which we assume that the effect of the magnetic flux density B_y can be ignored when modeling the eddy current damping force.

Figure 7.

Magnetic flux density distribution inside the upper copper plate along L1: (a) B_z and (b) B_y.

Figure 8.

Magnetic flux density distribution inside the upper copper plate along L2, L3, and L4: (a) B_z and (b) B_y.

Figure 9.

Magnetic flux density distribution along L3, L5, and L6: (a) B_z and (b) B_y.

Damping force analysis

The relation between the damping force (F) and the relative velocity (v) is analyzed using the finite element method with COMSOL Multiphysics as well, and the results are shown in Figure 10. It can be found that the damping force and the relative velocity have opposite signs and the magnitude of the damping force increases when the relative velocity increases. The damping coefficient (c) of the ECD is calculated from the curve shown in Figure 10, which turns out to be 82.8 N s/m.

Figure 10.

Relation between the damping force (F) and the relative velocity (v).

In this research, the damping coefficient is assumed to be a constant and is not changing with the frequency of the relative velocity. This assumption is evaluated using the magnetic Reynold number defined as ²⁶

$$R_m=\frac{l{v}_c}{\eta }$$ (11)

where l is the characteristic length scale, v_c is the characteristic velocity, and $\eta =1/({\sigma }_c{\mu }_c)$ is the magnetic diffusivity. l is selected to be half the thickness of the copper plate which is 2.5 mm. v_c is selected to be 0.1 m/s which is the maximum relative velocity. ${\sigma }_c=5.7\times {10}^7\mathrm{S}/\mathrm{m}$ and ${\mu }_c=1.257\times {10}^{-6}\mathrm{H}{\mathrm{m}}^{-1}$ . It can be calculated that $R_m=0.0179<<1$ . Therefore, the damping coefficient can be assumed to be a constant value without considering the frequency effect in this study.

Modeling of the tool tip dynamics

Tool tip dynamics modeling without ECDs

Figure 11 shows the mass–spring–damper model of the milling robot without ECD. In this model, the dynamics of the milling tool, the spindle, and the robot at the end of its arm are all assumed to be a single-degree-of-freedom (DOF) mass–spring–damper system for simplicity. They are connected in series. The x-axis is taken as an example for the model analysis in this section. F_x is the cutting force in the x-direction. x₁, x₂, and x₃ are the displacements of the milling tool tip, the spindle, and the robot end in the x-axis, respectively. m_i, k_i, and c_i (i = 1, 2, 3) are the modal mass, modal stiffness, and damping coefficient of the system, respectively. The subscripts 1, 2, and 3 refer to the milling tool, the spindle, and the robot, respectively.

Figure 11.

The mass–spring–damper model of the milling robot without ECD.

From Newton’s second law, the equations of motion for masses m₁, m₂, and m₃ can be derived as follows

$$\{\begin{array}{c}-m_1{\stackrel{··}{x}}_1-c_1({\stackrel{·}{x}}_1-{\stackrel{·}{x}}_2)-k_1(x_1-x_2)=F_x\hfill \\ \multicolumn{1}{c}{-m_2{\stackrel{··}{x}}_2-c_2({\stackrel{·}{x}}_2-{\stackrel{·}{x}}_3)-k_2(x_2-x_3)=c_1({\stackrel{·}{x}}_1-{\stackrel{·}{x}}_2)}\\ \multicolumn{1}{c}{+k_1(x_1-x_2)}\\ \multicolumn{1}{c}{-m_3{\stackrel{··}{x}}_3-c_3{\stackrel{·}{x}}_3-k_3{x}_3=c_2({\stackrel{·}{x}}_2-{\stackrel{·}{x}}_3)+k_2(x_2-x_3)}\end{array}$$ (12)

The state space model of the system can be further obtained as

$$\{\begin{array}{c}\stackrel{·}{q}=A_{q}q+B_q{F}_x\hfill \\ \multicolumn{1}{c}{x_1=C_{q}q+D_q{F}_x}\end{array}$$ (13)

where $q=[\begin{array}{cccccc}{\stackrel{·}{x}}_3\hfill & x_3\hfill & {\stackrel{·}{x}}_2\hfill & x_2\hfill & {\stackrel{·}{x}}_1\hfill & x_1\hfill \end{array}{]}^T$ is the state vector. The cutting force F_x is the system input, and the tool tip displacement x₁ is the system output. The state space matrices of the system are given as follows

$$A_q=\left[\begin{array}{cccccc}\frac{c_2-c_3}{m_3}\hfill & \frac{k_2-k_3}{m_3}\hfill & -\frac{c_2}{m_3}\hfill & -\frac{k_2}{m_3}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{c_2}{m_2}}& \frac{k_2}{m_2}\hfill & \frac{c_1-c_2}{m_2}\hfill & \frac{k_1-k_2}{m_2}\hfill & -\frac{c_1}{m_2}\hfill & -\frac{k_1}{m_2}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{c_1}{m_1}\hfill & \frac{k_1}{m_1}\hfill & -\frac{c_1}{m_1}\hfill & -\frac{k_1}{m_1}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]$$ (14a)

$$B_q={\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & -\frac{1}{m_1}\hfill & 0\hfill \end{array}\right]}^T$$ (14b)

$$C_q=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right],D_q=\left[0\right]$$ (14c)

Tool tip dynamics modeling with ECDs

When the ECDs are installed onto the milling spindle, the mass–spring–damper model of the milling robot is developed as shown in Figure 12. The two ECDs that are used to damp the vibration in the x-axis are equivalent to one mass–spring–damper model. x₄ represents the equivalent displacement of the moving part of the ECDs. m₄, k₄, and c₄ are the equivalent modal mass, modal stiffness, and damping coefficient of the ECDs, respectively.

Figure 12.

The mass–spring–damper model of the milling robot with ECDs.

From Newton’s second law, the equations of motion for masses m₁, m₂, m₃, and m₄ can be written as follows

$$\{\begin{array}{c}-m_1{\stackrel{··}{x}}_1-c_1({\stackrel{·}{x}}_1-{\stackrel{·}{x}}_2)-k_1(x_1-x_2)=F_x\hfill \\ \multicolumn{1}{c}{-m_2{\stackrel{··}{x}}_2-c_2({\stackrel{·}{x}}_2-{\stackrel{·}{x}}_3)-k_2(x_2-x_3)=c_1({\stackrel{·}{x}}_1-{\stackrel{·}{x}}_2)}\\ \multicolumn{1}{c}{+k_1(x_1-x_2)+c_4({\stackrel{·}{x}}_4-{\stackrel{·}{x}}_2)+k_4(x_4-x_2)}\\ \multicolumn{1}{c}{-m_3{\stackrel{··}{x}}_3-c_3{\stackrel{·}{x}}_3-k_3{x}_3=c_2({\stackrel{·}{x}}_2-{\stackrel{·}{x}}_3)+k_2(x_2-x_3)}\\ \multicolumn{1}{c}{-m_4{\stackrel{··}{x}}_4-c_4({\stackrel{·}{x}}_4-{\stackrel{·}{x}}_2)-k_4(x_4-x_2)=0}\end{array}$$ (15)

The state space model of the robotic milling system can be derived as

$$\{\begin{array}{c}\stackrel{·}{p}=A_{p}p+B_p{F}_x\hfill \\ \multicolumn{1}{c}{x_1=C_{p}p+D_p{F}_x}\end{array}$$ (16)

where $p=[\begin{array}{cccccccc}{\stackrel{·}{x}}_4\hfill & x_4\hfill & {\stackrel{·}{x}}_3\hfill & x_3\hfill & {\stackrel{·}{x}}_2\hfill & x_2\hfill & {\stackrel{·}{x}}_1\hfill & x_1\hfill \end{array}{]}^T$ is the state vector. The system input is the cutting force F_x and the system output is the tool tip displacement x₁. The state space matrices of the robotic milling system are given as follows

$$A_p=\left[\begin{array}{cccccccc}-\frac{c_4}{m_4}\hfill & -\frac{k_4}{m_4}\hfill & 0\hfill & 0\hfill & \frac{c_4}{m_4}\hfill & \frac{k_4}{m_4}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{c_2-c_3}{m_3}\hfill & \frac{k_2-k_3}{m_3}\hfill & -\frac{c_2}{m_3}\hfill & -\frac{k_2}{m_3}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-\frac{c_4}{m_2}}& -\frac{k_4}{m_2}\hfill & \frac{c_2}{m_2}\hfill & \frac{k_2}{m_2}\hfill & \frac{c_1+c_4-c_2}{m_2}\hfill & \frac{k_1+k_4-k_2}{m_2}\hfill & -\frac{c_1}{m_2}\hfill & -\frac{k_1}{m_2}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{c_1}{m_1}\hfill & \frac{k_1}{m_1}\hfill & -\frac{c_1}{m_1}\hfill & -\frac{k_1}{m_1}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]$$ (17a)

$$B_p={\left[\begin{array}{cccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -\frac{1}{m_1}\hfill & 0\hfill \end{array}\right]}^T$$ (17b)

$$C_p=\left[\begin{array}{cccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right],D_p=\left[0\right]$$ (17c)

Simulations

Analysis of the tool tip FRFs

The tool tip FRF in the x-axis is taken as an example for the analysis in this section. The modal parameters of the robotic milling system and the ECDs used in the simulation are given in Table 2. The FRF from the cutting force (F_x) to the tool tip displacement (x₁) is calculated for the cases with and without ECDs, as shown in Figure 13. The resonant frequency of the ECDs are tuned to 493 Hz which is the third modal frequency of the tool tip FRF caused by the milling tool itself (Figure 13). This is because the tooth passing frequency is normally high in the milling process, and the high-frequency mode of the system is easier to be excited to cause chatter.

Table 2.

Modal parameters of the robotic milling system and ECDs.

Parameters	Values
m 1	0.2 kg
c 1	15 N s/m
k 1	2 × 106 N/m
m 2	6 kg
c 2	180 N s/m
k 2	1 × 107 N/m
m 3	10 kg
c 3	30 N s/m
k 3	1 × 105 N/m
m 4	0.6 kg
c 4	165.6 N s/m
k 4	5.76 × 106 N/m

ECD: eddy current damper.

Figure 13.

FRF from the cutting force (F_x) to the tool tip displacement (x₁).

From Figure 13, it can be seen that when no ECD is attached to the spindle, the tool tip FRF has three dominant modal frequencies at 25, 135, and 493 Hz, respectively. The peak amplitudes of the FRFs in Figure 13 are given in Table 3. With ECDs installed, the second and the third peak of the tool tip FRF are reduced by 3.4 and 5.2 dB, respectively, while the first peak is slightly increased by 0.9 dB.

Table 3.

Peak amplitudes of the FRFs.

Cases	Peak amplitudes (dB)
	Mode 1	Mode 2	Mode 3
Without ECD	−73.6	−85.0	−92.5
With ECDs	−72.7	−88.4	−97.7

FRF: frequency response function; ECD: eddy current damper.

Since the robot stiffness will change with its postures, three different values are selected for the stiffness (k₃) at the end of the robot to test the robustness of the designed ECDs. The modal parameters (m₄, c₄, k₄) of the ECDs are still the same as the ones given in Table 2. The resulted tool tip FRFs are shown in Figure 14. For the cases without ECD, the first modal frequency changes to 13 Hz when the stiffness k₃ reduces from 1.0 × 10⁵ N/m (case 2) to 3.0 × 10⁴ N/m (case 1), and it increases to 39 Hz when the stiffness k₃ rises to 2.4 × 10⁵ N/m (case 3). With ECDs installed, the damping performance is quite consistent even when the stiffness (k₃) of the robot has changed significantly, which shows good robustness of the designed ECDs.

Figure 14.

FRF from the cutting force (F_x) to the tool tip displacement (x₁) when robot stiffness changes: case 1—k₃ = 3.0 × 10⁴ N/m; case 2—k₃ = 1.0 × 10⁵ N/m; and case 3—k₃ = 2.4 × 10⁵ N/m.

Chatter stability analysis

The milling tool is considered to have two orthogonal DOFs in the x- and z-axes with uncoupled dynamics and is assumed to have N number of teeth with a zero helix angle. The cutting forces excite the structure in the feed (x) and normal (z) directions, causing dynamic displacements x₁ and z₁ at the tool tip, respectively. The block diagram of the robotic milling process is shown in Figure 15. According to Altintas, ⁶ the total chip thickness can be expressed as

$$h({\varphi }_j)=h_s({\varphi }_j)+h_d({\varphi }_j)$$ (18a)

$$\{\begin{array}{c}{h}_s({\varphi }_j)=s_t\sin ({\varphi }_j)g({\varphi }_j)\hfill \\ \multicolumn{1}{c}{h_d({\varphi }_j)=[\mathrm{\Delta }{x}_1\sin ({\varphi }_j)+\mathrm{\Delta }{z}_1\cos ({\varphi }_j)]g({\varphi }_j)}\end{array}$$ (18b)

$$g({\varphi }_j)=\{\begin{array}{c}\begin{array}{cc}1,\hfill & {\varphi }_{\mathrm{st}}<{\varphi }_j<{\varphi }_{\mathrm{ex}}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}0,\hfill & {\varphi }_j<{\varphi }_{\mathrm{st}}\mathrm{or}{\varphi }_j>{\varphi }_{\mathrm{ex}}\hfill \end{array}}\end{array}$$ (18c)

where h_s is the static chip thickness; h_d is the dynamic chip thickness; s_t is the feed rate per tooth; Φ _j is the instantaneous angular immersion of the tooth j measured clockwise from the normal (z) axis; Φ _st and Φ _ex are the start and exit immersion angles of the cutter to and from the cut, respectively. $\mathrm{\Delta }{x}_1=x_1-x_{1p}$ and $\mathrm{\Delta }{z}_1=z_1-z_{1p}$ . ( $x_1$ , $z_1$ ) and ( $x_{1p}$ , $z_{1p}$ ) are the dynamic displacements of the tool tip at the present and previous tooth periods, respectively.

Figure 15.

Block diagram of the robotic milling process.

The tangential (F_tj) and radial (F_rj) cutting forces acting on the tooth j can be obtained as

$$F_{\mathrm{tj}}=K_t\mathrm{ah}({\varphi }_j),F_{\mathrm{rj}}=K_r{F}_{\mathrm{tj}}$$ (19)

where K_t and K_r are the cutting coefficients and a is the axial depth of cut. The cutting forces in the x- and z-axes can be expressed as

$$\{\begin{array}{c}{F}_{\mathrm{xj}}=-F_{\mathrm{tj}}\cos ({\varphi }_j)-F_{\mathrm{rj}}\sin ({\varphi }_j)\hfill \\ \multicolumn{1}{c}{F_{\mathrm{zj}}=F_{\mathrm{tj}}\sin ({\varphi }_j)-F_{\mathrm{rj}}\cos ({\varphi }_j)}\end{array}$$ (20)

By summing the cutting forces contributed by all teeth, the total milling forces acting on the cutter can be obtained as

$$F_x=\sum _{j=0}^{N-1}{F}_{\mathrm{xj}}({\varphi }_j),F_z=\sum _{j=0}^{N-1}{F}_{\mathrm{zj}}({\varphi }_j)$$ (21)

A slotting process is taken as an example for the cutting force simulation in this section. The parameter values used for calculating the cutting forces are given in Table 4. The cutting coefficients (K_t, K_r) of aluminum 7075 ⁶ are used in the simulation. Assuming that the robotic milling system has the same dynamics in the x- and z-axes, the modal parameters given in Table 2 are used in this section for modeling the dynamics of the robotic milling system in both axes. Since the changing stiffness (k₃) of the robot has little effect on the third peak of the tool tip FRF (Figure 14) which is easier to cause chatter in the high-speed milling process, a constant stiffness (k₃) is used in the chatter stability analysis. The chatter stability lobes are calculated for the cases without and with ECDs using the zero-order solution method. ⁶ The tool tip dynamics in both the x- and z-axes are considered when calculating the stability lobe diagrams. The results are shown in Figure 16. It can be found that the minimum critical stable depth of cut (a_lim) is only about 0.05 mm for the high spindle speed range above 20,000 r/min without ECD. However, when the ECDs are installed, this minimum critical stable depth of cut (a_lim) is improved twofold to 0.1 mm.

Table 4.

Values of parameters used in the cutting force calculation.

Parameters	Values
Tangential cutting force coefficient, Kt	796 MPa
Cutting force coefficient, Kr	0.212
Tooth number, N	2
Feed rate per tooth, st	0.01 mm
Start immersion angle, Φ st	0
Exit immersion angle, Φ ex	π rad

Figure 16.

Chatter stability lobe diagrams for the cases without and with ECDs.

The cutting forces and tool tip displacements at different spindle speeds and axial depths of cut are also simulated for the cases without and with ECDs. The saturation tool tip displacements are set to ±1 mm for both x₁ and z₁ in the simulation. Since the information of the cutting forces and tool tip displacements in one axis is enough to tell whether chatter happens or not, only the results in the x-axis are plotted. The results are shown in Figures 17 and 18. At 27,000 r/min spindle speed, the tool tip chatters at 0.06 mm axial depth of cut without ECD. Chatter is caused by the milling tool mode at 493 Hz. With ECDs being installed, the cutting process is stable at 0.095 mm axial depth of cut with no chatter happening. The small peak at 900 Hz in the fast Fourier analysis of x₁ is due to forced vibration caused by the cutting force F_x at the tooth passing frequency. At 40,000 r/min spindle speed, chatter happens at 0.07 mm axial depth of cut without ECD. It is also caused by the milling tool mode at 493 Hz. When ECDs are attached to the spindle, the cutting process is stable at 0.12 mm axial depth of cut, which is about 71.4% improvement. The small peak at 1333 Hz in the fast Fourier analysis of x₁ is also because of the forced vibration caused by the cutting force F_x at the tooth passing frequency.

Figure 17.

Cutting forces and tool tip displacements during the robotic milling process at the spindle speed of 27,000 r/min.

Figure 18.

Cutting forces and tool tip displacements during the robotic milling process at the spindle speed of 40,000 r/min.

Experiments

The experimental setup for vibration suppression in the robotic milling process is shown in Figure 19. A high-speed Jager spindle is installed to the end of a Comau NJ220 robot using an ATI QC160 tool changer. Two pairs of ECDs are clamped to the spindle with fixtures. The pair of ECDs installed facing the x-direction is working in parallel to damp the spindle vibrations in the z-direction. The two ECDs installed facing the z-direction are working together to suppress the x-direction spindle vibrations.

Figure 19.

Experimental setup.

The direct tool tip FRFs without ECD are first measured by giving a hammer blow at the tool tip and collecting the resulted tool tip acceleration signals in the same direction. The measured tool tip FRFs in both x- and z-axes are shown in Figure 20. The system has five dominant modes. The peaks around 690 and 992 Hz are caused by the modes of the milling tool. The peaks in the low-frequency region around 28 and 52 Hz are caused by the modes of the robot. The peak near 120 Hz is caused by the spindle mode. Since the milling tool mode around 690 Hz is the most flexible one and most likely to be excited during robotic milling, the frequencies of the ECDs are tuned to 690 Hz. When the ECDs are installed onto the spindle, the tool tip FRFs are measured again and the results are shown in Figure 20. It can be found that the peaks around 690 and 992 Hz are damped by 15.5% and 8.5%, respectively, in the x-axis, while they are damped by 22.1% and 12.4%, respectively, in the z-axis. The robot modes are not damped much. The spindle modes are slightly damped by 9.4% and 5.3% in the x- and z-axes, respectively.

Figure 20.

Measured tool tip FRFs without and with ECDs: (a) FRF from F_x to x₁ and (b) FRF from F_z to z₁.

The identified modal parameters for the measured tool tip FRFs without and with ECDs are given in Tables 5 and 6. It can be seen that the damping ratio of the spindle modes and the milling tool modes are increased when ECDs are installed. The chatter stability lobes are also calculated based on the identified modal parameters (Tables 5 and 6) and the process parameters (Table 4) using zero-order solution. The tool tip dynamics in both the x- and z-axes are considered when calculating the stability lobe diagrams. The results are shown in Figure 21. It can be found that the minimum critical stable depth of cut (a_lim) is improved from 0.085 to 0.105 mm in the high spindle–speed zone, which is about 23.5% increment.

Table 5.

Identified modal parameters for tool tip FRFs without ECD.

Modes	FRF from Fx to x1		FRF from Fz to z1
	Frequency (Hz)	Damping ratio (%)	Frequency (Hz)	Damping ratio (%)
Mode 1	27.5	0.80	28	0.95
Mode 2	52	2.14	51	2.75
Mode 3	123	5.23	120.5	5.37
Mode 4	693	2.55	686	3.18
Mode 5	993	1.67	992	1.94

FRF: frequency response function; ECD: eddy current damper.

Table 6.

Identified modal parameters for tool tip FRFs with ECDs.

Modes	FRF from Fx to x1		FRF from Fz to z1
	Frequency (Hz)	Damping ratio (%)	Frequency (Hz)	Damping ratio (%)
Mode 1	27	1.35	26.8	0.87
Mode 2	52	2.27	51.3	2.63
Mode 3	119	6.03	118	5.81
Mode 4	687.5	3.07	683	4.19
Mode 5	991	1.83	989	2.31

FRF: frequency response function; ECD: eddy current damper.

Figure 21.

Chatter stability lobe diagrams for the cases without and with ECDs, calculated with the identified modal parameters.

Robotic milling tests have been carried out to further verify the effectiveness of the designed ECDs. The cutting parameter configurations chosen in the cutting tests are labeled with a black dot in Figure 21. The spindle speed is selected to be 10,650 r/min in all the cutting tests. For the case without ECD, we first cut at 0.2 mm depth of cut (a). The process is stable and the workpiece surface is good. When we increase the depth of cut (a) to 0.3 mm, the process becomes unstable and the surface finish is poor, as shown in Figure 22. For the case with ECDs, we cut with a 0.4 mm depth of cut (a) first, and the cutting process is stable with a decent surface finish. However, when we raise the depth of cut (a) to 0.5 mm, chatter occurs and the surface finish is quite poor, as shown in Figure 22. From the cutting tests, we can see that the stable depth of cut is improved twofold at the selected spindle speed with ECDs.

Figure 22.

Robotic milling results without and with ECDs at the spindle speed of 10,650 r/min.

Conclusion

This article presents a novel ECD design for chatter suppression in the robotic milling process. The proposed ECDs are installed on the milling spindle to damp the tool tip vibrations. The structural design and the working principle of the ECD are introduced, and the finite element method is used to analyze the magnetic flux density distribution and magnetic force generation of the designed ECD. The tool tip dynamics without and with ECDs are modeled and simulated, and the peaks of the tool tip FRF caused by the spindle and milling tool modes are damped by 3.4 and 5.2 dB, respectively, in the simulation. Experiments are also carried out to further verify the damping performance of the designed ECDs. The results show that the peaks of the tool tip FRFs caused by the milling tool modes are damped by 15.5% and 8.5%, respectively, in the x-axis, while they are damped by 22.1% and 12.4%, respectively, in the z-axis. The minimum critical stable depth of cut is increased from 0.085 to 0.105 mm in the high spindle–speed region, which is about 23.5% improvement. Also, the stable depth of cut is improved twofold at 10,650 r/min, which is verified by cutting tests.