# Implicit Subspace Iteration to Improve the Stability Analysis in Grinding Processes
## 2 Materials and Methods  
The stability criterion for discretized time domain simulations is based on the Floquet theorem [24], as shown by [PERSON] et al. [11] for milling process and [PERSON] 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)\ =\ e^{\lambda t}\mathbf{p}(t) \tag{2}\]  
where \(\mathbf{p}(t)\) is also periodic with principal time period \(T\) and being \(\mathbf{A}(t)=\mathbf{A}(t+T)\) and \(\mathbf{B}(t)=\mathbf{B}(t+T)\). Then, the solutions are exponentially increasing or decreasing periodical functions and the decay is described by the characteristic exponent \(\lambda\). After enough iterations, the subspace **S** will converge to the dominant eigenvectors, and the eigenvalues of \(\lambda\) 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 (\(\tau_{\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. 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)\ =\ \Phi\mathbf{P}k_{wn}[\delta r_{w}(t-\tau)-\delta r_{w}(t)] \tag{11}\]  
Figure 1: Centerless infeed grinding process system. Normalizing to unit mass, \(\Phi\) 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.