Research on machining parameter optimization in finishing milling with multiple constraints

Abstract

High efficiency and high precision milling, as the eternal goal of CNC machining, needs to balance many constraints for selecting the most reasonable processing parameters. This paper presents an efficient machining parameter optimization method for finishing milling operation with multiple constraints. Firstly, under the multiple constraints of parameter feasible region, milling force, milling stability, roughness, and machining contour accuracy, a multi-variable parameter optimization model with machining efficiency as the objective is established. A four level cycle optimization strategy has been detailly described for solving the optimization problem, in which the feed per tooth is optimized by using the golden section method, and with the aid of the random vector search method, the spindle speed, radial, and axial depth cuts are both numerically iterated. The optimal machining parameter combination of the tooth number, feed per tooth, spindle speed, radial, and axial depth of cuts are achieved at last. Finally, the experimental verification results show that the proposed method can greatly improve the machining efficiency under chatter free condition and achieve an efficient finishing milling with consideration of the multiple constraints.

Introduction

The goal of CNC machining is to achieve high efficiency and high precision, that is, to maximize the material removal rate (MRR) under the premise of ensuring machining quality and other constraints. In the traditional methods, operators mostly do repeated experiments and use manual experience as the basis reference to select the machining parameters. Unfortunately, this approach is always time-consuming and uneconomical, and the selected parameters are mostly conservative, thus making it difficult to achieve a real high efficiency machining. ¹ There is no doubt that increasing machining parameters is the most effective way to improve the processing efficiency. However, it is easily prone to generate milling vibrations and chatter, ² even damage the machine and tools, and deteriorate the machining quality of parts. Therefore, an effective method is in urgent need of research to guide the selection of reasonable machining parameters for high efficiency and high precision.

In order to improve the MRR, many scholars have studied the parameter optimization problem taking surface roughness into account. Onwubolu ³ presented a new methodology to use the roughness as the one of bounds for the autonomous selection of cutting conditions in milling operations by using Tribes. Nair et al. ⁴ took surface roughness and MRR as the optimization targets, and systematically analyzed the experimental data. They found that when the roughness evaluation characteristics are selected differently, the results obtained by the four different optimization methods are not the same, which indicated that the optimization need to consider the influence of the parameter levels. Zerti et al. ⁵ analyzed the influence of different machining parameters on roughness, cutting force, and MRR in dry turning process, then established four optimization goals, such as quality, efficiency, energy consumption, and cutting force, and further obtained the optimal parameters by using neural networks. It is expected that the roughness is mainly affected by the feed speed. Qu et al. ⁶ studied the parameter optimization problem in thin-walled part milling process, and established the multi-objective optimization problem with three contradictory goals, including roughness, milling force, and MRR, lastly obtained the optimal parameter combination of spindle speed, feed per tooth, and axial depth of cut. Tebassi et al. ⁷ experimentally investigated the influence of cutting parameters on roughness, cutting force, machining efficiency, and energy consumption during the turning process, and proposed a multi-objective optimization method to obtain the smallest roughness, cutting force and machining efficiency, and maximum MRR. Dikshit et al. ⁸ carried out a mathematical model of surface roughness for milling process in terms of feed per tooth, cutting speed, axial, and radial depth of cuts. Here, a second-order full quadratic model was developed and the optimization using the teaching-learning-based methodology could take minimum effort and provide better results.

At the same time, obtaining high MRR has also been deeply studied by scholars in the case of ensuring the contour precision of machined surface, which is mainly affected by the vibration of the tool and workpiece during cutting process. Budak and Tekeli ⁹ studied the optimal selection method for the combination of radial and axial depth of cuts under stable milling condition. In the given analysis case, the milling time was reduced about 40%. Kurdi ¹⁰ carried out a multi-objective optimization problem with milling stability and surface location error as the objective function. When calculating the surface location error, a certain disturbance amount of the spindle speed was introduced to avoid the rapid change on the surface location error caused by the resonance region, and the robustness of the optimization model was enhanced. Based on the analysis of milling stability and surface location error, Insperger et al. ¹¹ discussed a method to optimize the spindle speed when the radial and axial depth of cuts were determined. The study pointed out that the optimal spindle speed was selected from the left side of the resonance zone, which provided an enlightening idea for the optimization of machining parameters. Based on milling dynamics, Merdol and Altintas ¹² divided the milling parameter optimization into two steps. Firstly, the spindle speed, radial, and axial depth of cuts were selected with consideration of milling stability, torque, and power limits. Secondly, the milling force, milling stability, torque, and power were checked along the whole path after generating the NC program and the feed speed and spindle speed were adjusted to avoid overruns. Wan et al. ¹³ proposed a parameter optimization method to ensure the surface dimensional accuracy based on statics analysis, in which the linear relationship between surface error and feed per tooth and radial depth of cut was established. Zhang and Ding ¹⁴ put forward a chatter free milling parameter optimization method, in which the tool vibration displacement was minimized and the MRR was maximized by optimizing the spindle speed and axial depth of cut. Hu et al. ¹⁵ constructed an objective function based on MRR and tool life, and proposed a milling parameter optimization model under the constraint of cutting stability. The expected value of MRR and tool life were used in the optimization model to realize the quantitative adjustment of the optimization parameters and the effective control of the optimization direction. Li et al. ¹⁶ optimized the milling parameters based on the dynamic displacement response for a ruled surface impeller machining. In the method, the machining parameters and maximum dynamic displacement of the cutter were used as training samples to train the BP neural network, and the maximum dynamic displacement of the cutter was taken as the optimization objective. Mokhtari et al. ¹⁷ improved the micro-milling tool geometry to find the maximum allowable cutting depth for the chatter-free machining, in which the objective function was defined as the maximization of the minimal cutting depth within the stability lobe diagram.

For the thin-walled part milling, Alan et al. ¹⁸ studied the variation of dynamic stiffness and milling stability of the parts due to material removal, and selected the axial depth of cut under its critical limits for improving the machining efficiency. Budak et al. ¹⁹ built up a prediction method of time-varying dynamics as the margin was removed layer by layer after in five-axis milling of blade parts, then established a milling stability domain diagram that changes with the cutting position, lastly by which a more reasonable spindle speed and axial depth of cut were selected. Germashev et al. ²⁰ analyzed the vibration response of the flexible thin-walled parts during the milling process and found that the dynamic action between the force and workpiece would be different at various spindle speeds, and a better spindle speed was selected by considering the vibration amplitude of workpiece. Recently, Ringgaard et al. ²¹ presented an optimization method to maximize the MRR without violating forced vibration and chatter stability constraints. The cost function was formulated using penalty terms which penalized the function as chatter and forced vibration constraints. At last, the optimal spindle speed, feed per tooth, and axial depth of cut were achieved in milling of thin-walled structures.

The above literature reviews show that although lots of research have been focused on the machining parameter optimization by many scholars, there are still some shortcomings in these existing methods. Firstly, some research does not consider the physical dynamics properties of the machining system, especially the dynamical engagement characteristics between the cutter and workpiece, including milling stability and machining accuracy. Secondly, only part of machining parameters, such as spindle speed and axial depth of cut, can be optimized, and other important machining parameters are ignored, which thus limits the applicability of these method. In machining process, the typical machining parameters have spindle speed, radial depth of cut, axial depth of cut, feed per tooth, and cutter tooth number, etc. Meanwhile, the constraints are multiple, such as milling stability, profile accuracy, and roughness of machined surface, etc. Hence, how to obtain the maximum machining efficiency under so many complex constraints is a difficult challenge. The purpose of this paper is to propose a reasonable and effective milling parameter optimization method under chatter free condition by establishing and solving a single objective, multiple constraints, and multiple variables problem. The structure of this paper is organized as follows: in Section 2, considering the physical constraints of the machining system, the mathematical model of the optimization problem is established; in Section 3, the algorithm for solving the optimization model is given; in Section 4, the experimental verification is carried out, and finally the conclusion of this paper is summarized.

Mathematical model of the optimization problem

Optimization objective

Figure 1 shows a typical finishing milling process. The thickness of the finishing allowance is usually small in the normal direction of the part surface. According to the distribution of the allowance, two kinds of tool path topology are usually used to remove it: depth-first and width-first. The goal of parameter optimization is to maximize the machining efficiency under the relevant constraints. Generally, in the condition of constant cutting parameters, the milling machining efficiency can be expressed by the MRR as following.

$$f_{\mathrm{MRR}}=N_t·f_t·n·a_e·a_p/1000(\mathrm{c}{\mathrm{m}}^3/min)$$

(1)

Where N_t is the cutter tooth number. f_t is the feed per tooth (mm/tooth). n is the spindle speed (rpm). a_e is the radial depth of cut (mm). a_p is the axial depth of cut (mm).

Figure 1.

Schematic diagram of finishing milling process.

Optimization variables

It can be seen from the above formula that the optimization objective function includes five variable parameters, such as tooth number and other four basic cutting parameters (spindle speed, radial depth of cut, axial depth of cut, and feed per tooth). Thus, in this paper, the optimization variable is defined as $\mathit{x}=[n,a_e,a_p,N_t,f_t{]}^T$ .

Multiple constraints

The selection of milling parameters would be restricted by many factors, such as the feasible region of machining parameters, cutting force, milling stability, workpiece contour accuracy, and roughness of machined surface. The quantitative model of these constraints will be established step by step in the following.

Parameter feasible region constraint

During the actual machining process, the determination of optimization variable will consider multiple constraints, in which the first one is the limitation of the optional range of parameters. Considering the machining ability of spindle and cutter, there exist some upper or lower limits for the selection of spindle speed, radial and axial depth of cuts, and feed per tooth. For instance, the spindle speed should be less than or equal to its maximum allowable value, the radial depth of cut should not exceed the cutter diameter, the axial depth of cut should be less than the length of cutting edge of cutter, the selection of feed per tooth needs to take the cutting ability of cutter into account, and the choose of tooth number should ensure the chip removal ability of cutter. In summary, the feasible region of parameter selection can be formed, which is also the most basic constraint for the optimization variable, as shown in the following formula.

$$\{\begin{array}{c}n\in [n^l,n^u]\hfill \\ a_e\in [a_e^l,a_e^u]\hfill \\ a_p\in [a_p^l,a_p^u]\hfill \\ f_t\in [f_t^l,f_t^u]\hfill \\ N_t\in [N_t^l,N_t^u](N_t\in Z^{+})\hfill \end{array}$$

(2)

Where the superscript ^u denotes the upper limit. The Superscript ^l denotes the lower limit.

In addition, when choosing the depth-first tool path topology, the line spacing in width direction should be uniform. While the line spacing in axial direction should be equispaced under the width-first milling situation. Hence, one can further obtain:

$$\{\begin{array}{c}\multicolumn{1}{c}{\mathrm{Depth}-\mathrm{first}:a_e=a_e^l=a_e^u=\frac{L_w}{N_w}(N_w\in Z^{+})}\\ \multicolumn{1}{c}{\mathrm{Width}-\mathrm{first}:a_p=a_p^l=a_p^u=\frac{L_h}{N_h}(N_h\in Z^{+})}\end{array}$$

(3)

Where $L_w$ and $L_h$ are the length of allowance in width and axial direction, respectively. $N_w$ and $N_h$ are the discrete quantity of line spacing in width and axial direction, respectively.

Milling force constraint

During machining process, milling force will be produced when the cutter and workpiece are engaging dynamically with each other. Excessive milling force will induce cutter deflection, cause impact, and even tool breakage. Therefore, it is necessary to restrict the milling force. As shown in Figure 2, according to the basic cutting force mechanism model, the tangential, radial, and axial cutting forces of the cutting point P of the ith tooth on the jth cutting unit at an arbitrarily rotation angle ${\varphi }_{i,j}(t)$ can be expressed as the sum of shearing and ploughing force components by the following formula.

$$\{\begin{array}{c}{\mathrm{dF}}_{t,i,j}(t)=(k_{\mathrm{ts}}{h}_{i,j}(t)\mathrm{db}+k_{\mathrm{tp}}\mathrm{db})W_{i,j}(t)\hfill \\ \multicolumn{1}{c}{{\mathrm{dF}}_{r,i,j}(t)=(k_{\mathrm{rs}}{h}_{i,j}(t)\mathrm{db}+k_{\mathrm{rp}}\mathrm{db})W_{i,j}(t)}\\ \multicolumn{1}{c}{{\mathrm{dF}}_{a,i,j}(t)=(k_{\mathrm{as}}{h}_{i,j}(t)\mathrm{db}+k_{\mathrm{ap}}\mathrm{db})W_{i,j}(t)}\end{array}$$

(4)

Where $k_{\mathrm{qs}}$ and $k_{\mathrm{qp}}(q=t,r,a)$ are the shearing and ploughing force coefficients in tangential, radial, and axial direction, respectively. $\mathrm{db}$ is the axial thickness of cutting unit. $W_{i,j}(t)$ is the window function for determining the cutting point whether to participate in cutting, which is given by

$$W_{i,j}(t)=\{\begin{array}{cc}1\hfill & {\theta }_{s,i,j}\le {\varphi }_{i,j}(t)\le {\theta }_{e,i,j}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{others}\hfill \end{array}$$

(5)

Where ${\theta }_{s,i,j}$ and ${\theta }_{e,i,j}$ are the start and end angle of the cutting point of the ith tooth on the jth cutting unit, respectively.

Figure 2.

Milling mechanics and dynamics model.

In equation (4), the instantaneous uncut chip thickness $h_{i,j}(t)$ of cutting point P can be decided by its two dimensional trajectory and workpiece geometry, which is equal to the line segment length $\overline{\mathrm{PT}}$ . As shown in Figure 2, point T is left on the machined surface after cutting by the previous m_ith cutter teeth at a certain delay time ${\tau }_{i,j}(t,m_i)$ in the past. At time $t-{\tau }_{i,j}(t,m_i)$ , the position of point T in the global reference coordinate system can be expressed as

$$\{\begin{array}{c}{T}_{x,i-m_i,j}=o_{\mathrm{rx}}(t-{\tau }_{i,j}(t,m_i))+R_{i-m_i,j}\sin ({\varphi }_{i-m_i}(t-{\tau }_{i,j}(t,m_i))\hfill \\ \multicolumn{1}{c}{T_{y,i-m_i,j}=o_{\mathrm{ry}}(t-{\tau }_{i,j}(t,m_i))+R_{i-m_i,j}\cos ({\varphi }_{i-m_i}(t-{\tau }_{i,j}(t,m_i))}\end{array}$$

(6)

Where $R_{i,j}$ is the actual cutting radius of the ith tooth on the jth cutting unit.

At the same time, the point T can also be derived by

$$\{\begin{array}{c}{T}_{x,i-m_i,j}=o_{\mathrm{rx}}(t)+(R_{i,j}-h_{i,j}(t,m_i))\sin ({\varphi }_{i,j}(t))\hfill \\ \multicolumn{1}{c}{T_{y,i-m_i,j}=o_{\mathrm{ry}}(t)+(R_{i,j}-h_{i,j}(t,m_i))\cos ({\varphi }_{i,j}(t))}\end{array}$$

(7)

Based on the cutter feed motion relationship, where the dynamic displacement response of the cutter is ignored, one can achieve

$$\{\begin{array}{c}{o}_{\mathrm{rx}}(t)-o_{\mathrm{rx}}(t-{\tau }_{i,j}(t,m_i))=f_{\mathrm{vx}}{\tau }_{i,j}(t,m_i)\hfill \\ \multicolumn{1}{c}{o_{\mathrm{ry}}(t)-o_{\mathrm{ry}}(t-{\tau }_{i,j}(t,m_i))=0}\end{array}$$

(8)

Where $f_{\mathrm{vx}}$ is the cutter feed speed.

Combining the above three equations, the delay time term can be derived as

$${\tau }_{i,j}(t,m_i)=\frac{R_{i-m_i,j}{\chi }_{i-m_i,i,j}}{R_{i-m_i,j}\frac{\pi n}{30}+f_{\mathrm{vx}}\cos ({\varphi }_{i,j}(t))}$$

(9)

Where ${\chi }_{i-m_i,i,j}$ is the tooth spacing angle between the (i−m_i)th and ith tooth of the jth cutting unit.

Furthermore, the available instantaneous uncut chip thickness can be obtained by

$$\begin{array}{c}\multicolumn{1}{c}{h_{i,j}(t,m_i)=f_{\mathrm{vx}}{\tau }_{i,j}(t,m_i)\sin ({\varphi }_{i,j}(t))}\\ +R_{i,j}-R_{i-m_i,j}\cos ({\chi }_{i-m_i,i,j}\hfill \\ -\frac{\pi n}{30}{\tau }_{i,j}(t,m_i))\hfill \end{array}$$

(10)

Then, the actual instantaneous uncut chip thickness is the minimum of all available chip thickness greater than zero, which can be mathematically expressed as

$$h_{i,j}(t)=max(0,min(h_{i,j}(t,m_i)))m_i=1,2,\cdots N_t$$

(11)

Finally, substituting equation (11) into (4), the cutting forces generated by all the cutting units participating in the cutting at the same time are summed in the range of effective cutting depth, and the total cutting force acting on the cutter is

$$\left[\begin{array}{c}{F}_{cx}(t)\\ F_{cy}(t)\\ F_{cz}(t)\end{array}\right]=\sum _i\sum _j\left[\begin{array}{ccc}-\cos ({\varphi }_{i,j}(t))& -\sin ({\varphi }_{i,j}(t))& 0\\ \sin ({\varphi }_{i,j}(t))& -\cos ({\varphi }_{i,j}(t))& 0\\ 0& 0& -1\end{array}\right]\left[\begin{array}{c}d{F}_{t,i,j}(t)\\ d{F}_{r,i,j}(t)\\ d{F}_{a,i,j}(t)\end{array}\right]$$

(12)

In this case, the resultant cutting force acting on the cutter is

$$F_c(\mathit{x},t)=\sqrt{F_{\mathrm{cx}}^2(t)+F_{c\mathrm{y}}^2(t)+F_{c\mathrm{z}}^2(t)}$$

(13)

Considering the impact resistance of the cutter, the maximum milling force should be limited within a certain range as following.

$$max(F_c(\mathit{x},t))\le F_{\mathrm{clim}}$$

(14)

Where $F_{\mathrm{clim}}$ is the allowable maximum milling force.

Milling stability constraint

In order to ensure the milling process in a stable state, it is necessary to consider the constraint of milling chatter, which often restricts the upper limit of milling parameter selection.

As shown in Figure 4, milling chatter is mainly caused by the regeneration mechanism of the flexible machining system. Considering the regenerative effect induced by the flexible cutter system in milling process, the transient chip thickness can be obtained as follows.

$$h_{i,j}^{~}(t)=(h_{i,j}(t)+[\sin ({\varphi }_{i,j}(t))\cos ({\varphi }_{i,j}(t))]\left[\begin{array}{c}{O}_{\mathrm{rx}}(t)-O_{\mathrm{rx}}(t-{\tau }_{i,j}(t))\hfill \\ \multicolumn{1}{c}{O_{\mathrm{ry}}(t)-O_{\mathrm{ry}}(t-{\tau }_{i,j}(t))}\end{array}\right])W_{i,j}(t)$$

(15)

Next, the milling dynamics model can be established by

$$\mathbf{M}\stackrel{··}{X}(t)+\mathbf{C}\stackrel{·}{X}(t)+\mathbf{K}X(t)=F_0(t)+K_c(t)[X(t)-X(t-\tau )]$$

(16)

Where $\mathbf{M}=\left[\begin{array}{cc}{m}_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m_y\hfill \end{array}\right],\mathbf{C}=\left[\begin{array}{cc}{c}_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& c_y\hfill \end{array}\right],\mathbf{K}=\left[\begin{array}{cc}{k}_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& k_y\hfill \end{array}\right]$ , in which $m_x,m_y,c_x,c_y$ and $k_x,k_y$ are the modal masses, modal dampings, and modal stiffnesses along the X and Y direction, respectively. $X(t)=[O_{\mathrm{rx}}(t),O_{\mathrm{ry}}(t){]}^T$ is the cutter vibration displacement vector. $F_0(t)=[F_{\mathrm{cx}}(t),F_{\mathrm{cy}}(t){]}^T$ is the nominal milling force. $\tau ={\overline{\tau }}_{i,j}(t)$ is the average delay time. The second term on the right of the above formula is the regenerative term of milling force, in which $K_c(t)$ is the coefficient matrix that can be derived as

$$\begin{array}{c}\multicolumn{1}{c}{K_c(t)=\mathrm{db}\sum _i\sum _j\left[\begin{array}{cc}-k_{\mathrm{ts}}\sin ({\varphi }_{i,j})\cos ({\varphi }_{i,j})-k_{\mathrm{rs}}\sin ({\varphi }_{i,j}{)}^2\hfill & -k_{\mathrm{ts}}\cos ({\varphi }_{i,j}{)}^2-k_{\mathrm{rs}}\sin ({\varphi }_{i,j})\cos ({\varphi }_{i,j})\hfill \\ \multicolumn{1}{c}{k_{\mathrm{ts}}\sin ({\varphi }_{i,j}{)}^2-k_{\mathrm{rs}}\sin ({\varphi }_{i,j})\cos ({\varphi }_{i,j})}& k_{\mathrm{ts}}\sin ({\varphi }_{i,j})\cos ({\varphi }_{i,j})-k_{\mathrm{rs}}\cos ({\varphi }_{i,j}{)}^2\hfill \end{array}\right]W_{i,j}(t)}\end{array}$$

(17)

Then, the system dynamics equation can be further written as following state space model.

$$\stackrel{·}{U}(t)=\mathbf{A}(t)U(t)+\mathbf{B}(t)U(t-\tau )+F(t)$$

(18)

Where $\mathbf{A}(t)=\left[\begin{array}{cc}{\mathbf{0}}_{2\times 2}& {\mathbf{I}}_{2\times 2}\\ {\mathbf{M}}^{-1}(K_c(t)-\mathbf{K})& -{\mathbf{M}}^{-1}\mathbf{C}\end{array}\right],\mathbf{B}(t)=\left[\begin{array}{cc}{\mathbf{0}}_{2\times 2}& {\mathbf{0}}_{2\times 2}\\ -{\mathbf{M}}^{-1}{K}_c(t)& {\mathbf{0}}_{2\times 2}\end{array}\right],F(t)={\left[{\mathbf{0}}_{2\times 1}{\mathbf{M}}^{-1}{F}_0\right]}^T,U={\left[X\stackrel{.}{X}\right]}^T$

In this paper, the idea of the SDM ²² is used to solve the above equation for stability prediction. Now, the cutter rotation period is needed to be discretized into a time series with a small equal interval of $\mathrm{\Delta }t=T/N_m$ . After the discretization of time and delay terms, the stability of the system would be evaluated from a nonlinear problem as a linear problem. In one tiny period of time, the solution of the above equation can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{U_{k+1}={\mathbf{P}}_k{U}_k+0.5{\mathbf{R}}_k{U}_{k-{\tau }_k+1}}\\ +0.5{\mathbf{R}}_k{U}_{k-{\tau }_k}+{\mathbf{Q}}_k{F}_k\hfill \end{array}$$

(19)

Where ${\mathbf{P}}_k=e^{{\mathbf{A}}_k\mathrm{\Delta }t},{\mathbf{R}}_k=(e^{{\mathbf{A}}_k\mathrm{\Delta }t}-\mathbf{I}){\mathbf{A}}_k^{-1}{\mathbf{B}}_k,{\mathbf{Q}}_k=(e^{{\mathbf{A}}_k\mathrm{\Delta }t}-\mathbf{I}){\mathbf{A}}_k^{-1}$ ${\mathbf{A}}_k=\frac{1}{\mathrm{\Delta }t}{\int }_{t_k}^{t_{k+1}}\mathbf{A}(t)\mathrm{dt},{\mathbf{B}}_k=\frac{1}{\mathrm{\Delta }t}{\int }_{t_k}^{t_{k+1}}\mathbf{B}(t)\mathrm{dt}$ . ${\tau }_k=\mathrm{int}(0.5+\tau /\mathrm{\Delta }t)$ , in which int is the rounding down operation.

Here, the above equation can be further transformed into a discrete map, as shown below.

$${\mathbf{V}}_{k+1}={\mathbf{Z}}_k{\mathbf{V}}_k+{\mathbf{E}}_k$$

(20)

Where ${\mathbf{V}}_k=[O_{\mathrm{rx},k},O_{\mathrm{ry},k},{\stackrel{·}{O}}_{\mathrm{rx},k},{\stackrel{·}{O}}_{\mathrm{ry},k},{\stackrel{·}{O}}_{\mathrm{rx},k-1},{\stackrel{·}{O}}_{\mathrm{ry},k-1},\cdots ,{\stackrel{·}{O}}_{\mathrm{rx},k-N_m},{\stackrel{·}{O}}_{\mathrm{ry},k-N_m}]$ ${\mathbf{Z}}_k=\left[\begin{array}{ccccccc}{\mathbf{P}}_k\hfill & \cdots \hfill & 0.5{\mathbf{R}}_k\hfill & 0.5{\mathbf{R}}_k\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\mathbf{I}}& \cdots \hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & \mathbf{I}\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & \mathbf{I}\hfill & 0\hfill \end{array}\right]$ . ${\mathbf{E}}_k={\left[{\mathbf{Q}}_k{F}_k{0}_{1\times 2N_m}\right]}^T$ .

The state transition relationship of the system in one cutter rotation period can be expressed as

$${\mathbf{V}}_{N_m+1}=\mathbf{\Phi }{\mathbf{V}}_1+\mathbf{G}$$

(21)

Where $\mathbf{\Phi }={\mathbf{Z}}_{N_m}{\mathbf{Z}}_{N_m-1}...{\mathbf{Z}}_1$ is the state transformation matrix, and $\mathbf{G}={\mathbf{E}}_{N_m}+\sum _{i=1}^{N_m}{\mathbf{Z}}_{i+1}{\mathbf{Z}}_{i+2}\cdots {\mathbf{Z}}_{N_m}{\mathbf{E}}_i$ .

According to the Floquet theory, the system is stable only if the modulus of all eigenvalues of the state transformation matrix is less than or equal to one. Hence, the constraint of milling stability can be expressed as

$$|{\lambda }_{\max }(\mathit{x})|\le 1$$

(22)

Where $|{\lambda }_{\max }(\mathit{x})|$ is the maximal modulus of all eigenvalues of the transformation matrix.

Machining quality constraint

When the milling process is stable, the dynamical response of the cutter can be obtained by calculating the system fixed point. Let $\mathbf{V}*={\mathbf{V}}_{N_m+1}={\mathbf{V}}_1$ in equation (21), the system response of the dynamic milling system becomes

$$\mathbf{V}*=(\mathbf{I}-\mathbf{\Phi }{)}^{-1}\mathbf{G}$$

(23)

In this case, the transient displacement response of the cutter can be substituted into the following equation to obtain the instantaneous trajectory surface of the cutting point P.

$$\left[\begin{array}{c}{P}_{x,i,j}\hfill \\ \multicolumn{1}{c}{P_{y,i,j}}\end{array}\right]=\left[\begin{array}{c}{o}_{\mathrm{rx}}(t)\hfill \\ \multicolumn{1}{c}{o_{\mathrm{ry}}(t)}\end{array}\right]+\left[\begin{array}{c}{R}_{i,j}\sin ({\varphi }_{i,j}(t))\hfill \\ \multicolumn{1}{c}{R_{i,j}\cos ({\varphi }_{i,j}(t))}\end{array}\right]$$

(24)

Where ${[o_{\mathrm{rx}},o_{\mathrm{ry}}]}^T=[(x_{k-N_m},x_{k-N_m+1},\cdots x_k);(y_{k-N_m},y_{k-N_m+1},\cdots ,y_k){]}^T$ .

As shown in Figure 3, the machined surface contour $S_{u,v}^m$ can be obtained by finding the intersection of the actual tool tip trajectory envelope surface and the workpiece entity, then the machining error can be expressed as the absolute average value of the deviation between the actual machining surface and the ideal design surface $S_{u,v}^d$ , as shown in the following formula.

$$E_s=\frac{1}{N_u{N}_v}\sum _u\sum _v|S_{u,v}^m-S_{u,v}^d|$$

(25)

Figure 3.

Definition diagram of the machining quality constraint.

The constraint condition of the surface contour error can be achieved as

$$E_s(\mathit{x})\le E_s^d$$

(26)

Where $E_s^d$ is the allowable value of contour error.

Then, the roughness would be calculated by

$$R_a=\frac{1}{N_u}\sum _{u*}|S_{u*,v}^m-{\overline{S}}_{u*,v}^m|$$

(27)

Correspondingly, the constraint condition of the roughness is defined as

$$R_a(\mathit{x})\le R_a^d$$

(28)

Where $R_a^d$ is the allowable value of roughness.

Mathematical model of the optimization problem

Finally, considering all above multiple constraints together, the mathematical model of the machining parameter optimization problem can be formed. It becomes

$$\begin{array}{c}min-f_{\mathrm{MRR}}(\mathit{x})\hfill \\ s.t.\{\begin{array}{c}\mathit{x}=[n,a_e,a_p,N_t,f_t{]}^T\hfill \\ \multicolumn{1}{c}{n^l\le n\le n^u}\\ \multicolumn{1}{c}{a_e^l\le a_e\le a_e^u}\\ \multicolumn{1}{c}{a_p^l\le a_p\le a_p^u}\\ \multicolumn{1}{c}{f_t^l\le f_t\le f_t^u}\\ \multicolumn{1}{c}{N_t^l\le N_t\le N_t^u}\\ \multicolumn{1}{c}{max(F_c(\mathit{x},t))\le F_{\mathrm{clim}}}\\ \multicolumn{1}{c}{|{\lambda }_{\max }(\mathit{x})|\le 1}\\ \multicolumn{1}{c}{E_s(\mathit{x})\le E_s^d}\\ \multicolumn{1}{c}{R_a(\mathit{x})\le R_a^d}\end{array}\hfill \end{array}$$

(29)

Solution of the optimization model

Uncertainty treatment of model parameters

Because there is a certain degree of uncertainty for the model parameters of the actual dynamic system, which will cause some deviation of the constraint critical value between the actual condition and the real one. Therefore, it is necessary to give a certain safety margin to the constraint conditions. At this time, the above formula is improved into

$$\begin{array}{c}min-f_{\mathrm{MRR}}(\mathit{x})\hfill \\ s.t.\{\begin{array}{c}\mathit{x}=[n,a_e,a_p,N_t,f_t{]}^T\hfill \\ \multicolumn{1}{c}{n^l\le n\le n^u}\\ \multicolumn{1}{c}{a_e^l\le a_e\le a_e^u}\\ \multicolumn{1}{c}{a_p^l\le a_p\le a_p^u}\\ \multicolumn{1}{c}{f_t^l\le f_t\le f_t^u}\\ \multicolumn{1}{c}{N_t^l\le N_t\le N_t^u}\\ \multicolumn{1}{c}{max(F_c(\mathit{x},t))\le {\lambda }_F{F}_{\mathrm{clim}}}\\ \multicolumn{1}{c}{|{\lambda }_{\max }(\mathit{x})|\le {\lambda }_S}\\ \multicolumn{1}{c}{E_s(\mathit{x})\le {\lambda }_E{E}_s^d}\\ \multicolumn{1}{c}{R_a(\mathit{x})\le {\lambda }_R{R}_a^d}\end{array}\hfill \end{array}$$

(30)

Where ${\lambda }_U(U=F,S,E,R)$ is a real value less than one. The smaller the value is, the greater the safety margin of the constraint has. The range of its values suggested in this paper is 0.7–0.98.

f_t Optimization based on the golden section method

According to the prior knowledge, the roughness of milling machined surface is closely related to the feed per tooth, and is less affected by other machining parameters. Consequently, in order to improve the machining efficiency, it is often expected to increase the feed per tooth as much as possible on the premise of satisfying the roughness constraint. The roughness constraint can be dropped mathematically from equation (30) as following.

$$|R_a(\mathit{x})-{\lambda }_R{R}_a^d|\le {\epsilon }_{R_a}(\mathit{x}=[n^l,a_e^u,a_p^u,N_t^l,f_t{]}^T)$$

(31)

It can be judged that the above function is convex, and its solving can be applied by the golden section method. The optimization procedure is as follows:

• Step (1): Initialize the optimization variable $\mathit{x}=[n^l,a_e^u,a_p^u,N_t^l,f_t{]}^T$ , set the iteration number $k=1$ . Preset the search interval of $f_t$ as $a_1=f_t^l,b_1=f_t^u$ and the convergence accuracy ${\epsilon }_{R_a}$ . Based on the principle of golden section method, let the golden section ratio $\eta =0.618$ .

• Step (2): Calculate ${\alpha }_1=b_1-\eta (b_1-a_1)$ and ${\beta }_1=a_1+\eta (b_1-a_1)$ , and assign ${\alpha }_1$ and ${\beta }_1$ to $f_t$ respectively, and calculate the objective function $|R_a(\mathit{x})-{\lambda }_R{R}_a^d|$ .

• Step (3): Let $k=k+1$ . Reduce the search interval according to the principle of interval elimination: if $|R_a(\mathit{x}({\alpha }_k))-{\lambda }_R{R}_a^d|\ge |R_a(\mathit{x}({\beta }_k))-{\lambda }_R{R}_a^d|$ , then let $a_{k+1}={\alpha }_k,b_{k+1}=b_k$ , ${\alpha }_{k+1}={\beta }_k$ , and ${\beta }_{k+1}=a_{k+1}+\eta (b_{k+1}-$ $a_{k+1})$ . Otherwise, let $a_{k+1}=a_k,b_{k+1}={\beta }_k$ , ${\beta }_{k+1}={\alpha }_k$ , and ${\alpha }_{k+1}=b_{k+1}-\eta (b_{k+1}-a_{k+1})$ .

• Step (4): If satisfying $|R_a(\mathit{x})-{\lambda }_R{R}_a^d|\le {\epsilon }_{R_a}$ , let $f_t^{*}=(a_{k+1}+b_{k+1})/2$ and terminate the iteration. Otherwise, go to step (3).

[n,a_e,a_p]^T Optimization based on the random vector search method

After the $f_t$ optimization is finished, the problem in equation (30) can be further described as following.

$$\begin{array}{c}min-f_{\mathrm{MRR}}(\mathit{x})\hfill \\ s.t.\{\begin{array}{c}\mathit{x}=[n,a_e,a_p,N_t^l,f_t^{*}{]}^T\hfill \\ \multicolumn{1}{c}{g_1=n^l-n\le 0}\\ \multicolumn{1}{c}{g_2=n-n^u\le 0}\\ \multicolumn{1}{c}{g_3=a_e^l-a_e\le 0}\\ \multicolumn{1}{c}{g_4=a_e-a_e^u\le 0}\\ \multicolumn{1}{c}{g_5=a_p^l-a_p\le 0}\\ \multicolumn{1}{c}{g_6=a_p-a_p^u\le 0}\\ \multicolumn{1}{c}{g_7=max(F_c(\mathit{x},t))-{\lambda }_F{F}_{\mathrm{clim}}\le 0}\\ \multicolumn{1}{c}{g_8=|{\lambda }_{max}(\mathit{x})|-{\lambda }_S\le 0}\\ \multicolumn{1}{c}{g_9=E_s(\mathit{x})-{\lambda }_E{E}_s^d\le 0}\end{array}\hfill \end{array}$$

(32)

In order to solve the multiple constraints optimization problem, the random vector search method is employed here. As shown in Figure 4, the basic idea of this method is to generate multiple initial samples randomly in the feasible region, then for the sth sample ${\mathit{x}}_s$ , let it search forward with a random vector. In its the kth iteration, there would generate a random vector $V_{s,k,q}$ at the qth exploratory search, and we can obtain a new sample in the (k+1)th iteration by ${\mathit{x}}_{s,k+1}={\mathit{x}}_{s,k}+V_{s,k,q}$ . Then one need judge whether all constraints are satisfied and the objective function value has decreased. If satisfied, let $k=k+1$ and iterate forward. Otherwise, let $q=q+1$ , continue to generate random vectors $V_{s,k,q+1}$ at ${\mathit{x}}_{s,k}$ , and keep searching. The whole optimization procedure of multi-sample random vector search is as follows:

• Step (1): Determine the optimization objective function $-f_{\mathrm{MRR}}(\mathit{x})$ , variable $\mathit{x}$ , and constraints $g_j\left(\mathit{x}\right)$ . Randomly generate initial samples of optimization variable, and only retain the individual which meet the constraints. Finally we can form a sample set $\left\{{\mathit{x}}_s\right\}$ with a number of $N_s$ individuals.

• Step (2): Extract the sth sample ${\mathit{x}}_s$ from the sample set.

• Step (3): Further extract the kth variable ${\mathit{x}}_{s,k}=[n_{s,k},a_{e,s,k},a_{p,s,k}]$ .

• Step (4): On its qth search, generate a random vector as following.

$$\{\begin{array}{c}\delta n_{s,k,q}={\alpha }_{s,k,q}^1(n^u-n^l)\hfill \\ \delta a_{e,s,k,q}={\alpha }_{s,k,q}^2(a_e^u-a_e^l)\hfill \\ \delta a_{p,s,k,q}={\alpha }_{s,k,q}^3(a_p^u-a_p^l)\hfill \\ V_{s,k,q}=[\delta n_{s,k,q},\delta a_{e,s,k,q},\delta a_{p,s,k,q}]\hfill \end{array}$$

(33)

Where ${\alpha }_{s,k,q}^r(r=1,2,3)$ is the random coefficient, which is taken the number of interval [−1, 1].

• Step (5): Calculate the (k+1)th value of the sample ${\mathit{x}}_s$ by

$$\{\begin{array}{c}{\mathit{x}}_{s,k+1}={\mathit{x}}_{s,k}+V_{s,k,q}\hfill \\ \multicolumn{1}{c}{n_{s,k+1}=\min (\max (\lfloor n_{s,k+1}\rfloor ,n^l),n^u))}\\ \multicolumn{1}{c}{a_{e,s,k+1}=\min (\max (\lfloor 10a_{e,s,k+1}\rfloor /lfloor;\lfloor 10a_{e,s,k+1}\rfloor 10,a_e^l),a_e^u))10,a_e^l),a_e^u))}\\ \multicolumn{1}{c}{a_{p,s,k+1}=\min (\max (\lfloor 10a_{p,s,k+1}\rfloor /lfloor;\min (\max (\lfloor 10a_{p,s,k+1}\rfloor 10,a_p^l),a_p^u))10,a_p^l),a_p^u))}\end{array}$$

(34)

Where ⌊⌋ is the rounding down operation. After this operation, one can make sure that the spindle speed is integer, and both the radial and axial depth of cuts have only one decimal.

• Step (6): Substitute the obtained parameters ${\mathit{x}}_{s,k+1}$ into equation (32), and further make a judgment by

$$\{\begin{array}{c}\mathrm{if}k+1<N_k\hfill \\ \multicolumn{1}{c}{\mathrm{if}{g}_j\left(\mathit{x}|{\mathit{x}}_{s,k+1}\right)\le 0\mathrm{and}{f}_{\mathrm{MRR}}({\mathit{x}}_{s,k+1})-f_{\mathrm{MRR}}({\mathit{x}}_{s,k})>0}\\ \multicolumn{1}{c}{{\mathit{x}}_{s,k}={\mathit{x}}_{s,k+1},k=k+1,\mathrm{goto}\mathrm{step}(3);}\\ \multicolumn{1}{c}{\mathrm{else}}\\ \multicolumn{1}{c}{{\mathit{x}}_{s,k}={\mathit{x}}_{s,k},q=q+1,\mathrm{goto}\mathrm{step}(4);}\\ \multicolumn{1}{c}{\mathrm{else}}\\ \multicolumn{1}{c}{\mathrm{stop};}\\ \multicolumn{1}{c}{\mathrm{end}}\end{array}$$

Then, one can obtain the optimal result ${\mathit{x}}_s^{*}={\mathit{x}}_{s,N_k}$ of the sth sample, where $N_k$ is the maximal search number.

• Step (7): Let $s=s+1$ . If $s<N_s$ , go to step (2). Otherwise, finish the iteration. Finally, compare the objective function of all samples $-f_{\mathrm{MRR}}({\mathit{x}}_s^{*})$ to find the minimal one, and let this one be as the optimal individual ${\mathit{x}}^{*}={\min (-f_{\mathrm{MRR}}({\mathit{x}}_s^{*}))|}_{{\mathit{x}}_s^{*}}$ , which corresponds to the optimization result.

Figure 4.

Schematic diagram of random vector search.

Overall optimization procedure

Figure 5 summarizes the overall optimization flowchart. At the beginning of optimization, the tool geometry parameters, dynamics parameters, cutter runout, and specific cutting force coefficient are initialized, and the depth-first or width-first tool path topology is selected. Later, the machining parameters are optimized in a four layer cycle iteration. The first layer is to optimize the discrete number of allowance under depth-first or width-first situation. The second layer is to optimize the tooth number of cutter. The third layer is to optimize the feed per tooth by using the golden section method. The fourth layer is to optimize the spindle speed, radial, and axial depth of cuts combination with the aid of the random vector search method. After the four layer iteration calculation, the machining parameters with the maximum milling efficiency are selected as the optimal results.

Figure 5.

Overall optimization flowchart.

Experimental verification

Experimental setup

In order to verify the validity of the proposed machining parameter optimization method, the paper firstly carried out a few parameter optimization cases that performed under the given boundary conditions, then conducted a series of milling experiments on the five-axis high-speed machining center DMG DMU50 to verify the theoretical analysis results. As shown in Figure 6, the workpiece material was aluminum alloy 7050 with a machining allowance of L_w = 3 mm, L_h = 50 mm, and L_l = 110 mm. The machining tool path topology was selected as depth-first. In the tests, a series of cutters with different tooth number, a diameter of 16 mm, a helix angle of 45°, and an overhung length of 85 mm were employed. The cutter had runout with offset value 8 µm and its angle 25.6°, and inclination value 1.1 × 10⁻⁴° and its angle 46.5°. The milling process had no coolant, and the specific cutting force coefficients of the cutter-workpiece combination were calibrated as k_ts = 999 N/mm², k_rs = 327 N/mm², k_as = 357 N/mm², k_tp = 30 N/mm, k_rp = 32 N/mm, k_ap = 0.4 N/mm. Along the X and Y directions in machining coordinate system, the modal masses of the tool tip were set as 0.407 and 0.655 kg, the damping ratios were set as 0.035 and 0.028, and the natural frequencies were set as 959.4 and 922.1 Hz, respectively, through the modal impacting and identification experiment. During the milling process, the transient cutting force signal had been measured by a 9265B Kistler dynamometer with a sampling frequency of 10 kHz. After machining, the Taylor Hobson aspheric profile tester was used to measure the surface roughness topography of the machined surface with a resolution of 0.125 µm, and the global classic SR575 three-coordinate measuring machine was used to test the contour accuracy of the part with a resolution of 1.8 µm.

Figure 6.

Experimental setup.

Verification cases

In terms of the constraints, the parameter feasible region of cutter tooth is set as $[N_t^l,N_t^u]=[1,4]$ , and that of milling parameters are preset to $[n^l,n^u]=[5,10]\mathrm{krpm}$ , $[a_e^l,a_e^u]=[0.01,3]\mathrm{mm}$ , $[a_p^l,a_p^u]=[0.01,50]\mathrm{mm}$ , and $[f_t^l,f_t^u]=[0.01,0.4]\mathrm{mm}/\mathrm{tooth}$ , respectively. The milling type is selected as down milling. Additionally, the allowable maximum milling force is 5.0 kN, and the allowable values of the machining contour error and roughness are assigned to 0.04 mm and 0.4 µm, respectively. The safety coefficients of milling force, stability, machining contour error and roughness are chosen as 0.95 and 0.92, 0.95 and 0.9, respectively. As depth-first type is adopted in the machining tool path topology, the discrete quantity of width direction is selected as $N_w=1,2,3$ . The default sample number and iteration step are set as 30 and 60, respectively.

In the total optimization procedure, the feed per tooth is firstly optimized under the constraint of roughness, and it can be concluded that the optimal feed per tooth for the single-tooth cutter, two-tooth cutter, three-tooth cutter, and four-tooth cutter are 0.275, 0.15, 0.10, 0.075 mm/tooth, respectively. Then, the optimization of other milling parameters for the cutters with different tooth number are further performed. Figure 7 demonstrates the MRR iteration process of all samples under different tooth numbers and radial depth of cuts. Among them, Figure 7(a)–(d) are the MRR variation trend of all 30 samples in total 60 iteration steps for the single-tooth, two-tooth, three-tooth, and four-tooth cutter with a radial depth of cut of 3 mm (N_w = 1), respectively. Figure 7(e) and (f) are the MRR variation trend of all 30 samples in total 60 iteration steps for the three-tooth cutter with a radial depth of cut of 1 mm (N_w = 3) and 1.5 mm (N_w = 2), respectively. As can be seen that the MRR of all samples keeps increasing during the iteration process. After the iteration completed, the MRR of the 28th individual of the three-tooth cutter with a radial depth of cut 3 mm (N_w = 1) in Figure 7(c) reaches to the global maximum value, which is equal to 131.2 cm³/min. It shows that this individual is the best individual of all samples. From the comparison between Figure 7(c) and (d), one can still realize that more tooth number may not obtain a higher machining efficiency, especially under the multiple constraint conditions.

Figure 7.

MRR iteration process of all samples under different tooth numbers and radial depth of cuts: (a) MRR iteration result (N_t = 1, a_e = 3 mm), (b) MRR iteration result (N_t = 2, a_e = 3 mm), (c) MRR iteration result (N_t = 3, a_e = 3 mm), (d) MRR iteration result (N_t = 4, a_e = 3 mm), (e) MRR iteration result (N_t = 3, a_e = 1 mm), and (f) MRR iteration result (N_t = 3, a_e = 1.5 mm).

Figure 8 detailly illustrates the evolution process of the above optimal individual in the whole iteration calculation. Among them, Figure 8(a) shows the initial and final spindle speed and axial depth of cut values of the total 30 samples, as well as its iterative values of the 28th individual. It is obvious that all samples are satisfied the requirement of milling stability. Besides, the distribution of the samples before optimization is scattered, while has a stronger aggregation after optimizing. Figure 8(b) further plots the variation of spindle speed and axial depth of cut of the 28th individual during total iteration, and the optimized spindle speed and axial depth of cut are obtained as 9411 rpm and 15.6 mm, respectively. Besides, Figure 8(c)–(e) show the variation of the machining contour error, maximum milling force, and maximum eigenvalue of the milling stability transfer matrix with the iteration number of the 28th individual. What we can see more is that the above three physical parameters always satisfy its corresponding constraint conditions. The machining contour error is very close to its allowance value in Figure 8(e), which indicates that the optimization iteration has well been done. Figure 8(f) compares the change of MRR before and after optimization. Among them, the average MRR of all samples after optimization has increased by 216.7% compared to that of before optimization, and the maximum MRR of all samples after optimization has increased about 53.6% than before. More importantly, the optimized MRR of the best individual is 659.2% higher than that of initial, indicating that the parameter optimization process can greatly improve the machining efficiency.

Figure 8.

Evolution process of the optimal individual parameters: (a) iteration result of n and a_p, (b) iteration process of n and a_p, (c) iteration process of max(F_c), (d) iteration process of λ_max, (e) iteration process of E_s, and (f) MRR comparison.

The paper makes further efforts on the milling experiments by using the optimal machining parameters $[n,a_e,a_p,N_t,f_t]=[9411\mathrm{rpm},3\mathrm{mm},15.6\mathrm{mm},3,0.1\mathrm{mm}/\mathrm{tooth}]$ for verification. Under the existing machining allowance, the cutter needs feed along the radial direction for one time, and feed along the axial direction for three times. The axial depth of cut in the root cleaning is 3.2 mm.The measured physical parameters of the optimal individual are shown in Figure 9. Among them, Figure 9(a) is a comparison diagram of the predicted and measured roughness values of the machined surface. Here, its measured value is about 0.327 µm, which perfectly satisfies the roughness constraint. Figure 9(b) brings insight into the measured results of the milling force, and the maximum value of the resultant force is 743.3 N, which also meets its constraints. As shown in Figure 9(c), it is obvious that the frequency spectrum of measured milling force only contains the spindle rotation frequency and its high order harmonic components, without any chatter frequency. The result indicates that the milling process is stable. Lastly, Figure 9(d) shows the comparison between the predicted and the measured value of the machined surface profile. Here, the machining contour error value is below zero, which means that overcut phenomenon occurs at this time. The actual measured value of the contour error is 0.046 mm, which is slightly larger than the constraint value of 0.038 mm. It is mainly caused by cutter vibration during milling process and the parameter uncertainty in the dynamics model. Taking into account the existence of these disturbance factors, the optimization result is accurate and effective.

Figure 9.

Measured values of the optimal individual parameters: (a) roughness comparison, (b) measured milling force, (c) frequency spectrum of milling force, and (d) machining error comparison.

Conclusion

In this paper, an efficient parameter optimization method for finish milling machining with chatter free has been proposed. Through the theoretical research and experimental verification, the results show that the method can greatly improve the milling efficiency and realize chatter free finish machining under multiple constraints. In the optimization model, multiple constraints such as parameter feasible region, milling force, milling stability, machining contour accuracy, and roughness are considered, and multiple variables such as tooth number, feed per tooth, spindle speed, radial, and axial depth of cuts are both studied. Concerning the uncertainty of model parameters, safety margin coefficients are used to enhance the adaptability of the model. As the roughness is mainly affected by the feed per tooth, it is proposed to use the golden section method to solve the single variable optimization problem. Furthermore, when solving the optimization problem with the remaining variables, a numerical random vector search optimization method is proposed, which can ensure a fast iteration and realize multi-variable optimization of spindle speed, radial, and axial depth of cuts. Based on the proposed four level cycle optimization strategy, the paper can determine the optimal individual in all samples and select the machining parameter combination with highest efficiency at last.