*Article* **Implicit Subspace Iteration to Improve the Stability Analysis in Grinding Processes**

**Jorge Alvarez 1,\*, Mikel Zatarain 1, David Barrenetxea 1, Jose Ignacio Marquinez <sup>1</sup> and Borja Izquierdo <sup>2</sup>**


Received: 9 October 2020; Accepted: 18 November 2020; Published: 19 November 2020

#### **Featured Application: The calculation time reduction of the new methodology for obtaining the stability maps in grinding processes contributes to the industrialization of chatter avoidance technologies.**

**Abstract:** An alternative method is devised for calculating dynamic stability maps in cylindrical and centerless infeed grinding processes. The method is based on the application of the Floquet theorem by repeated time integrations. Without the need of building the transition matrix, this is the most efficient calculation in terms of computation effort compared to previously presented time-domain stability analysis methods (semi-discretization or time-domain simulations). In the analyzed cases, subspace iteration has been up to 130 times faster. One of the advantages of these time-domain methods to the detriment of frequency domain ones is that they can analyze the stability of regenerative chatter with the application of variable workpiece speed, a well-known technique to avoid chatter vibrations in grinding processes so the optimal combination of amplitude and frequency can be selected. Subspace iteration methods also deal with this analysis, providing an efficient solution between 27 and 47 times faster than the abovementioned methods. Validation of this method has been carried out by comparing its accuracy with previous published methods such as semi-discretization, frequency and time-domain simulations, obtaining good correlation in the results of the dynamic stability maps and the instability reduction ratio maps due to the application of variable speed.

**Keywords:** stability; machining; chatter; grinding

#### **1. Introduction**

Grinding is a finishing process characterized by its ability to achieve workpieces with smooth surfaces and fine tolerances. One of the most important limitations in grinding processes is the self-excited vibration or chatter. The appearance of this regenerative effect leads to lower form accuracy and worse surface finish of the ground parts, loss of productivity, grinding machine damage or even wheel breakage.

Chatter in grinding occurs when the grinding forces change during the process due to the variations of the chip thickness in the presence of vibrations, then the corresponding delayed responses to the forces increase the undulations over the workpiece surface. The generative effect can appear on both the workpiece and the grinding wheel, adding complexity to the dynamic behavior. The undulations generated on the workpiece surface grow quite rapidly, whereas those generated in the grinding wheel surface grow more gradually [1]. Besides, when the vibration amplitude builds up to a certain limit, those undulations on the grinding wheel surface are removed through the dressing process. For these reasons, this paper is only focused on the workpiece regenerative chatter.

Different approaches have been developed to analyze this phenomenon and devise practical methods for suppressing these vibrations [2]. The first discussion of regenerative stability of a grinding process was done by Hahn using the Nyquist criteria [3]. Then, Snoeys and Brown [4] performed a similar analysis by means of a doubly regenerative closed-loop block diagram and obtained the characteristic equation of the process and a stability criterion. Thompson proposed an alternative approach obtaining solutions of double regenerative chatter for the stability boundary [5] and chatter growth [6]. Later, different time-domain simulation approaches were applied to grinding stability analysis and regenerative chatter prediction [7,8], considering also non-linear effects [9,10].

Alvarez et al. applied the semi-discretization method proposed by Insperger et al. for milling [11] and turning [12]—later used for selection of variable pitch for chatter suppression in face milling [13] to infeed cylindrical grinding process [14] and traverse cylindrical grinding process [15]—demonstrating the potential of this approach to analyze the application of variable speed to chatter avoidance with an efficient computing effort comparing to previous time-domain simulations. This periodically variable speed is one of the most extended techniques for chatter suppression in machining processes in general [12,16,17], and in grinding processes in particular. Inasaki [18] was the first to simulate the workpiece sinusoidal speed variation effect, concluding that chatter can be avoided with short periods and large amplitudes of variation. Barrenetxea et al. [19] and Alvarez et al. [20] demonstrated the application of this technique for infeed and through-feed centerless grinding processes, respectively, by devising a time-domain dynamic approach for generating stability maps whose axes were the amplitude and frequency of the sinusoidal variation. The analysis concluded that the sinusoidal shape is the optimal speed variation signal. Those maps were validated experimentally, but the most important limitation of that method was the very long time required for obtaining a complete stability map due to the large amount of simulations needed.

Therefore, the goal of this paper is to obtain a method for computation of the instability lobes in the time domain with much lower computing effort compared to traditional ones. The method is based on the implicit subspace iteration presented by Zatarain et al. [21,22] for milling process analysis. The improvement comparing with the semi-discretization method is achieved by removing the need for calculating the complete transition matrix whose eigenvalues give the stability degree of parameter combinations.

This method is applied both in cylindrical and centerless infeed grinding processes and validated by comparing the results with previously published ones by Alvarez et al. [14] and Barrenetxea et al. [19,23], respectively, analyzing the accuracy and the calculation time reduction. New formulation has been included in the method due to the geometrical instabilities that appear during the centerless grinding processes.

#### **2. Materials and Methods**

The stability criterion for discretized time domain simulations is based on the Floquet theorem [24], as shown by Insperger et al. [11] for milling process and Alvarez et al. for grinding process [14]. This theorem states that when the coefficient matrices are periodic, the solutions of the homogeneous time periodic delay Equation: **.**

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

have the general form

$$\mathbf{x}(t) \,:= \, \varepsilon^{\lambda t} \mathbf{p}(t) \tag{2}$$

where **p**(*t*) is also periodic with principal time period *T* and being **A**(*t*) = **A**(*t*+*T*) and **B**(*t*) = **B**(*t*+*T*). Then, the solutions are exponentially increasing or decreasing periodical functions and the decay is described by the characteristic exponent λ.

After discretization, the solution at a certain moment and the solution one period later are related by the transition matrix **Z:**

$$\mathbf{q}\_{i+m} = \mathbf{Z}\mathbf{q}\_i\tag{3}$$

being **q**<sup>i</sup> the multidimensional state at time *t*<sup>i</sup> and *m* = Δ*t*/*T*. The eigenvectors of the transition matrix **Z** are the discretized counterpart of the principal solutions of (1) and the eigenvalues μ of **Z**, that are called characteristic multipliers, are related to the characteristic exponents [25]

$$\lambda \, \lambda \, \, = \, \frac{1}{T} (\ln \|\mu\| + i(\arg(\mu) + 2k\pi)) \tag{4}$$

where *k* corresponds to the harmonics involved.

The stability condition is that there are no eigenvalues with a magnitude larger than 1. Therefore, the interesting eigenvalues are those of highest magnitude, as they will characterize the stability of the process.

The formation of the transition matrix **Z** is usually a costly procedure. Zatarain et al. proposed a method based on implicit subspace iteration with the aim of avoiding the calculation of the transition matrix as well as the resolution of the corresponding eigenvalue problem [21]. Therefore, the computational cost can be reduced considerably.

Power iteration is the simplest way to obtain the highest magnitude eigenvalue and associated eigenvector of a matrix but is a relatively slow method. Subspace iteration consists of performing the power iteration of some vectors **S**<sup>0</sup> representing a guess of the dominant eigenvectors.

$$\mathbf{S}\_0 = \begin{bmatrix} \mathbf{s}\_{0,1}, \mathbf{s}\_{0,2}, \mathbf{s}\_{0,3}, \dots, \mathbf{s}\_{0,n} \end{bmatrix} \tag{5}$$

Initially there is no way to consider that the vectors in **S** are the eigenvectors corresponding to the largest magnitude eigenvalues, nor that they contain a large projection on these eigenvectors. Nevertheless, after each time period, the projection of all the vectors on the dominant eigenvectors will increase in proportion to their corresponding eigenvalues.

The use of a group of vectors (subspace **S**) instead of a single vector and the orthogonalization of these vectors with respect to the matrix **S** are the basis of the subspace iteration method and give rise to a reduction in the calculation time. Therefore, the procedure will consist of multiplying the subspace by the transition matrix Z once again and proceeding to the orthogonalization of the subspace after each or after several multiplications. The implicit subspace iteration method [21] follows that methodology but without calculating the transition matrix.

The multiplication of the subspace **S** by transition matrix **Z** is the first step and can be substituted by the time integration of the vectors in **S** by a fast procedure

$$\mathbf{V}\_{\mathrm{i}} = \begin{bmatrix} \mathbf{v}\_{\mathrm{i},1}, \mathbf{v}\_{\mathrm{i},2}, \mathbf{v}\_{\mathrm{i},3}, \dots, \mathbf{v}\_{\mathrm{i},\mathrm{u}} \end{bmatrix} = \mathbf{Z} \mathbf{S}\_{\mathrm{i}} \tag{6}$$

Now it is necessary to proceed to the orthogonalization of the vectors in **S**. The proposed method performs the orthogonalization, calculating an **H** matrix as

$$\mathbf{H}\_{i} = \left(\mathbf{S}\_{i}^{T}\mathbf{S}\_{i}\right)^{-1}\mathbf{S}\_{i}^{T}\mathbf{V}\_{i} \tag{7}$$

and then performing an eigenvalue decomposition on it as

$$\mathbf{H}\_{l} = \mathbf{G}\_{l}\mathbf{\hat{\lambda}}\_{l}\mathbf{G}\_{l}^{-1} \tag{8}$$

Finally, the calculation of the new subspace is obtained as

$$\mathbf{S}\_{i+1} = \mathbf{G}\_i \mathbf{V}\_i \tag{9}$$

After enough iterations, the subspace **S** will converge to the dominant eigenvectors, and the eigenvalues of λ will be the dominant eigenvalues.

With this method, it is possible not to form the transition matrix **Z**, but just to perform the time integration in one period (corresponding to a workpiece revolution) and obtain the new set of vectors **V**. The transition matrix is not required either to perform the orthogonalization of the subspace vectors.

The method of semi-discretization for calculating the eigenvalues requires calculating the transition matrix, which is highly time consuming, and then obtaining its eigenvalues, also time consuming. The implicit subspace iteration method avoids calculating the transition matrix, and, instead of performing multiplications of the subspace vectors by the transition matrix, calculates their time integration along a time period. As the transition matrix is not calculated and the eigenvalue systems to solve are very small compared to the size of the transition matrix, the calculation time is reduced very significantly. If we start with any initial solution for the state over a length of time identified by the maximum time delay (τmax), we can calculate the state at period *T* by multiplying the initial solution by the transition matrix or by performing the time integration over that period. The result of both methods must be the same, except for discrepancies given by numerical errors. Moreover, efficient time integration algorithms can result in much faster solutions, depending on the size of the state, than the time required to form the transition matrix.

The time integration along a period is a process that must be performed many times in this method. To reduce the calculation time, "step matrices" relating the state at each moment with the state one time step later are pre-calculated.

The formulation for the calculation of the step matrices in infeed grinding processes is shown next for the centerless configuration (Figure 1). In the case of cylindrical configuration, the procedure is the same just subtracting the geometrical parameters of the formulation.

**Figure 1.** Centerless infeed grinding process system.

Then, the centerless infeed grinding process is defined by an Equation of motion for a multiple degree of freedom system (corresponding to considered vibration modes):

$$\mathbf{M}\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + \mathbf{K}\mathbf{x}(t) \ = \mathbf{F}(t) \tag{10}$$

where **M**, **C** and **K** are the mass, viscous damping and stiffness matrices of the system, and **F**(*t*) is the grinding force in the normal direction, which follows the next expression:

$$\mathbf{F}(t) = \ \mathbf{\dot{q}} \mathbf{P} k\_{\text{inv}} [\delta r\_w(t - \tau) - \delta r\_w(t)] \tag{11}$$

Normalizing to unit mass, φ is a matrix with the modal displacements of considered modes in the normal direction of the contacts between the workpiece and the grinding wheel, the regulating wheel and the support blade. *kwn* is the cutting stiffness and δ*rw(t-*τ*)* and δ*rw(t)* are the radius defects at the cutting point in the previous and the current workpiece revolutions, respectively, being τ the time between revolutions. **P** is a dimensionless vector relating the forces at the contact points with the normal force at the cutting point [26]. Thus, the considered direction of motion of the process system is corresponding to the normal force **F**(*t*), defined as **x**(*t*).

In centerless grinding processes, the radius defect depends on a geometric displacement of the workpiece due to roundness errors passing through the contact points with the blade and regulating wheel, as well as on the static defection of the system and the vibrations generated in the process. The radius defect at the current revolution can be expressed as:

$$
\delta r\_w(t) = \mathbf{C} \Phi x(t) + \frac{\mathbf{F}(t)}{k\_I} + g\_b \delta r\_w(t - \tau\_b) - g\_r \delta r\_w(t - \tau\_r) \tag{12}
$$

where *kr* is the inverse of the residual flexibility. It can be calculated by subtracting the flexibility of considered vibration modes from the process equivalent stiffness *keq*, which is the sum of the system static flexibility and the contact flexibility between the workpiece and the grinding wheel, the regulating wheel and the support blade.

$$\frac{1}{k\_r} = \frac{1}{k\_{aq}} - \sum\_{i=1}^{N\_m} \frac{\mathbf{C} \phi\_i \phi\_i^{\prime} \mathbf{P}}{\omega\_i^2} \tag{13}$$

where ω*<sup>i</sup>* are the natural frequencies of considered modes. There are different experimental methods to obtain the equivalent stiffness [27]. δ*rw(t-*τ*b)* and δ*rw(t-*τ*r)* are the radius defects at the contact points with the blade and regulating wheel (represented in Figure 1 as δ*<sup>b</sup>* and δ*r*), and *gb* and *gr* are two geometrical parameters depending on the contact angles of the blade and regulating wheel with the workpiece (ϕ*<sup>b</sup>* and ϕ*r*), obtained from the geometric configuration of the process represented by the workpiece height *h* and the blade angle θ [28]. **C** is a dimensionless vector representing the displacement at the cutting point due to the displacements at the contact points [26].

Replacing Equations (11) and (12) in Equation (10) and renaming:

$$\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + [\mathbf{K} + \mathbf{V}k\_l]\mathbf{x}(t) \ = \mathbf{q}\Phi \mathbf{P}k\_l[\delta + \mathbf{g}\_b \delta\_{\mathbf{b}} - \mathbf{g}\_r \delta\_r] \tag{14}$$

being **<sup>V</sup>** <sup>=</sup> **<sup>C</sup>***T*φφ*T***P**, <sup>δ</sup>*<sup>i</sup>* <sup>=</sup> *rw*(*<sup>t</sup>* <sup>−</sup> <sup>τ</sup>), <sup>δ</sup>*<sup>b</sup>* <sup>=</sup> *rw*(*<sup>t</sup>* <sup>−</sup> <sup>τ</sup>*b*), <sup>δ</sup>*<sup>r</sup>* <sup>=</sup> *rw*(*<sup>t</sup>* <sup>−</sup> <sup>τ</sup>*r*) and *kt* <sup>=</sup> *kwn* 1+ *kwn kr* .

Following the same procedure of integrating and rearranging in state-space mode as [14] for cylindrical grinding process, the next expression is obtained:

$$\begin{array}{rcl} \mathbf{p}\_{i+1} &=& (\mathbf{A}\_i - \mathbf{B}\_i)\mathbf{p}\_i + \mathbf{G}\_{0,i} \Big( \delta\_{i-n} + g\_b \delta\_{b,i-n1} - g\_r \delta\_{r,i-n2} \Big) \\ &+ \mathbf{G}\_{1,i} \Big( \delta\_{i-n+1} + g\_b \delta\_{b,i-n1+1} - g\_r \delta\_{r,i-n2+1} \Big) \end{array} \tag{15}$$

$$\mathbf{B}\_{l} = \begin{bmatrix} \frac{\Delta t^{2}}{2} & \frac{\Delta t^{3}}{6} \\ \Delta t & \frac{\Delta t^{2}}{2} \end{bmatrix} \begin{bmatrix} \mathbf{V}k\_{l} & \mathbf{0} \\ \mathbf{0} & \mathbf{V}k\_{l} \end{bmatrix} \tag{16}$$

$$\mathbf{G}\_{0,i} = \begin{bmatrix} \frac{\Delta t^2}{2} \\ \frac{\Delta t}{2} \end{bmatrix} \blacksquare \Phi \mathbf{P} k\_t \; , \; \mathbf{G}\_{1,i} = \begin{bmatrix} \frac{\Delta t^2}{6} \\ \frac{\Delta t}{2} \end{bmatrix} \blacksquare \Phi \mathbf{P} k\_t \tag{17}$$

Matrix **A** is constant for all the calculations. It is obtained as [21]. **B**i, **G**0,i and **G**1,i vary and they are calculated for each discretized segment. This set of matrices represents the step matrices of the process.

Once all the step matrices are calculated, the procedure to obtain the state after one period begins with the definition of the initial state of the process by including the initial discretized radius defect and the initial dynamic modal solution, presented in the following mixed system state:

$$\mathbf{q}\_n = \begin{pmatrix} \delta\_1 & \delta\_2 & \dots & \delta\_{n-2} & \delta\_{n-1} & \mathbf{x}\_n & \dot{\mathbf{x}}\_n \end{pmatrix} \tag{18}$$

For the first segment:

$$\mathbf{p}\_1 = \mathbf{A}\mathbf{p}\_n\tag{19}$$

$$
\delta\_1 = \mathbf{T} \mathbf{p}\_1 + \mathcal{W} \delta\_{1-n} + G\_b \delta\_{1-n1} - G\_r \delta\_{1-n2} \tag{20}
$$

$$\mathbf{T} = \frac{\mathbf{C}^T \Phi}{\left(1 + \frac{k\_{\text{min}}}{k\_r}\right)} , \; W = \frac{k\_l}{k\_r} \tag{21}$$

Reaching the last segment:

$$\delta\_{i+1} = \mathbf{T} \mathbf{p}\_{i+1} + \mathcal{W} \delta\_{i+1-n} + G\_b \delta\_{i+1-n1} - G\_r \delta\_{i+1-n2} \tag{22}$$

$$\begin{array}{rcl} \mathbf{p}\_{i+1} &=& (\mathbf{A}\_i - \mathbf{B}\_i)\mathbf{p}\_i + \mathbf{G}\_{0,i} \big( \delta\_{i-n} + \mathcal{g}\_b \delta\_{b,i-n1} - \mathcal{g}\_r \delta\_{r,i-n2} \big) \\ &+ \mathbf{G}\_{1,i} \big( \delta\_{i-n+1} + \mathcal{g}\_b \delta\_{b,i-n1+1} - \mathcal{g}\_r \delta\_{r,i-n2+1} \big) \end{array} \tag{23}$$

The calculation of the product of the complete state vector using the step matrices, as proposed here, leads to a lower requirement of mathematical operations compared to the product by the transition matrix. Besides, the calculation of the transition matrix is the operation which requires the highest computational effort compared to the calculation of the step matrices when obtaining the stability maps for grinding processes.

Another important issue is the calculation time dependence on the size of the state, represented by the discretization segments. The requirement is that the integration time, equal to the length of the segment divided by the workpiece rotational speed, must be below 0.1 times the corresponding period related to the maximum natural frequency considered.

The criterion of stability is the same as the one used in the semi-discretization method. Then, if any eigenvalue of the relation between the initial and final states has a magnitude larger than 1, then the system is unstable, whereas when all magnitudes are lower than 1, the system is stable.

Besides dynamic instabilities, there are also geometric instabilities in centerless grinding processes due to their special geometric configuration [26], named geometric lobing. They are included in the model by the geometric parameters corresponding to the contact points between the workpiece and the regulating wheel and the blade support. For that reason, it is necessary to include different lobe numbers in the initial radius defect of the process. Then, the evolution of this lobe for each geometric configuration can be analyzed and included in the stability diagrams as a result of the subspace iteration method. Cylindrical infeed grinding processes do not give rise to that phenomenon.

As in the semi-discretization method in [14], with the developed method to analyze the stability in grinding processes there is also the possibility to include the continuous workpiece speed variation in the calculation.

The condition is that the speed variation period is an integer multiple of the workpiece speed period. Then, the calculation of the step matrices must be done over the speed variation period and the eigenvectors analysis is carried out between the first and the last system states.

Therefore, the same stability maps presented in [14,19] can be obtained, whose axes are the amplitude and frequency of the workpiece speed variation.

#### **3. Results**

#### *3.1. Application of the Method to Cylindrical Grinding Processes*

First application of the method is carried out for cylindrical infeed grinding process by comparing the results with the data proposed in Alvarez et al. [14]. In this analysis, a validation of the accuracy of the new method was done and a comparison of the computational effort of both semidiscretization and subspace iteration methods is presented.

The process conditions for stability analysis in [14] are shown in Table 1.


**Table 1.** Process conditions for stability analysis in cylindrical infeed grinding [14].

Figure 2a shows the stability map obtained with the subspace iteration method for the process conditions of Table 1, which represents the stability degree of each combination of workpiece speed and ground length. Values below the stability limit of 1 are stable combinations and values over the stability limit are unstable combinations. Almost the same stability map as the one obtained with the semi-discretization method is achieved. Figure 2b shows isolines with different degrees of stability. Dash lines correspond to subspace iteration calculations, while solid lines correspond to semidiscretization calculations. Good correlation between isoline values can be noticed.

**Figure 2.** (**a**) Subspace iteration stability map; (**b**) comparison of stability isolines between semi-discretization and subspace iteration methods.

Necessary computational effort for achieving the stability map has been reduced significantly. The time required for obtaining the map with the semi-discretization method is 11,441 s, while the time with the subspace iteration method is 256 s (45 times faster), considering that the map corresponds to a <sup>11</sup> <sup>×</sup> 11 matrix of stability values. For this comparison, an Intel ® Core TM i7-4610M CPU 3.00 GHz with RAM 16 GB has been used.

One important aspect is the dependency of the number of discretized points on the calculation time. As mentioned above, the integration time is calculated based on the considered maximum natural frequency. Figure 3 shows the comparison of the calculation time of both semi-discretization and subspace iteration methods related to that frequency. As it can be noticed, the higher the frequency the higher the difference between both calculations. For a natural frequency of 100 Hz, the subspace iteration method is six times faster, being 138 times faster for a frequency of 700 Hz. Therefore, the subspace iteration method is even more suitable for systems with higher natural frequencies.

**Figure 3.** Comparison of time for calculation of stability maps between semi-discretization and subspace iteration methods for cylindrical grinding.

The subspace iteration method is also valid to calculate the stability maps with continuous workpiece speed variations, whose axes are the frequency and amplitude of the sinusoidal variation. Figure 4 shows a comparison of the method with that of the semi-discretization method following the same procedure. Stability maps are calculated with process conditions of [14].

**Figure 4.** Comparison of instability reduction ratio maps with the application of continuous workpiece speed variation, obtained with both methods.

Isolines of instability reduction ratio are presented for both methods and good correlation can be noticed. Black boxes show the theoretical values of isolines while white ellipses show the experimental values.

The time required for obtaining the map with the semi-discretization method is 46,314 s, while the time with the subspace iteration method is 2773 s (16 times faster), considering that the map corresponds to a 16 × 26 matrix of instability reduction ratio values.

Calculation time for achieving this type of map is increased due to the necessity of calculating the step matrices over the speed variation period. Therefore, the higher the speed variation frequency, the lower the calculation time. Figure 5 shows the time required for frequencies of 0.5 and 3 Hz, where the matrices must be recalculated six and two times, respectively. In this case, for a natural frequency of 100 Hz, the subspace iteration method is 2.5 times faster, being 47 times faster for a frequency of 700 Hz.

**Figure 5.** Comparison of time for calculation of stability maps between semi-discretization and subspace iteration methods for cylindrical grinding with the application of continuous workpiece speed variation.

#### *3.2. Application of the Method to Centerless Grinding Processes*

Firstly, the application of the subspace iteration method to centerless grinding processes is validated by comparing it to the work by Barrenetxea et al. [23], where stability maps are calculated and validated in the frequency domain. Figure 6a shows the map obtained in the frequency domain and Figure 6b shows the one by the subspace iteration method. Good correlation between both can be noticed. The main difference is that the stability limit for frequency domain is equal to zero with negative unstable conditions, while subspace iteration is equal to 1 with unstable conditions over that value.

**Figure 6.** Dynamic stability map for centerless grinding by (**a**) frequency domain method [21]; (**b**) subspace iteration method.

In the case of centerless grinding processes, aside from the work done by Barrenetxea et al. [19], there is no bibliography regarding the obtainment of stability maps with continuous workpiece speed variation. Then, the subspace iteration method is compared to the time domain simulation procedure by Barrenetxea et al. [19] in terms of calculation accuracy and computational time consumption. Figure 7 shows the map obtained in [19] by the time simulation approach while Figure 8 shows the one obtained by the subspace iteration method. Good correlation can be noticed again, validating the accuracy of the proposed method.

**Figure 7.** Instability reduction ratio for centerless grinding by time simulation method [19].

**Figure 8.** Instability reduction ratio for centerless grinding by subspace iteration method.

Calculation time for achieving this type of instability reduction ratio in centerless grinding with the application of continuous workpiece speed variation is represented in the Figure 9. In this case, the time-domain simulations are independent of the variation frequency, since the simulation is done during all the workpiece revolutions of the whole grinding cycle. The time required for achieving the time-domain 11 × 11 map of Figure 7 is 63,525 s, while the time for achieving the subspace iteration map of Figure 8 is 2541 s, being 25 times faster.

**Figure 9.** Comparison of time for calculation of stability maps between time-domain and subspace iteration methods for centerless grinding with the application of continuous workpiece speed variation.

#### *3.3. Conclusions*

The results obtained from the comparison between subspace iteration and semi-discretization methods in cylindrical grinding process show that:


Regarding the results obtained comparing subspace iteration and time simulation methods in centerless grinding process:


#### **4. Discussion**

The subspace iteration method has been applied successfully to the generation of stability maps for cylindrical and centerless infeed grinding processes. The method has also been validated to calculate the instability reduction ratio maps with the application of continuous workpiece speed variation in cylindrical and centerless infeed grinding processes.

Accuracy of this method has been validated by comparing it with previously published methods such as semi-discretization, frequency and time-domain simulations, obtaining almost the same results and validating the subspace iteration method in grinding processes.

Reduction of calculation time related to the semi-discretization method is achieved in analyzed cases because the transition matrix is not explicitly calculated and, instead of calculating all the eigenvectors of the transition matrix, an iterative process calculating only the dominant eigenvalues is used. In the comparison with the semi-discretization method for cylindrical infeed grinding, subspace iteration is up to 138 times faster for stability maps and 47 times faster for instability reduction ratio maps. In the comparison with time-domain for centerless infeed grinding, subspace iteration is 27 times faster for instability reduction ratio maps.

The influence of the number of segments related to the maximum natural frequency of the grinding system has been analyzed, concluding that subspace iteration has a better performance for a higher number of segments than semi-discretization.

**Author Contributions:** Conceptualization, J.A. and M.Z.; Data curation, J.A. and B.I.; Formal analysis, M.Z. and J.I.M.; Investigation, J.A. and M.Z.; Methodology, D.B. and J.I.M.; Project administration, D.B.; Resources, J.A., D.B. and J.I.M.; Software, J.A.; Supervision, M.Z.; Validation, J.A., D.B. and B.I.; Visualization, J.A. and B.I.; Writing—original draft, J.A. and D.B.; Writing—review & editing, J.A. and M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **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/).
