*Article* **Uncharted Stable Peninsula for Multivariable Milling Tools by High-Order Homotopy Perturbation Method**

**Jose de la Luz Sosa 1, Daniel Olvera-Trejo 1,\*, Gorka Urbikain 2,\*, Oscar Martinez-Romero 1, Alex Elías-Zúñiga <sup>1</sup> and Luis Norberto López de Lacalle <sup>2</sup>**


Received: 9 October 2020; Accepted: 3 November 2020; Published: 6 November 2020

**Abstract:** In this work, a new method for solving a delay differential equation (DDE) with multiple delays is presented by using second- and third-order polynomials to approximate the delayed terms using the enhanced homotopy perturbation method (EMHPM). To study the proposed method performance in terms of convergency and computational cost in comparison with the first-order EMHPM, semi-discretization and full-discretization methods, a delay differential equation that model the cutting milling operation process was used. To further assess the accuracy of the proposed method, a milling process with a multivariable cutter is examined in order to find the stability boundaries. Then, theoretical predictions are computed from the corresponding DDE finding uncharted stable zones at high axial depths of cut. Time-domain simulations based on continuous wavelet transform (CWT) scalograms, power spectral density (PSD) charts and Poincaré maps (PM) were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool.

**Keywords:** chatter; multivariable tool; stable peninsula; homotopy perturbation method

#### **1. Introduction**

There are many phenomena in different fields of science and engineering where the physical response of a variable involves not only the value at time *t* but also the effects that occur in an earlier state *t* − τ. Thus, delay systems appear in many engineering problems, such as in the shimmy effect (wheel vibration) [1], vehicle traffic models [2], feedback stabilization problems [3] and in the regenerative vibration of machine-tools better known as chatter [4]. In cases where the net force depends on the current values and some past values (history) such as position and speed, the system dynamic behavior can be modeled using a differential delay equation (DDE).

It is well-known that during a milling process, unstable vibrations also known as self-excited vibration or chatter may occur. Chatter reduces the machining efficiency due to low material removal rate by reducing the workload and affects surface quality, shortens tool life and accelerates tool wear. Researchers are studying several ways to overcome this limitation. Kuljanic et al. [5] studied the incorporation of a chatter detection system based on multiple sensors to milling operations for industrial conditions, Zhuo et al. [6] used a method based on fractal dimension for the flank milling of a thin-walled blade, which can reflect the chatter severity level through the morphological change in signal. Paul and Morales [7], to mitigate chatter, presented an active controller based on the technique of discrete time sliding mode control (DSMC) blended with the type-2 fuzzy logic system. Moreover, Peng et al. [8] presented a method based on a dynamic cutting force simulation model and a machine learning approach based on statistical learning theory to predict and avoid the cutting chatter. In addition, to control and suppress chatter vibrations, the use of piezoelectric actuators embedded in the tool holder [9], electromagnetic actuators integrated into the spindle system [10] and tunable clamping table [11] has been analyzed. In the milling process, the use of variable pitch cutters has demonstrated to improve productivity [12]. Different from the uniform pitch cutter, when a variable pitch cutter is used the dynamics model of cutting vibration changes from DDEs with a single delay to DDEs with multiple delays [13]. A common technique offline to predict unstable vibrations is the so-called stability lobes of the DDE based on Floquet theory [14], in which a curve describes the limit of stable vibration under feasible range values of cutting parameters.

The stability analysis of the milling process with multiple delays has been studied through different methods. Among all these methods, those with variable pitch tools play a critically important role [15]. Slavicek [16] was the one who first demonstrated the effectiveness of variable pitch cutters in suppressing vibrations in the milling process, he assumed a rectilinear tool motion for cutting teeth, and applied the theory of orthogonal stability to the irregular pitch of the tooth, by assuming an alternating step variation then, he obtained an expression of the stability limit as a function of the step angle variation. Budak [17,18] proposed an analytical method for nonconstant pitch milling cutters from a design point of view, showing for some applications how this variable effect helps to reduce self-excited vibrations, so he found that chatter stability can be improved significantly even at slow cutting speeds by properly designing the pitch angles. Altintas et al. [19] used the frequency domain method to analyze the milling stability of the variable pitch cutter and introduced a method to select the optimal pitch angles. Olgac and Sipahi proposed a mathematical approach, the cluster treatment of characteristic roots (CTCR), which optimizes the design of variable pitch cutters [20]. Jin et al. [21] presented an improved semi-discretization algorithm to predict the stability lobes for variable pitch cutters, which were verified and compared with previous works such as the Altintas analytical method (zero-order method) [19]. Comak and Budak [22] showed the optimal design of a tool for milling operations with variable geometry to widen the stability zones using the semi-discretization method, validating it experimentally. They also used a design methodology to determine the optimal pitch angle geometry for a given cutting condition, allowing increased stability.

Zatarain et al. [23] extended the multifrequency solution proposed by Budak and Altintas [24] to include the helix effect, they pointed out that the variation of the helix angle plays an important role in stability graphs due to repetitive vibrations driven by impact (flip), they found that the flip lobes became closed curves that are separated by horizontal lines where the depth of cut is equal to a multiple of the helix pitch. A similar phenomenon was confirmed using the semi-discretization method (SDM) in [25], meanwhile, B.R. Patel et al. [26] considered the influence of the helix angle of the tool to obtain an analytical force model, they found that isolated islands of instability can occur in the milling processes, which are induced by the helix angle of the tool and lead to separate regions of period-doubling and quasi-period behavior. Sims et al. [27] by using an adapted and time-averaged version of the SDM analyzed both the influence on the variation of the helix angle and the pitch angle of the tool to improve the prediction of vibrations and estimate predictions of surface errors. They used the semi-discretization method, the time-averaged semi-discretization method and the temporal finite element method to predict vibration stability for variable helix and variable pitch milling tools. Turner et al. [28] modeled and compared stability for variable pitch and helix angle cutters, demonstrating that variable helix angle tools can have higher stability and productivity.

Yusoff and Sims in [29] combined SDM with differential evolution to optimize variable helix milling tools to minimize vibration, their analysis predicted total vibration mitigation using the optimized variable helix milling tool at low radial immersion. Furthermore, Dombovari and Stepan [30] introduced a general mechanical model based on SDM to predict the linear stability of specialty cutters with optional continuous variation of the helix angle. Using an extended second-order SDM, Zhan et al. [15] predicted the stability lobe diagrams for tools with variable pitch angles. Meanwhile, Huang et al. [31] conducted a stability analysis for milling operations with variable pitch mills at variable speed, while Cai et al. [32] proposed an integrated process machine model based on the computer graphics method to simulate the milling process of a variable pitch cutter.

On the other hand, Olvera and Elías-Zuñiga in [33] led to the development of the enhanced multistage homotopy perturbation method (EMHPM) to solve differential delay equations (DDEs) with constant and variable coefficients and then this EMHPM was applied to predict the stability of a multivariate milling tool in which they consider the helix angle and the pitch angle variation of the cutting edges [34]. Based on the Laplace formulation, Sims [35] studied the stability of milling operations with a variable helix angle. Using the multi-frequency solution, Otto et al. [36] derived a dynamic process model where the non-linear shear force and the runout effect are included for milling with non-uniform pitch and variable helix tools. Niu et al. [37] found that runout can significantly increase the stability limits regardless of spindle speed ranges, while Olvera et al. [38] in a study for a thin-walled workpiece demonstrated that by considering the effects of the runout, the helix angle and characterization dependent on the cutting speed, more precise stability boundaries are achieved.

To demonstrate that one of the effective ways to suppress vibration in milling operations is to use tools with variable pitch and helix angle, Wang et al. [12] proposed an improved semi-discretization method based on Floquet- s theory. Since the delay between each cutting edge varies along with the axial depth of the tool in milling, they discretized the cutting tool in some axial layers to simplify the calculation. Iglesias et al. [39] presented a method to find the optimal angles between the inserts, and the stability diagrams were obtained through the iterative brute force (BF) method, which consists of an iterative maximization of stability through the semi-discretization method. They conclude that, if an optimal selection of the angle between the inserts is possible then, the material removal rate can be improved up to three times. Gou et al. [40] proposed an effective optimization method for the variable helical cutter introducing an index called "suppression factor" to measure stability quantitatively.

Therefore, in the present work, the EMHPM developed in [33] and extended for analysis of multivariable tools in [34], is now expanded to solve the dynamics of the machining process in milling in which the approximation to the delay is performed with polynomials of degree two and three. In order to study the proposed method performance in terms of convergency and computational cost, a multivariable milling tool with a variable pitch cutter and helix angle is used to determine milling process in stability domains.

This paper is summarized as follows. Section 2 focuses on the development of second- and third-order EMHPM for stability analysis of DDE. Section 3 studies the application of the second- and third-order EMHPM on the milling equation to demonstrate its improvement in the convergence rate. Section 4 is focused on the use of the third-order EMHPM to compute the stability analysis in milling for multivariable tools, and theoretical predictions with time-domain simulations are performed. Finally, some conclusions are drawn.

#### **2. Enhanced Multistage Homotopy Perturbation Method**

#### *2.1. Second-Order EMHPM*

Olvera et al. enhanced in [33] the multistage homotopy perturbation method (MHPM) proposed by Hashim [41]. The EMHPM considers the general case in which the nonlinear equation contains terms of the independent variable. This method is also useful to solve an n-dimensional DDE in the state-space form **.**

$$
\dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x} + \mathbf{B}(t)\mathbf{x}(t-\pi) \tag{1}
$$

where **A**(*t* + τ) = **A**(*t*), **B**(*t* + τ) = **B**(*t*), **x**(*t*), is the state vector, and τ is the time delay. Equation (1) can be written equivalently as: **.**

$$\dot{\mathbf{x}}\_{i}(T) - \mathbf{A}\_{l}\mathbf{x}\_{i}(T) \approx \mathbf{B}\_{l}\mathbf{x}\_{i}^{\top}(T) \tag{2}$$

where **x***i*(*T*) indicates the *m*-order solution for Equation (1) that satisfies the initial conditions **x***i*(0) = **x***i*−1, **A***<sup>t</sup>* and **B***<sup>t</sup>* are the periodic matrix whose values vary with time *t*. In [42], Puma et al. applied the first-order EMHPM to estimate the delayed term **x**<sup>τ</sup> *<sup>i</sup>* (*T*) in Equation (2), in which the period [*t*<sup>0</sup> − τ, *t*0] was discretized in *N* equally spaced discrete state values, and the function that describes the delayed term **x**<sup>τ</sup> *<sup>i</sup>* (*T*) in the delayed interval [*ti*−*N*, *ti*−*N*+1] was approximated as a first-order polynomial representation. Defining **x***<sup>i</sup>* ≡ **x***i*(*T*) to simplify the notation, Equation (2) can be written as

$$\dot{\mathbf{x}}\_{i}(T) = \mathbf{A}\_{t}\mathbf{x}\_{i}(T) + \mathbf{B}\_{t} \Big(\mathbf{x}\_{i-N} + \frac{N-1}{\tau}(\mathbf{x}\_{i-N+1} - \mathbf{x}\_{i-N})T\Big) \tag{3}$$

Figure 1a shows the representation of the approximation of the delayed term with the first-order polynomial. In the second-order EMHPM, to approximate the function that describes the delayed term **x**<sup>τ</sup> *<sup>i</sup>* (*T*) in Equation (2), the Lagrange equation is used, making use of the discrete values **x***i*−*N*, **x***i*−*N*+1, **x***i*−*N*+<sup>2</sup> as follows:

$$f\_n(\mathbf{x}) = \sum\_{i=0}^n Li(\mathbf{x}) f(\mathbf{x}\_i), \ L\_i(\mathbf{x}) = \prod\_{i=0, i \neq k}^n \frac{\mathbf{x} - \mathbf{x}\_i}{\mathbf{x}\_k - \mathbf{x}\_i} \tag{4}$$

to achieve a second-degree polynomial approximation, we have from Equation (4) that

$$P\_2(\mathbf{x}) = \frac{(\mathbf{x} - \Delta t)(\mathbf{x} - 2\Delta t)}{(0 - \Delta t)(0 - 2\Delta t)} f(\mathbf{x}\_{i-N}) + \frac{(\mathbf{x} - 0)(\mathbf{x} - 2\Delta t)}{(\Delta t - 0)(\Delta t - 2\Delta t)} f(\mathbf{x}\_{i-N+1}) + \frac{(\mathbf{x} - 0)(\mathbf{x} - \Delta t)}{(2\Delta t - 0)(2\Delta t - \Delta t)} f(\mathbf{x}\_{i-N+2}) \tag{5}$$

(**a**) (**b**)

**Figure 1.** (**a**) Scheme for the approximation of the delayed term by a first-order (solid black line), second-order (dashed blue line) and third-order (dotted red line) polynomial: (**b**) zoom in the time interval [*ti*−*N*+1, *ti*−*N*+3].

Substituting *x* = *T* and Δ*t* = (*N* − 1)/τ, we obtain the function that describes the delayed interval as:

$$\mathbf{x}\_{i-N+1}(T) \approx \mathbf{x}\_{i-N} + \left(\frac{N-1}{\tau}\right)T\left(-\frac{3}{2}\mathbf{x}\_{i-N} + 2\mathbf{x}\_{i-N+1} - \frac{1}{2}\mathbf{x}\_{i-N+2}\right) + \left(\frac{N-1}{\tau}\right)^2 \frac{T^2}{2} (\mathbf{x}\_{i-N} - 2\mathbf{x}\_{i-N+1} + \mathbf{x}\_{i-N+2}) \tag{6}$$

When the delay is approximated by a second-degree polynomial it is called second-order EMHPM and should not be confused with the order of solution *m* and which is determined by the last deformation taking into account the approximated solution. Notice that a polynomial of the second-degree requires three points. Likewise four points in the case of a third-degree polynomial as shown in Figure 1.

The procedure to calculate the second-order EMHPM solution is based on the EMHPM procedure described in [33]. The solution for second-order EMHPM is recursively expressed of **X***ik*(*T*) as

$$\mathbf{X}\_{\rm ik} = \mathbf{X}\_{\rm ik}^{a} + \mathbf{X}\_{\rm ik}^{b} + \mathbf{X}\_{\rm ik'}^{c}, k = 1, 2, 3 \dots \tag{7}$$

where

$$\mathbf{X}\_{i0}^{a} = \mathbf{x}\_{i-1} \mathbf{X}\_{i0}^{b} = \mathbf{X}\_{i0}^{c} = 0 \tag{8}$$

and

$$\begin{aligned} \mathbf{X}\_{\bar{k}}^{a} &= \frac{T}{k} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{a} + \mathbf{g}(k) \mathbf{B}\_{l} \mathbf{x}\_{i-N} \Big) \\ \mathbf{X}\_{\bar{k}}^{b} &= \frac{T}{k+1} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{b} + \mathbf{g}(k) \Big( \frac{N-1}{\tau} \Big) \mathbf{B}\_{l} T \Big( -\frac{3}{2} \mathbf{x}\_{i-N} + 2 \mathbf{x}\_{i-N+1} - \frac{1}{2} \mathbf{x}\_{i-N+2} \Big) \Big) \\ \mathbf{X}\_{\bar{k}}^{c} &= \frac{T}{k+2} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{c} + \mathbf{g}(k) \Big( \frac{N-1}{\tau} \Big)^{2} \mathbf{B}\_{l} \frac{T^{2}}{2} (\mathbf{x}\_{i-N} - 2 \mathbf{x}\_{i-N+1} + \mathbf{x}\_{i-N+2}) \Big) \end{aligned} \tag{9}$$

So, the solution of Equation (1) is obtained by adding each of the approximations **X***ik* of Equation (7).

$$\mathbf{x}\_{l}(T) \approx \sum\_{k=0}^{m} \mathbf{x}\_{ik}(T) \tag{10}$$

#### *2.2. Third-Order EMHPM Solution*

For the polynomial representation of the third-degree, the function that describes the delayed term **x**τ *<sup>i</sup>* (*T*) is approximated by a polynomial of order three, then Equation (4) of the Lagrange interpolator is used accordingly. In this case, it is necessary to employ the **x***i*−*N*, **x***i*−*N*+1, **x***i*−*N*+2, **x***i*−*N*+<sup>3</sup> discrete values. Following the same procedure described in Section 2.1, the function that describes the delayed interval is given as:

$$\begin{aligned} \mathbf{x}\_{i}^{\varepsilon}(T) &= \mathbf{x}\_{i-N+1}(T) \approx \mathbf{x}\_{i-N} + \left(\frac{\mathbf{x}\_{i-1}}{\mathfrak{r}}\right) \mathbf{T} \Big( -\frac{11}{6} \mathbf{x}\_{i-N} + 3 \mathbf{x}\_{i-N+1} - \frac{3}{2} \mathbf{x}\_{i-N+2} + \frac{1}{3} \mathbf{x}\_{i-N+3} \Big) + \\ & \quad \left(\frac{N-1}{\mathfrak{r}}\right)^{2} \frac{\mathbf{T}^{2}}{\mathfrak{T}} \Big( 2\mathbf{x}\_{i-N} - 5 \mathbf{x}\_{i-N+1} + 4 \mathbf{x}\_{i-N+2} - \mathbf{x}\_{i-N+3} \Big) + \left(\frac{N-1}{\mathfrak{r}}\right)^{3} \frac{\mathbf{T}^{3}}{\mathfrak{G}} \Big( -\mathbf{x}\_{i-N} + 3 \mathbf{x}\_{i-N+1} - 3 \mathbf{x}\_{i-N+2} + \mathbf{x}\_{i-N+3} \Big) \end{aligned} \tag{11}$$

Following the EMHPM procedure, the recursive solution of Equation (1) **X***ik*(*T*) is expressed as

$$\mathbf{X}\_{\overline{i}\overline{k}} = \mathbf{X}\_{\overline{i}\overline{k}}^{a} + \mathbf{X}\_{\overline{i}\overline{k}}^{b} + \mathbf{X}\_{\overline{i}\overline{k}}^{c} + \mathbf{X}\_{\overline{i}\overline{k}}^{d}, k = 1, 2, 3 \dots \tag{12}$$

where

$$\mathbf{X}\_{i0}^{a} = \mathbf{x}\_{i-1,\*} \mathbf{X}\_{i0}^{b} = \mathbf{X}\_{i0}^{c} = \mathbf{X}\_{i0}^{d} = 0 \tag{13}$$

and

$$\begin{aligned} \mathbf{X}\_{\mathrm{ik}}^{a} &= \frac{T}{k} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{a} + \mathbf{g}(k) \mathbf{B}\_{l} \mathbf{x}\_{i-N} \Big) \\ \mathbf{X}\_{\mathrm{ik}}^{b} &= \frac{T}{k+1} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{b} + \mathbf{g}(k) \Big( \frac{N-1}{\tau} \Big) \mathbf{B}\_{l} T \Big( -\frac{11}{6} \mathbf{x}\_{i-N} + \ 3 \mathbf{x}\_{i-N+1} - \frac{3}{2} \mathbf{x}\_{i-N+2} + \frac{1}{3} \mathbf{x}\_{i-N+3} \Big) \Big) \\ \mathbf{X}\_{\mathrm{ik}}^{c} &= \frac{T}{k+2} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{c} + \mathbf{g}(k) \Big( \frac{N-1}{\tau} \Big) \mathbf{B}\_{l} \frac{T^{2}}{2} \Big( 2 \mathbf{x}\_{i-N} - 5 \mathbf{x}\_{i-N+1} + 4 \mathbf{x}\_{i-N+2} - \mathbf{x}\_{i-N+3} \Big) \Big) \\ \mathbf{X}\_{\mathrm{ik}}^{d} &= \frac{T}{k+3} \Big( \mathbf{A}\_{l} \mathbf{X}\_{i(k-1)}^{d} + \mathbf{g}(k) \Big( \frac{N-1}{\tau} \Big)^{3} \mathbf{B}\_{l} \frac{T^{3}}{6} \Big( -\mathbf{x}\_{i-N} + 3 \mathbf{x}\_{i-N+1} - 3 \mathbf{x}\_{i-N+2} + \mathbf{x}\_{i-N+3} \Big) \Big) \end{aligned} \tag{14}$$

The approximate solution of Equation (1) can be obtained by substituting Equation (12) into Equation (10) adding each of the approximations **X***ik*.

#### *2.3. Stability Analysis*

To calculate the stability of the differential Equation (1) using the second-order EMHPM, the solution of Equation (10) for second-order EMHPM must be rewritten by grouping each of the discrete values **x**i, **x***i*−*N*+2, **x***i*−*N*+1, **x***i*−*N*, resulting in

$$\mathbf{x}\_{i}(T) \approx \mathbf{P}\_{i}(T)\mathbf{x}\_{i-1} + \mathbf{Q}\_{i}^{\prime}(T)\mathbf{x}\_{i-N+2} + \mathbf{Q}\_{i}(T)\mathbf{x}\_{i-N+1} + \mathbf{R}\_{i}(T)\mathbf{x}\_{i-N} \tag{15}$$

where

$$\begin{aligned} \mathbf{P}\_{i}(T) &= \sum\_{k=0}^{m} \frac{1}{k!} \mathbf{A}\_{i}^{k} T^{k}, \\ \mathbf{Q}\_{i}^{\prime}(T) &= \sum\_{k=1}^{m} \left( \frac{1}{(k+2)!} \binom{N-1}{\tau}^{2} \mathbf{A}\_{i}^{k-1} \mathbf{B}\_{l} T^{k+2} - \frac{1}{2(k+1)!} \binom{N-1}{\tau} \mathbf{A}\_{l}^{k-1} \mathbf{B}\_{l} T^{k+1} \right) \\ \mathbf{Q}\_{i}(T) &= \sum\_{k=1}^{m} \frac{1}{(k+1)!} \binom{N-1}{\tau} \mathbf{A}\_{l}^{k-1} \mathbf{B}\_{l} T^{k+1} - 2 \mathbf{Q}\_{i}^{\prime} i \\ \mathbf{R}\_{i}(T) &= \sum\_{k=1}^{m} \frac{1}{k!} \mathbf{A}\_{l}^{k-1} \mathbf{B}\_{l} T^{k} - \mathbf{Q}\_{i}^{\prime} i - \mathbf{Q}\_{i} \end{aligned} \tag{16}$$

Similarly, to compute the stability lobes for the third-order EMHPM, the solution of the differential Equation (1) for third-order EMHPM is rewritten as

$$\mathbf{x}\_{i}(T) \approx \mathbf{P}\_{i}(T)\mathbf{x}\_{i-1} + \mathbf{Q}\mathbf{y}\_{i}''(T)\mathbf{x}\_{i-N+3} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nkern-1.5mu\mathbf{y}\_{i} + \mathbf{Q}\nmid}$$

where

**<sup>P</sup>***i*(*T*) = *<sup>m</sup> k*=0 1 *k*! **A***k <sup>t</sup>Tk*, **Q"** *<sup>i</sup>*(*T*) = *<sup>m</sup> k*=1 <sup>1</sup> (*k*+1)! - *N*−1 τ **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>1</sup> - 1 3 <sup>−</sup> <sup>1</sup> (*k*+2)! - *N*−1 τ 2 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>2</sup> <sup>+</sup> <sup>1</sup> (*k*+3)! - *N*−1 τ 3 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>3</sup> **Q**- *<sup>i</sup>*(*T*) = *<sup>m</sup> k*=1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 (*k*+1)! - *N*−1 τ **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>1</sup> <sup>+</sup> <sup>1</sup> (*k*+2)! - *N*−1 τ 2 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>2</sup> - −7 2 + <sup>1</sup> (*k*+3)! - *N*−1 τ 3 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>3</sup> - 9 2 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ <sup>−</sup> <sup>15</sup> <sup>2</sup> **<sup>Q</sup>**- - *i* **<sup>Q</sup>***i*(*T*) = *<sup>m</sup> k*=1 <sup>1</sup> (*k*+1)! - *N*−1 τ **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tTk*+<sup>1</sup> <sup>−</sup> <sup>3</sup>**Q**- - *<sup>i</sup>* <sup>−</sup> <sup>2</sup>**Q**- *i* **<sup>R</sup>***i*(*T*) = *<sup>m</sup> k*=1 1 *k*! **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***tT<sup>k</sup>* <sup>−</sup> **<sup>Q</sup>**- - *<sup>i</sup>* <sup>−</sup> **<sup>Q</sup>**- *<sup>i</sup>* − **Q***<sup>i</sup>* (18)

The approximate solution obtained from Equation (17) was used to define a discrete map following the procedure described in [43]:

$$\mathbf{w}\_{l} = \mathbf{D}\_{l}\mathbf{w}\_{i-1} \tag{19}$$

where **w***i*−<sup>1</sup> is a vector with dimension equal to the total number of states (displacement and velocity) for all *N* discrete intervals:

$$\mathbf{w}\_{i-1} = \begin{bmatrix} \mathbf{x}\_{(i-1)}, \dot{\mathbf{x}}\_{(i-1)}, \mathbf{x}\_{(i-2)}, \dots, \mathbf{x}\_{(i-N)} \end{bmatrix}^T \tag{20}$$

**D***<sup>i</sup>* is a coefficient matrix and for the third-order EMHPM it has the form:

$$\mathbf{D}\_{i} = \begin{bmatrix} \mathbf{P} & 0 & 0 & 0 & \cdots & 0 & \mathbf{Q}^{\prime} \mathbf{i} & \mathbf{Q}^{\prime} \mathbf{i} & \mathbf{Q}\_{i} & \mathbf{R}\_{i} \\ \mathbf{I} & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbf{I} & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \mathbf{I} & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & \mathbf{I} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & \mathbf{I} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & \mathbf{I} & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \mathbf{I} & 0 \end{bmatrix} \tag{21}$$

It is important to point out that in the case of the second-order EMHPM, the matrix **D**<sup>i</sup> is like the matrix of the third-order EMHPM without the matrix **Q**-- *i* .

Then, the Floquet transition matrix **Φ** is calculated over the main period τ = (*N* − 1)/ Δ*t*, coupling each of the discrete maps **D***i*, *i* = 1, 2, ... , (*N* − 1), to obtain:

$$\Phi = \mathbf{D}\_{N-1} \mathbf{D}\_{N-2} \dots \mathbf{D}\_2 \mathbf{D}\_1 \tag{22}$$

Thus, the stability of Equation (1) is determined by calculating the eigenvalues of the transition matrix given by Equation (22). The eigenvalues of the transition matrix are actually the Floquet multipliers which are the exponents of each complex exponential functions that describe the motion of Equation (1). If the modulus of greatest magnitude is greater than or equal to one, it implies that the system will behave in an unstable way and the amplitude of the vibration will increase exponentially, otherwise it will have a stable behavior.

#### **3. Numeric Solution of the Milling Equation**

#### *3.1. Dynamic Model to the Milling Equation*

To validate the proposed EMHPM methods, the numerical solution of the delay differential equation analyzed by Olvera et al., in [33] was calculated, which describes the dynamic model of the milling process in one degree of freedom (DOF):

$$
\ddot{\mathbf{x}}(t) + 2\zeta a\_n \dot{\mathbf{x}}(t) + a\_n^2 \mathbf{x}(t) = -\frac{a\_p h\_{\mathbf{x}\mathbf{x}}(t)}{m\_m} (\mathbf{x}(t) - \mathbf{x}(t-\tau))\tag{23}
$$

where ζ is the modal damping ratio, ω*<sup>n</sup>* is the natural frequency of the workpiece, *ap* is the axial depth of cut, *mm* is the modal mass, τ represents the time delay corresponding to the hitting period between each tooth of the tool and *hxx*(*t*) is the specific cutting force in the *x*-direction due to flexibility in *x*-direction, which was calculated depending on the position of the tool

$$h\_{\rm xx}(t) = \sum\_{i z=1}^{z\_{\rm ref}} g(\phi\_{i\bar{z}}(t)) \sin \phi\_{i\bar{z}}(t) \left( K\_{\rm tc} \cos \phi\_{i\bar{z}}(t) + K\_{\rm nc} \sin \phi\_{i\bar{z}}(t) \right) \tag{24}$$

*zn* is the number of edges of the tool, *Ktc* and *Knc* are the average specific cut coefficients in the tangential and normal direction, respectively, and φ*iz*(*t*) is the angular position of each left edge described by

$$
\phi\_{\rm iz} = (2\pi n/60)t + 2\pi iz/z\_n \tag{25}
$$

where *n* is the spindle speed in revolution per minute (rpm). The function *g*(φ*iz*(*t*)) is a window function, which has the value of one when the current edge *iz* is cutting material, otherwise it takes the value zero.

In up-milling <sup>φ</sup>*st* = 0 and <sup>φ</sup>*ex* = cos−1(<sup>1</sup> <sup>−</sup> <sup>2</sup>*ad*), conversely in down-milling <sup>φ</sup>*st* = cos−1(2*ad* <sup>−</sup> <sup>1</sup>) and φ*ex* = π, *ad* is the radial immersion ratio of the cut and φ*st* and φ*ex* are the angular positions where each edge enters and leaves the workpiece.

The second- and third-order EMHPM is applied to obtain the solution of Equation (23) and it is compared with the solution given by the first-order EMHPM [33]; for a regular tool the matrix **A** and **B** are represented as:

$$\mathbf{A}\_{t} = \begin{bmatrix} 0 & 1 \\ -\omega\_{n}^{2} - \frac{a\_{p}l\_{\mathrm{tr}}(t)}{m\_{\mathrm{m}}} & -2\zeta\omega\_{n} \end{bmatrix}, \mathbf{B}\_{t} = \begin{bmatrix} 0 & 0 \\ \frac{a\_{p}l\_{\mathrm{tr}}(t)}{m\_{\mathrm{m}}} & 0 \end{bmatrix} \tag{26}$$

**A***<sup>t</sup>* and **B***<sup>t</sup>* correspond to the periodic matrix evaluated at time *t*. For demonstration purposes, time-domain simulations were computed for a full-immersion down-milling operation. We used the parameters employed by Insperger et al., in [43] where the stability lobes were also calculated. The modal parameters *fn* = 922 Hz, ω*<sup>n</sup>* = 5793 rad/s, ζ = 0.011 and *mm* = 0.03993 kg corresponds to a single degree of freedom. The tangential and normal cutting coefficients are *Ktc* = 6 <sup>×</sup> 108 N/m<sup>2</sup> and *Knc* = 2 <sup>×</sup> 108 N/m2 respectively for an end-mill with *zn* = 2. The time-domain solution was computed using the EMHPM considering *N* = 76 discrete intervals and *m* = 7. Two sets of cutting conditions were chosen for a fixed spindle speed value of *n* = 12,000 rpm where the axial depth of cut of *ap* = 1.5 mm corresponds to a stable cutting operation while that for an unstable operation *ap* = 3 mm was chosen. In Figure 2 we plot the second- and third-order EMHPM solutions and compare it with the first-order EMHPM and the dde23 routine in Matlab, which is used to integrate DDE.

**Figure 2.** Numerical comparison of the enhanced homotopy perturbation method (EMHPM) solutions of the milling equation, Equation (23), with the dde23 MATLAB routine. (**a**) Stable milling operation with *ap* = 1.5 mm, *ad* = 1 and *n* = 12000 rpm and (**b**) unstable milling operation with *ap* = 3 mm, *ad* = 1 and *n* = 12000 rpm.

#### *3.2. Numerical Comparison between Methods*

In order to observe the rate of convergence of the first-, second- and third-order EMHPM, we chose the stable case with cutting conditions *ap* = 1.5 mm, *ad* = 1 and *n* = 12,000 rpm presented in Figure 3a, and the unstable case with cutting conditions *ap* = 3 mm, *ad* = 1 and *n* = 12,000 rpm showed in Figure 3b. The rate of convergence was analyzed by computing the absolute error between the solution with *N* discrete intervals and a converged solution. All methods were compared against itself using the solution provided with *N* = 200 discrete intervals, which are considered the converged solution. In Figure 3a it is observed that the convergence is better for the second- and third-order than the first-order, however, the difference of convergence between second- and third-order with the parameters used was negligible. On the other hand, Figure 3b shows that for few discrete intervals the third-order EMHPM had the fastest convergence in comparison with the second- and the first-order EMHPM. However, the second-order and third-order curves behaved very similarly after *N* = 50 discrete intervals. It is important to mention that for a typical stability solution in the ranges of spindle speed 5000–10,000 rpm, *N* = 40 discrete intervals will be enough to have accurate predictions.

(**a**)

**Figure 3.** Convergence rate of absolute error between first-, second- and third-order EMHPM for down-milling operation. Cutting parameters for (**a**) *ap* = 1.5 mm, *ad* = 1 and *n* = 12,000 rpm and (**b**) *ap* = 3 mm, *ad* = 1 and *n* = 12,000 rpm.

Since the rate of convergence was proved for time-domain simulations, we next explored the convergence of the methods applied to the stability analysis. The stability lobes computed with the second- and third-order EMHPM for regular milling tools were compared with its predecessor for radial immersion value of *ad* =1 and the other parameters indicated above as it was used in [44]. Figure 4 shows the stability diagrams for spindle speed in the range 2000–3000 rev/min where the precision of the method was compromised due to the higher value of the time delay. While the shaded gray area represents the stability lobes computed with *N* = 200 discrete intervals in all subfigures, in each subfigure solid black lines draw the stability frontier for a specific discrete interval and using the first-, second- or third-order EMHPM. In Figure 4 the first, second and third column represents the solution for the first-, the second- and the third-order EMHPM respectively, while the first and the second row was for *N* = 60 and *N* = 100 discrete intervals, respectively. It is observed that the error achieved in the third-order EMHPM was less than those attained for the first-order and second-order EMHPM solutions. This confirms that the third-order EMHPM had the highest rate of convergence.

**Figure 4.** Stability diagrams for down-milling operation. The first (**a**–**c**) and second (**d**–**f**) rows of subfigures correspond to N = 60 and N = 100 discrete intervals, respectively. First (**a**,**d**), second (**b**,**e**) and third (**c**,**f**) columns correspond to the first-, second- and third-order EMHPM, respectively. Subfigure (**g**) shows a comparison between the third-order EMHPM (black line), the SDM (dot line) and the FDM (dash line) with N = 60 discrete intervals.

The results were also compared, in Figure 4g, with the semi-discretization method (SDM) presented by Insperger and Stépán in [43] (dot line) and with the full-discretization method (FDM) presented by Ding et al. in [45] (dash line). It is observed that the EMHPM converged faster than the SDM and FDM. Table 1 list the results for computation times for a different number of discretized intervals together with the absolute error between stability frontiers. Notice that the solution obtained with the EMHPM with *N* = 60 discrete intervals was faster than the SDM and the FDM and even the error was less in the solution by the EMHPM. For *N* = 100 discrete intervals the computation time for the FDM was similar to the solution obtained with the EMHPM but it was demonstrated that the solution by second- and third-order EMHPM requires a smaller number of discretized intervals to converge to the solution and using a smaller amount of computation time. Notice the error for the SDM was not calculated since there was no stability frontier in some values of spindle speed.


**Table 1.** Comparison of convergence for different methods for down-milling operation with *ad* = 1

It is noticeable that exists a significant improvement in the rate of convergence from first-order to second-order and third-order EMHPM, however, the difference between the second- and the third-order EMHPM is negligible if the number of discretized intervals increase. There is no best method between second- and third-order EMHPM in terms of rate of convergence and computation time since the precision depends on the nature of the studied problem. However, it is easy to prove that a higher-order approximation (fourth- and fifth-order EMHPM) could drastically increase the computational time without a significant improvement in the solution.

#### **4. Stability Analysis of Multivariable Milling Tools**

The EMHPM can be generalized for stability analysis of DDEs having multiple delays. A multivariable tool contains some of the following characteristics: uneven pitch between teeth, and/or at least one helix angle with a different value from the others. This analysis was developed by Compeán et al., in [34] by using the first-order EMHPM, where the methodology for the characterization of the cutting coefficients for a multivariable tool was discussed, and the dynamic behavior was studied from the productivity point of view. Since the angular spacing at the beginning of the edge is different between teeth (pitch) and the different values of helix angles of the edges between adjacent teeth, the angular spacing between teeth at a specific height changes continuously, which produces an infinite number of delays. A common approach to deal with the DDE with an infinite number of delay is to discretize the tool by cutting disks in the axial direction with a thickness Δ*adsk* to induce a DDE with a finite number of delays. A single disk still has the same number of flutes (discrete flutes) and considering that the maximum delay in the process is the period of rotation of the tool or the spindle rotation period τ*T*, then, it can be discretized in *N* − 1 intervals.

The angular position between two adjacent teeth in each cutting disk changes according to the axial position of the referred disk and is related to the expression ψ = *k*β*ap*, where *k*<sup>β</sup> = 2 tan β/2*D*. Here *D* is the diameter of the tool and ψ represents the cutting-edge offset angle due to the helix angle. A certain interval can be associated with a discrete time delay of each tooth *iz* and disk *l* using the following formulation

$$N\_{i\bar{z},l} = round\left( (N-1)\frac{\delta\phi\_{i\bar{z},l}}{2\pi} \right) \tag{27}$$

where δφ*iz*,*<sup>l</sup>* is the angular pitch between consecutive teeth for each disk, the round function converts the argument to the nearest integer. In Equation (27) *Niz*, is a table (matrix) of dimension *iz* × *l*. Since this procedure could generate several delayed terms and some of them with the same value of discrete time delay due to the discretization scheme, it is required to collect all the different (non-repeated) discrete time delays *dn* from *Niz*,*l*.

Thus, without loss of generality, the DDE with multiple delays can be written as

$$\dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t) + \sum\_{d=\min(d\_n)}^{\max(d\_n)} \mathbf{B}^d \mathbf{x}(t-\tau) \tag{28}$$

where **x** is the vector of states, **A**(*t* + τ*T*) = **A**(τ*T*), **B***d*(*t* + τ*T*) = **B***d*(τ*T*) and τ*<sup>T</sup>* is the period of rotation of the spindle. Following the EMHPM procedure, Equation (28) can be written equivalently by intervals as:

$$\dot{\mathbf{x}}\_{i}(T) - \mathbf{A}\_{t}\mathbf{x}\_{i}(T) \approx \sum\_{d=\min(d\_{n})}^{\max(d\_{n})} \mathbf{B}\_{t}^{d} \mathbf{x}\_{i}^{\pi\_{d}}(T) \tag{29}$$

being **x***i*(*T*) the solution by intervals of order m for Equation (28) in the *i* − *th* interval that satisfies the initial condition **<sup>x</sup>***i*(0) <sup>=</sup> <sup>x</sup>*i*−1, the matrices **At** and **<sup>B</sup>***<sup>d</sup> <sup>t</sup>* represent the values of the matrices **A**(*t*) and **B***d*(*t*) evaluated at time *t* respectively.

#### *4.1. Third-Order EMHPM for Multivariable Milling Tool*

To approximate the term associated with the delayed terms **x** τ*d <sup>i</sup>* (*T*) of Equation (29), the interval of the period τ*T*, [*t*<sup>0</sup> − τ*T*, *t*0] is discretized in *N* − 1 intervals that can be equal size. For simplicity, intervals of equal size Δ*t* = τ*T*/(*N* − 1) are chosen. Then it is assumed that the function **x** τ*d <sup>i</sup>* (*T*), which is defined in the interval [*ti*−*d*−1, *ti*−*d*+2], for the third-order EMHPM has the representation of the form:

$$\begin{aligned} \mathbf{x}\_{i}^{\tau\_{d}}(T) &= \mathbf{x}\_{i-d}(T) \approx \mathbf{x}\_{i-d-1} + \left(\frac{N-1}{\tau\_{T}}\right) \mathbf{T} \left(-\frac{11}{6} (\mathbf{x}\_{i-d-1}) + 3(\mathbf{x}\_{i-d}) - \frac{3}{2}(\mathbf{x}\_{i-d+1}) + \frac{1}{3}(\mathbf{x}\_{i-d+2})\right) + \\ &\left(\frac{N-1}{\tau\_{T}}\right)^{2} \frac{T^{2}}{2} (2(\mathbf{x}\_{i-d-1}) - 5(\mathbf{x}\_{i-d}) + 4(\mathbf{x}\_{i-d+1}) - (\mathbf{x}\_{i-d+2})) + \\ &\left(\frac{N-1}{\tau\_{T}}\right)^{3} \frac{T^{3}}{6} (-\mathbf{x}\_{i-d-1} + 3(\mathbf{x}\_{i-d}) - 3(\mathbf{x}\_{i-d+1}) + \mathbf{x}\_{i-d+2}) \end{aligned} \tag{30}$$

Defining **x***<sup>i</sup>* ≡ **x***i*(*Ti*) to simplify the notation, and substituting Equation (30) in Equation (29), the following equation is obtained:

$$\mathbf{x}\_{i}(T) = \mathbf{A}\mathbf{x}\_{i}(T) + \sum\_{d=\min\{d\_{i}\}}^{\max\{d\_{i}\}} \begin{pmatrix} \mathbf{B}\_{l}^{\ d}\mathbf{x}\_{i-d-1} + \mathbf{B}\_{l}^{\ d}T \Big(\frac{\mathbf{N}-1}{\tau\_{T}}\Big) \Big(-\frac{11}{6}\left(\mathbf{x}\_{i-d-1}\right) + 3\left(\mathbf{x}\_{i-d}\right) - \frac{3}{2}\left(\mathbf{x}\_{i-d+1}\right) + \frac{1}{3}\left(\mathbf{x}\_{i-d+2}\right)\Big) + \\\ \mathbf{B}\_{l}^{\ d}\frac{T^{2}}{2} \binom{\left(\mathbf{N}-1\right)^{2}}{\tau\_{T}} \Big(2\left(\mathbf{x}\_{i-d-1}\right) - 5\left(\mathbf{x}\_{i-d}\right) + 4\left(\mathbf{x}\_{i-d+1}\right) - \mathbf{x}\_{i-d+2}\Big) + \\\ \mathbf{B}\_{l}^{\ d}\frac{T^{3}}{6} \binom{\left(\mathbf{N}-1\right)^{3}}{\tau\_{T}} \Big(-\mathbf{x}\_{i-d-1} + 3\left(\mathbf{x}\_{i-d}\right) - 3\left(\mathbf{x}\_{i-d+1}\right) + \mathbf{x}\_{i-d+2}\Big) \end{pmatrix} \tag{31}$$

where

$$\mathbf{A}\_{l} = \begin{bmatrix} 0 & 1 \\ -\alpha\_{n}^{2} - \frac{\Delta t\_{\rm{disk}}}{m\_{n}} \sum\_{d=\min(d\_{n})}^{\max(d\_{n})} h\_{yy}^{d} & -2\zeta a \nu\_{n} \\\end{bmatrix} \cdot \mathbf{B}\_{l}^{d} = \begin{bmatrix} 0 & 0 \\\frac{\Delta t\_{\rm{disk}}}{m\_{n}} \sum\_{d=\min(d\_{n})}^{\max(d\_{n})} h\_{yy}^{d} & 0 \\\end{bmatrix} \tag{32}$$

here, *hyy* is the specific cutting force in the *y*-direction due to flexibility in *y*-direction, which is used for thin wall machining. This force was calculated depending on the position of the tool via the following equation:

$$h\_{yy}(t) = \sum\_{iz=1}^{z\_{\text{fl}}} \mathcal{g}(\phi\_{i\bar{z}}(t)) \cos \phi\_{i\bar{z}}(t) (-K\_{\text{fc}} \sin \phi\_{i\bar{z}}(t) + K\_{\text{nc}} \cos \phi\_{i\bar{z}}(t)) \tag{33}$$

#### Then, solving Equation (31) yields

**X***i*<sup>0</sup> = **x***i*−<sup>1</sup> **X***i*<sup>1</sup> = *T***A***t*(**x***i*−1) + max (*dn*) *d*=min(*dn*) ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *TBt <sup>d</sup>*(**x***i*−*d*−1) <sup>+</sup> *Bt d T*<sup>2</sup> 2 - *N*−1 τ*T* -<sup>−</sup><sup>11</sup> <sup>6</sup> (**x***i*−*d*−1) <sup>+</sup> <sup>3</sup>(**x***i*−*d*) <sup>−</sup><sup>3</sup> <sup>2</sup> (**x***i*−*d*+1) <sup>+</sup> <sup>1</sup> <sup>3</sup> (**x***i*−*d*+2) +*Bt d T*<sup>3</sup> 6 - *N*−1 τ*T* 2 (2(**x***i*−*d*−1) − 5(**x***i*−*d*) + 4(**x***i*−*d*+1) − **x***i*−*d*+2) +*Bt d T*<sup>4</sup> 24 - *N*−1 τ*T* 3 (−**x***i*−*d*−<sup>1</sup> + 3(**x***i*−*d*) − 3(**x***i*−*d*+1) + **x***i*−*d*+2) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ **<sup>X</sup>***i*<sup>2</sup> <sup>=</sup> *<sup>T</sup>*<sup>2</sup> <sup>2</sup> **A***<sup>t</sup>* <sup>2</sup>(**x***i*−1) + **<sup>A</sup>***<sup>t</sup>* max (*dn*) *d*=min(*dn*) ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *T*2 <sup>2</sup> *Bt <sup>d</sup>*(**x***i*−*d*−1) <sup>+</sup> *Bt d T*<sup>3</sup> 6 - *N*−1 τ*T* -<sup>−</sup><sup>11</sup> <sup>6</sup> (**x***i*−*d*−1) <sup>+</sup> <sup>3</sup>(**x***i*−*d*) <sup>−</sup><sup>3</sup> <sup>2</sup> (**x***i*−*d*+1) <sup>+</sup> <sup>1</sup> <sup>3</sup> (**x***i*−*d*+2) + *Bt d T*<sup>4</sup> 24 - *N*−1 τ*T* 2 (2(**x***i*−*d*−1) − 5(**x***i*−*d*) + 4(**x***i*−*d*+1) − **x***i*−*d*+2)+ *Bt d T*<sup>5</sup> 120 - *N*−1 τ*T* 3 (−**x***i*−*d*−<sup>1</sup> + 3(**x***i*−*d*) − 3(**x***i*−*d*+1) + **x***i*−*d*+2) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ . . . (34)

$$\mathbf{X}\_{k} = \frac{1}{k!} \mathbf{T}^{k} \mathbf{A}\_{i}^{k} (\mathbf{x}\_{i-1}) + \mathbf{A}\_{i}^{k-1} \sum\_{d=\min(d,j)}^{\max(d,j)} \begin{pmatrix} \frac{1}{k!} \mathbf{T}^{k} \mathbf{B}\_{i}^{d} (\mathbf{x}\_{i-d-1}) + \frac{1}{(k+1)!} \mathbf{T}^{k+1} \mathbf{B}\_{i}^{d} \left(\frac{\underline{N-1}}{\underline{r}\_{I}}\right) \left(-\frac{11}{6} (\mathbf{x}\_{i-d-1}) + 3(\mathbf{x}\_{i-d}) - \frac{3}{2} (\mathbf{x}\_{i-d+1}) + \frac{1}{(k+1)!} \mathbf{B}\_{i}^{d} (\mathbf{x}\_{i-1})\right) \\\\ \frac{1}{k!} (\mathbf{x}\_{i-d+2}) + \frac{1}{(k+2)!} \mathbf{T}^{k+2} \mathbf{B}\_{i}^{d} \left(\frac{\underline{N-1}}{\underline{r}\_{I}}\right)^{2} (2(\mathbf{x}\_{i-d-1}) - 5(\mathbf{x}\_{i-d}) + 4(\mathbf{x}\_{i-d+1}) - \mathbf{x}\_{i-d+2}) + \\ \frac{1}{(k+3)!} \mathbf{T}^{k+3} \mathbf{B}\_{i}^{d} \left(\frac{\underline{N-1}}{\underline{r}\_{I}}\right)^{3} (-\mathbf{x}\_{i-d-1} + 3(\mathbf{x}\_{i-d}) - 3(\mathbf{x}\_{i-d+1}) + \mathbf{x}\_{i-d+2}) \end{pmatrix}$$

Notice that Equation (34) can be written recursively as

$$\mathbf{X}\_{\vec{\text{ik}}} = \mathbf{X}\_{\vec{\text{ik}}}^{a} + \mathbf{X}\_{\vec{\text{ik}}}^{b} + \mathbf{X}\_{\vec{\text{ik}}}^{c} + \mathbf{X}\_{\vec{\text{ik}}}^{d}, k = 1, 2, 3 \dots \tag{35}$$

where **X***<sup>a</sup> <sup>i</sup>*<sup>0</sup> <sup>=</sup> **<sup>x</sup>***i*−1, **<sup>X</sup>***<sup>b</sup> <sup>i</sup>*0, <sup>=</sup> **<sup>X</sup>***<sup>c</sup> <sup>i</sup>*<sup>0</sup> <sup>=</sup> **<sup>X</sup>***<sup>d</sup> <sup>i</sup>*<sup>0</sup> = 0 and

**X***a ik* <sup>=</sup> *<sup>T</sup> k* ⎛ ⎜⎜⎜⎜⎝ **A***t***X***<sup>a</sup> <sup>i</sup>*(*k*−1) <sup>+</sup> *<sup>g</sup>*(*k*) max (*dn*) *d*=min(*dn*) **B***d <sup>t</sup>* **x***i*−*d*−<sup>1</sup> ⎞ ⎟⎟⎟⎟⎠ **X***b ik* <sup>=</sup> *<sup>T</sup> k*+1 ⎛ ⎜⎜⎜⎜⎝ **A***t***X***<sup>b</sup> <sup>i</sup>*(*k*−1) <sup>+</sup> *<sup>g</sup>*(*k*) max (*dn*) *d*=min(*dn*) - *N*−1 τ **B***d t T* - <sup>−</sup><sup>11</sup> <sup>6</sup> (**x***i*−*d*−1) <sup>+</sup> <sup>3</sup>(**x***i*−*d*) <sup>−</sup><sup>3</sup> <sup>2</sup> (**x***i*−*d*+1) <sup>+</sup> <sup>1</sup> <sup>3</sup> (**x***i*−*d*+2) ⎞ ⎟⎟⎟⎟⎠ **X***c ik* <sup>=</sup> *<sup>T</sup> k*+2 ⎛ ⎜⎜⎜⎜⎝ **A***t***X***<sup>c</sup> <sup>i</sup>*(*k*−1) <sup>+</sup> *<sup>g</sup>*(*k*) max (*dn*) *d*=min(*dn*) - *N*−1 τ 2 **B***d t T*2 <sup>2</sup> (2(**x***i*−*d*−1) − 5(**x***i*−*d*) + 4(**x***i*−*d*+1) − **x***i*−*d*+2) ⎞ ⎟⎟⎟⎟⎠ **X***d ik* <sup>=</sup> *<sup>T</sup> k*+3 ⎛ ⎜⎜⎜⎜⎝ **A***t***X***<sup>d</sup> <sup>i</sup>*(*k*−1) <sup>+</sup> *<sup>g</sup>*(*k*) max (*dn*) *d*=min(*dn*) - *N*−1 τ 3 **B***d t T*3 <sup>6</sup> (−**x***i*−*d*−<sup>1</sup> + 3(**x***i*−*d*) − 3(**x***i*−*d*+1) + **x***i*−*d*+2) ⎞ ⎟⎟⎟⎟⎠ (36)

the solution of order *m* for Equation (31) was obtained by adding each of the approximations *k* of Equation (35). Similar to Equation (15), to obtain the stability graphs the solution of Equation (35) is rewritten by grouping the discrete states, which results in:

$$\mathbf{x}\_{i}(T) \approx \mathbf{P}\_{i}(T)\mathbf{x}\_{(i-1)} + \sum\_{d=\min(d\_{n})}^{\max(d\_{n})} \left( \mathbf{Q}\_{i}^{"d}(T)\mathbf{x}\_{i-d+2} + \mathbf{Q}\_{i}^{"d}(T)\mathbf{x}\_{i-d+1} + \mathbf{Q}\_{i}^{d}(T)\mathbf{x}\_{i-d} + \mathbf{R}\_{i}^{d}(T)\mathbf{x}\_{i-d-1} \right) \tag{37}$$

where

**<sup>P</sup>***i*(*T*) = *<sup>m</sup> k*=0 1 *k*! **A***k <sup>t</sup>Tk*, **Q "***<sup>d</sup> <sup>i</sup>* (*T*) = *<sup>m</sup> k*=1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 (*k*+1)! - *N*−1 τ*T* **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>1</sup> - 1 3 <sup>−</sup> <sup>1</sup> (*k*+2)! - *N*−1 τ*T* 2 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+2<sup>+</sup> 1 (*k*+3)! - *N*−1 τ*T* 3 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>3</sup> ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ **Q** - *d <sup>i</sup>* (*T*) = *<sup>m</sup> k*=1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 (*k*+1)! - *N*−1 τ*T* **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>1</sup> <sup>+</sup> <sup>1</sup> (*k*+2)! - *N*−1 τ*T* 2 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>2</sup> - −7 2 + 1 (*k*+3)! - *N*−1 τ*T* 3 **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>3</sup> - 9 2 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ <sup>−</sup> <sup>15</sup> <sup>2</sup> **Q"***<sup>d</sup> i* **Q***d <sup>i</sup>* (*T*) = *<sup>m</sup> k*=1 <sup>1</sup> (*k*+1)! - *N*−1 τ*T* **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup>k*+<sup>1</sup> <sup>−</sup> <sup>3</sup>**Q"***<sup>d</sup> <sup>i</sup>* − 2**Q** - *d i* **R***d <sup>i</sup>* (*T*) = *<sup>m</sup> k*=1 1 *k*! **<sup>A</sup>***k*−<sup>1</sup> *<sup>t</sup>* **<sup>B</sup>***<sup>d</sup> <sup>t</sup> <sup>T</sup><sup>k</sup>* <sup>−</sup> **<sup>Q</sup>"***<sup>d</sup> <sup>i</sup>* − **Q** - *d <sup>i</sup>* <sup>−</sup> **<sup>Q</sup>***<sup>d</sup> i* (38)

The approximate solution obtained from Equation (37) was used to define a discrete map:

$$\mathbf{w}\_{i} = \mathbf{D}\_{i}\mathbf{w}\_{i-1} \tag{39}$$

where **w***i*−<sup>1</sup> is a vector represented by Equation (20) and **D***<sup>i</sup>* is a coefficient matrix given by


The transition matrix **Φ** over the period τ*<sup>T</sup>* = (*N* − 1)/ Δ*t* was determined by coupling each solution **x***<sup>i</sup>* through the discrete map *Di*, *i* = 1, 2, ... , (*N* − 1). However, the computational cost can be reduced by computing only the transition matrix up to the maximum delayed term without losing precision in the calculation of the eigenvalues:

$$\Phi = \mathbf{D}\_{\text{Nmax}} \mathbf{D}\_{\text{Nmax}-1} \dots \mathbf{D}\_2 \mathbf{D}\_1 \tag{41}$$

Thus, the stability graphs of Equation (28) were determined by computing the eigenvalues of the transition matrix of Equation (41). The results obtained from the EMHPM were corroborated with the stability lobes in the study of multivariate tools [27].

#### *4.2. Experimental Characterization of One Degree of Freedom Milling Equation and Cutting Force Model*

#### 4.2.1. Experimental Modal Analysis

An experimental workpiece was assembled with a 7075T6 aluminum block of 101 mm × 172 mm supported by two thin plates (walls) with a thickness of 4.5 mm. This assembly mimics a DOF as described in Equation (23). The workpiece assembly was rigidly fixed to the workbench of a Makino F3 machining center. For modal analysis, tap testing was performed using a 352C68 PCB Piezotronics accelerometer and an impact hammer model 9722A500. The signals were acquired with a Polytec VIB-E-220 data acquisition card and processed with VibSoft signal analyzer software as shown in Figure 5a. Using the CutPro 8 software, the modal parameters were fitted resulting the values ζ = 0.068, *mm* = 3.8 kg, *fm* = 132 Hz and ω*<sup>n</sup>* = 829 rad/s.

**Figure 5.** Scheme of the experimental setup for (**a**) the modal analysis and (**b**) cutting forces characterization.

#### 4.2.2. Experimental Determination of Cutting Coefficients

The force model in Equation (33) was used to predict the cutting force magnitude for a given depth of cut. It is based on a mechanistic approach that assumes a relationship between forces and the uncut chip thickness by means of the cutting coefficients. The cutting force model was established by introducing cutting (shearing) and edge coefficients for the tangential and normal directions of the milling tool. The characterization procedure assumed the linear relationship between the averaged experimental cutting forces *F* and the feed rate *fz* in *x*- and *y*- directions. This relationship is established as follows:

$$
\overline{F} = f\_{\overline{z}} \overline{F}\_{\overline{c}} + \overline{F}\_{\overline{c}} \tag{42}
$$

Here, *F <sup>c</sup>* and *F <sup>e</sup>* are the cutting shear and edge components, respectively. The experimental forces at each feed rate are measured, and the cutting-edge components *F <sup>c</sup>* and *F <sup>e</sup>* were evaluated

$$K\_{tc} = 4 \frac{F\_{yc}}{z\_{n}a\_{p}}, \; K\_{nc} = -4 \frac{\overline{F}\_{xc}}{z\_{n}a\_{p}} \tag{43}$$

A multivariable cutter provided by a local toolmaker was characterized by using Equation (43) and the experimental setup shown in Figure 5b. Table 2 summarizes the main geometric characteristics of the multivariable tool. A total of five cuttings were performed for full radial immersion in aluminum 7075T6 during dry machining. The forces were recorded by using a dynamometer 9257B Kistler and the spindle speed was set at 3000 rpm based on the dynamometer's natural frequency to avoid the amplification of milling forces. The force signals were acquired using a VibSoft-20 acquisition card at a sample rate of 48 kHz and processed in a custom-made MATLAB app to remove drift and noise. Cutting forces data were collected for the axial depth of cut of 2 mm and four values of feed per

tooth 0.05, 0.10, 0.015 and 0.20, so the resulting cutting coefficients *Ktc* for the tooth 1, 2, 3 and 4 were 1215 <sup>×</sup> 106, 1369 <sup>×</sup> 106, 897 <sup>×</sup> 10<sup>6</sup> and 1799 <sup>×</sup> 106 N/m<sup>2</sup> respectively, while that the coefficients *Knc* for the tooth 1, 2, 3 and 4 resulted 272 <sup>×</sup> 106, 520 <sup>×</sup> <sup>10</sup>6, 801 <sup>×</sup> <sup>10</sup><sup>6</sup> and 859 <sup>×</sup> 106 <sup>N</sup>/m2 respectively.


**Table 2.** Main geometric parameters of multivariable tool.

#### *4.3. Stability Analysis of 1 DOF Milling with a Multivariable Tool*

The stability lobes computed for the multivariable tool using the third-order EMHPM with a mesh of 400 × 200 (*n* × *ap*) are shown in Figure 6 together with stability lobes for a regular tool (angles of 90◦ and helix angles of 30◦ for all flutes). An approximation of order *m* = 7 was used with *N* = 241 and *ad* = 1 mm. Notice from Figure 6 that the stable zone obtained for the multivariable tool was significantly larger, meaning that the critical depth of cut was higher in most spindle speeds, which allowed having more global productivity. It is also observed in the range of spindle speed between 2000 and 3000 rpm, a stable peninsula formed with axial depth ranging from 11 to 20 mm or higher values of critical depth of cut. For instance, for the multivariable cutter at 2500 rpm, the critical depth of cut *ap* was 2.17 mm, however it became stable again as shown in Figure 6 for the interval values between 11 and 20 mm. To validate this unexpected behavior, we performed several time-domain simulations using the third-order EMHPM solution described by Equation (37).

**Figure 6.** (**a**) Comparison of stability lobes for regular (black solid line) and multivariable (red solid line) cutters by using the third-order EMHPM and (**b**) zoom in on chosen cutting conditions for time-domain simulations. The selected points are marked as follows: unstable (cross mark), stable (circle mark) and transition (plus mark) cutting conditions.

Furthermore, the simulated vibrations for the chosen cutting conditions were analyzed using the continuous wavelet transform (CWT), the power spectral density (PSD) and Poincaré maps (PM). The CWT is a time-frequency representation of a signal that offers the capability to observe how frequencies evolve in time. The scalograms display the absolute value of CWT of the simulated vibration and therefore, they were used to detect chatter phenomena that appeared when milling with a multivariable tool. The PSD is based on the Fourier transform that provides the transformation from the time-domain to the frequency-domain. Additionally, PSD is defined as the squared value of the signal and describes the power of a signal or time series distributed over different frequencies [46]. Moreover, a PM represents points in phase space, which are sampled every spindle rotation [47]. The frequencies *f* of the CWT and PM were normalized *fn* = *f*/ *fh* according to the spindle frequency *fh*. When milling with a regular milling tool the excitation frequency *fe* is equal to *zn* times frequencies of the spindle speed *fh* but in a multivariable tool, there are several excitation frequencies since the angular spacing between teeth change as a function of the axial depth of cut.

Figure 7 illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions A, B and C for the axial depths of cut of 1.0, 1.7 and 1.7 mm respectively. Figure 7a–c refers to the vibrations of the cutting conditions A marked in Figure 6, using a regular tool. The scalogram in Figure 7a identifies point A as a stable cutting since normalized cutting frequencies present a dominant value of *fn* = 3.2, which corresponds to the natural frequency *fm* = 132 Hz. This is also confirmed by the PSD analysis shown in Figure 7b. The PM illustrated in Figure 7c shows a vibration that decreased with time and sampled data concentrated in the center confirmed a typical stable case. When the axial depth of cut was increased to 1.7 mm, the stability diagram predicted unstable cutting conditions according to the stability lobes for the regular tool. This case is denoted with cutting conditions B and the corresponding scalogram (shown in Figure 7d) illustrated how the intensity of the dominant frequency increased with time even when the excitation frequency was the same as the case in A.

**Figure 7.** Analysis of cutting conditions A, B and C. Continuous wavelet transform (CWT) scalograms: (**a**,**d**,**g**); power spectral density (PSD): (**b**,**e**,**h**) and Poincaré maps (PM): (**c**,**f**,**i**) corresponds to the cutting conditions A, B and C respectively.

The PM diagram shown in Figure 7f exhibited a vibration far from zero. In fact, the PM diagram shows that the vibration amplitude grows exponentially because our equation of motion did not consider nonlinear effects such as those that appeared when the tool lost contact with the workpiece. Both cutting conditions A and B agreed with the stability boundaries in Figure 6. Now, the cutting conditions B were used but with a multivariable tool, which was referred to as cutting conditions C. The CWT plotted in Figure 7g described completely different results since there were no single dominant frequencies in comparison with cutting conditions A, but appeared several frequencies around *fn* = 3.2 and close to *fn* = 1 that reduced in intensity with time, suggesting a stable cutting. Figure 7i illustrates how the vibration amplitude approached to zero when using a multivariable tool in contrast to the PM obtained for the regular tool and exhibited in Figure 7f. This can be explained by observing that there were several excitation frequencies due to the irregular pitch and helix angles that break a single excitation frequency avoiding regenerative chatter phenomena.

Figure 8 illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions D, E, F and G for the axial depths of cut of 2.3, 3.0, 8.55 and 18 mm respectively. Notice that a stable case C was already validated when the axial depth was 1.7 mm in Figure 7g–i that corresponded to cutting conditions under the stability boundaries shown in Figure 6. For case D, a transient cutting condition was chosen very close to the critical axial depth of the cut. It is interesting to point out that transition cutting conditions in the CWT scalogram shown in Figure 8a not only shows frequencies with higher intensity in comparison with the stable case B, but also presents shifted frequencies that varied in intensity every single revolution. This shifting suggests a marginally stable cutting condition that was confirmed by the PM illustrated in Figure 8c, where circular trajectories were described close to the center point.

**Figure 8.** *Cont.*

**Figure 8.** Analysis of cutting conditions D, E, F and G. CWT scalogram:s (**a**,**d**,**g**,**j**); PSD: (**b**,**e**,**h**,**k**) and PM: (**c**,**f**,**i**,**l**) corresponds to the cutting conditions D, E, F and G respectively.

Unstable vibrations that appeared for case E were because of the intensity of frequencies increased exponentially with time, see Figure 8d. Notice that other frequencies arose with time close to the values of *fn* = 0.5 and *fn* = 1.5. These frenquencies also occurred for cutting conditions D, which is an indication of the appearance of chatter phenomena. In contrast to Figure 7i for a stable case, Figure 8f exhibited few trajectories because the vibration amplitude was out of the range selected (±1 mm). The qualitative and quantitative dynamic behaviour due to cutting conditions F, and illustrated in Figure 8g-i, were classified as transition cutting behaviour. Here, a more severe shifting in frequencies was observed in the scalogram (Figure 8g). From Figure 8g, it is seen that drastic shifting occurred in the time domain in the range of normalized frequencies from 3.5 to 6. It was also evident in the PM showed in Figure 8i, that the amplitude of vibration remained below 1 mm during several revolutions of the tool but the amplitude of vibration never aproached to the center point, in contrast to the stable cutting condition C shown in Figure 7i in which the oscillation aplitudes aproached the center.

An interesting dynamic behaviour was observed in the milling cutting process when the cutting conditions were selected in the middle of the stable peninsula, above unstable cutting conditions such as E cutting conditions. The axial depth of the cut was increased from the unstable axial depth of cut of 3–18 mm, 6 times higher of the stable cutting condition C and 2 times higher than the unstable condition E. Since the vibration quickly decreased in a few revolutions no dominant frequencies appeared in the CWT and PSD failed to clearly identify a dominant frequency since the vibration amplitude decreased to zero after few revolutions, as confirmed by the PM shown in Figure 8l.

Figure 9 shows the normalized excitation frequencies that the multivariable tool produced for a fixed spindle speed of 2500 rpm. The total number of disks of 50 μm of thickness was grouped in sets of each millimeter in the axial direction. The waterfall plot in Figure 9 explains that a stable peninsula was formed above 11 mm because the workpiece was excited with several frequencies simultaneously. For instance, for a milling operation with the axial depth of cut of 1 mm (stable cutting), 80 discrete disks were cut with four normalized excitation frequencies values (3.3, 3.6, 4.5 and 5.1). On the other hand, when milling at 18 mm (stable cutting), there were 14 normalized excitation frequencies (3.30,

3.35, 3.39, 3.44, 3.49, 3.54, 3.60, 4.55, 4.64, 4.73, 4.82, 4.92, 5.02 and 5.13), most of them with at least 115 discrete disks.

**Figure 9.** The number of discrete disks and discrete excitation frequencies as a function of the axial depth of cut for the multivariable tool.

#### **5. Conclusions**

In this work, quadratic and cubic polynomials were used to approximate the delayed terms of delay differential equations. Numerical simulations showed that using second- and third-order EMHPM improved the convergence rate and required less computational time when compared to the first-order EMHPM, and to semi-discretization and full-discretization methods, since fewer approximations or less discrete intervals were needed to reduce the computation time.

To further assess the applicability of the proposed method, the third-order EMHPM was used for determining the stability bounds in one-degree-of-freedom milling operation with a multivariable tool, demonstrating that the stability zone improved in comparison with a regular tool. For instance, at 2500 rpm the critical axial depth of the cut was 1.3 mm using the regular milling tool. However, using the multivariable tool, the critical axial depth of the cut was increased until 2.17 mm but more interesting, a stable zone appeared above 8.55 mm.

The CWT scalograms, PSD charts and PM were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool. Numerical solutions confirmed the system dynamics behavior predicted by the third-order EMHPM.

Based on the above results, this paper provided evidence the third-order EMHPM could be used to study dynamic phenomena that appeared at higher axial depths of cut due to the multivariable design of the tool, which broke the excitation frequencies at a lower depth of cut.

**Author Contributions:** Conceptualization, J.d.l.L.S. and D.O.-T.; Methodology, J.d.l.L.S.; Resources, O.M.-R., A.E.-Z. and L.N.L.d.L.; Supervision, D.O.-T., A.E.-Z. and L.N.L.d.L.; Validation, J.d.l.L.S., D.O.-T. and G.U.; Writing—original draft, J.d.l.L.S. and D.O.-T.; Writing—review and editing, O.M.-R., G.U. and A.E.-Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Tecnológico de Monterrey through the Research Group of Nanotechnology for Devices Design, and by the Consejo Nacional de Ciencia y Tecnología de México (Conacyt), Project Numbers 242269, 255837, 296176, and National Lab in Additive Manufacturing, 3D Digitizing and Computed Tomography (MADiT) LN299129.

**Acknowledgments:** The authors acknowledge the technical assistance from Oscar Escalera at Tecnologico de Monterrey.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
