# Improving Stability Prediction in Peripheral
## 2 Materials and Methods  
### Milling Stability Prediction with the EMHP Method  
In order to predict stability in machining operations, in this paper we focus on the manufacturing of a thin-walled workpiece through peripheral milling operations. The workpiece has a single degree of freedom (flexibility only in the \(y\)-direction according to the machine tool standard axis nomenclature), with a dynamic model of the following form:  
\[\ddot{y}(t)+2\zeta\omega_{n}\dot{y}(t)+\omega^{2}{}_{n}y(t)=\frac{F_{y}}{m_{m}}( y(t)-y(t-\tau)) \tag{1}\]  
Here, \(m_{m}\) represents the modal mass, \(\zeta\) is the damping ratio, \(\omega_{n}\) is the natural angular frequency, and \(F_{y}\) is the cutting force over the workpiece in the \(y\)-direction, which is given by:  
\[F_{y}(t)=\sum_{j=1}^{z_{n}}g\big{(}\phi_{j}(t)\big{)}\cos\big{(}\phi_{j}(t) \big{)}\left(-K_{tc}\sin\left(\big{(}\phi_{j}(t)\big{)}\right)+K_{rc}\cos\left( \big{(}\phi_{j}(t)\big{)}\right)\right) \tag{2}\]  
where \(z_{n}\) is the number of teeth, \(K_{tc}\) and \(K_{rc}\) are the tangential and the normal cutting force coefficients, and \(\phi_{j}\left(t\right)\) is defined as:  
\[\phi_{j}(t)=(2\pi n/60)t+j2\pi/z_{n} \tag{3}\]  
where \(n\) is the spindle speed in rpm. The function \(g(\phi_{j}(t))\) is a screen function, and is equal to 1 if the tooth \(j\) is in the cut, and is equal to 0 if \(j\) is out of the cut. To obtain the stability lobes of Equation (1), we proposed a solution procedure based on the EMHP algorithm [18], which relies on a sequence of subintervals that provide approximate solutions requiring less CPU time in comparison with other methods found in the literature. In this methodology, we rewrote Equation (1) as:  
\[\ddot{y}_{i}(T)+2\zeta\omega_{n}\dot{y}_{i}(T)+\omega_{n}y_{i}(T)\approx \frac{F_{yt}}{m_{m}}(y_{i}(T)-y_{i}^{\tau}(T)) \tag{4}\]  
where \(F_{yt}=F_{y}(t)\), and \(y_{i}(T)\) denotes the approximate solution of order \(m\) in the \(i\)-th sub-interval that satisfies the initial conditions \(y_{i}(0)=y_{i-1}\), \(\dot{y}_{i}(0)=\dot{y}_{i-1}\). The Equation (1) can be written in state space, in accordance with the following matrix-form representation:  
\[\dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t)+\mathbf{B}(t)(\mathbf{x}(t)- \mathbf{x}(t-\tau)) \tag{5}\]  
where \(x=\big{[}x,\dot{x}\big{]}^{T}\), \(\mathbf{A}(t-\tau)=\mathbf{A}(t)\), \(\mathbf{B}(t-\tau)=\mathbf{B}(t)\), and \(\tau\) is the time delay. By following the homotopic procedure, the equivalent form of Equation (5) can be written as:  
\[\dot{\mathbf{x}}_{i}(T)-\mathbf{A}_{i}\mathbf{x}_{i}(T)\approx\mathbf{B}_{i} \mathbf{x}_{i}^{\tau}(T) \tag{6}\]  
where \(\mathbf{x}_{i}(T)\) denotes the \(m\) order solution of Equation (6) in the \(i\)-th sub-interval that satisfies the initial conditions \(\mathbf{x}_{i}(0)=\mathbf{x}_{i-1}\), \(\mathbf{A}_{i}\), and \(\mathbf{B}_{i}\) represents the values of the periodic coefficients at the time \(t\). In order to approximate the delay term \(\mathbf{x}_{i}^{\tau}(T)\) in Equation (6), the period [\(t_{0}-t_{0}\)] is discretized into \(N\) equally-spaced points, but the method does not accept strictly-spaced sub-intervals. Thus, each sub-interval has a time span equal to \(\Delta t=\tau/(N-1)\). Here, we assumed that the function \(\mathbf{x}_{i}^{\tau}(T)\) in the delay sub-interval [\(t_{i-N}\), \(t_{i-N+1}\)] had a first-order polynomial representation, of the form:  
\[\mathbf{x}_{i}^{\tau}(T)=\mathbf{x}_{i-N+1}(T)\approx\mathbf{x}_{i-N}+\frac{N -1}{\tau}(\mathbf{x}_{i-N+1}-\mathbf{x}_{i-N})T \tag{7}\]  
To simplify the notation, we let \(\mathbf{x}_{i}\equiv\mathbf{x}_{i}(T_{i})\).