## 3 Stability Analysis and Experimental Validation  
### Stability Analysis for Thin-Floor Machining  
In finish milling operations of thin-floor components where the dominant mode is oriented in the \(z\)-direction, the stability model can be simplified considering only the Frequency Response Function (FRF) of the thin floor, so the system becomes a one degree of freedom (1 dof) model as described by [PERSON] et al. [30]. In this case, the dynamic chip thickness is determined from:  
\[h\Big{(}\phi_{j}(t)\Big{)}=\cos\kappa\Delta z \tag{8}\]  
where \(\phi_{j}(t)\) is the angular position of the cutting edge \(j\) and \(\Delta z=(z(t)-z(t-\tau))\). The equation of motion that describes the dynamics of the thin floor is given as:  
\[\ddot{z}(t)+2\zeta\omega_{n}\dot{z}(t)+\omega_{n}^{2}z(t)=\frac{F_{z}(t)}{m_{ m}}(z(t)-z(t-\tau)) \tag{9}\]  
where \(m_{m}\) is the modal mass, \(\zeta\) describes the damping ratio, \(\omega_{n}\) is the angular natural frequency, and \(F_{z}\) is the cutting force in the \(z\)-direction, which is given as:  
\[F_{z}(t)=\sum\
olimits_{j=0}^{N_{t}-1}g\Big{(}\phi_{j}(t)\Big{)}h\Big{(}\phi_ {j}(t)\Big{)}\Big{(}-K_{rc}\cos^{2}\kappa-K_{ac}\sin\kappa\cos\kappa\Big{)}. \tag{10}\]  
Here, \(g\Big{(}\phi_{j}(t)\Big{)}\) is a switching function that returns the value of one when the tooth cutter is cutting and the value of zero otherwise. The matrix representation of Equation (9) becomes:  
\[\dot{\mathbf{z}}(t)=\mathbf{A}(t)\mathbf{z}(t)+\mathbf{B}(t)\mathbf{z}(t-\tau) \tag{11}\]  
where \(\mathbf{z}=\left[z,\dot{z}\right]^{T}\), \(\mathbf{A}(t+\tau)=\mathbf{A}(t)\), \(\mathbf{B}(t+\tau)=\mathbf{B}(t)\), and \(\tau\) is the time delay. Following the Enhanced Multistage Homotopy Perturbation Method (EMHPM) described in [34], Equation (11) is written in an equivalent form as:  
\[\dot{\mathbf{z}}_{i}(T)-\mathbf{A}_{i}\mathbf{z}_{i}(T)\approx\mathbf{B}_{i} \mathbf{z}_{i}^{\tau}(T) \tag{12}\]  
Figure 4: Comparison between experiment data and the proposed force model at 3000 rev/min and the cutting conditions: (**a**) \(a_{p}=0.5\) mm,\(f_{z}=0.2\) mm/tooth; (**b**) \(a_{p}=1.5\) mm,\(f_{z}=0.15\) mm/tooth. where \(\mathbf{z}_{i}(T)\) indicates the m-order solution for Equation (11) that satisfies the initial conditions \(\mathbf{z}_{i}(0)=\mathbf{z}_{i-1}\), \(\mathbf{A}_{t}\) and \(\mathbf{B}_{t}\) are the periodic matrix whose values vary with time \(t\). The period \([t_{0}-\tau,t_{0}]\) is discretized in \(N\) points to approximate the term \(\mathbf{z}_{i}^{\tau}(T)\) to the delay in Equation (12). The function \(\mathbf{z}_{i}^{\tau}(T)\) could be approximate as a first-polynomial representation:  
\[\mathbf{z}_{i}^{\tau}(T)\approx\mathbf{z}_{i-N}+\frac{N-1}{\tau}(\mathbf{z}_{ i-N+1}+\mathbf{z}_{i-N})T \tag{13}\]  
Next, we define \(\mathbf{z}_{i}\equiv\mathbf{z}_{i}(T_{i})\) and then, Equation (12) is substituted into Equation (13) to get:  
\[\dot{\mathbf{z}}_{i}=\mathbf{A}_{t}\mathbf{z}_{i}+\mathbf{B}_{t}\mathbf{z}_{i -N}-\frac{N-1}{\tau}\mathbf{B}_{t}\mathbf{z}_{i-N}T+\frac{N-1}{\tau}\mathbf{B} _{t}\mathbf{z}_{i-N+1}T \tag{14}\]  
where  
\[\mathbf{A}_{t}=\left[ \begin{array}{cc}0&1\\ -{\omega_{n}}^{2}+\frac{F_{z}(t)}{m_{m}}&-2\zeta\omega_{n}\end{array} \right] \mathbf{B}_{t}=\left[ \begin{array}{cc}0&0\\ -\frac{F_{z}(t)}{m_{m}}&0\end{array} \right], \tag{15}\]  
To compute the stability lobes of Equation (14), the procedure described in [26] is followed. Then, it is assumed that Equation (14) can be written as a function of discretized states:  
\[\mathbf{z}_{i}(T)\approx\mathbf{P}_{i}(T)\mathbf{z}_{i-1}+\mathbf{Q}_{i}(T) \mathbf{z}_{i-N}+\mathbf{R}_{i}(T)\mathbf{z}_{i-N+1} \tag{16}\]  
where  
\[\mathbf{P}_{i}(T)=\sum\limits_{k=0}^{m}\frac{1}{k}\mathbf{A}_{t}^{k}T^{k},\]  
\[\mathbf{Q}_{i}(T)=\left\{ \begin{array}{cc}\sum\limits_{k=1}^{m}\frac{N-1}{( k+1)(\tau)}\mathbf{A}_{t}^{k-1}\mathbf{B}_{t}T^{k+1}&m\geq 1\\ 0&m=0\end{array} \right. \tag{17}\]  
\[\mathbf{R}_{i}(T)=\left\{