2.3.2. Algebraic Identifier with Variable Operation Velocity

In this section, the rotor-bearing system velocity is considered as a linear ramp of excitation. The mathematical model of the system is defined by Equation (3). In order to develop the parameter identifier, this equation is rewritten as follows

$$\left[\mathbf{M}\right]\left\{\ddot{\delta}\right\} + \left[\mathbf{C}\_1 + \dot{\phi}\mathbf{C}\_2\right]\left\{\dot{\delta}\right\} + \left[K\_1 + K\_2 + \ddot{\phi}K\_3\right]\{\delta\} = \dot{\phi}^2 F\_1(\phi) + \ddot{\phi}F\_2(\phi) \tag{16}$$

By multiplying Equation (16) by *t* <sup>2</sup> and integrating the result twice with respect to time, the following is obtained

$$\int\_{}^{\langle 2 \rangle} \left[ [M] \left\{ \bar{\delta} \right\} + \left[ \mathbb{C}\_1 + \dot{\phi} \mathbb{C}\_2 \right] \left\{ \dot{\delta} \right\} + \left[ K\_1 + K\_2 + \bar{\phi} K\_3 \right] \{ \delta \} \right] t^2 = \int\_{}^{\langle 2 \rangle} \left\{ \dot{\phi}^2 \mathbb{F}\_1(\phi) + \bar{\phi} \mathbb{F}\_2(\phi) \right\} t^2 \tag{17}$$

where 0 (2) *ϕ*(*t*) are iterated time-integrals of the form 0 *<sup>t</sup>* 0 0 *<sup>σ</sup>*<sup>1</sup> <sup>0</sup> ··· <sup>0</sup> *<sup>σ</sup>n*−<sup>1</sup> <sup>0</sup> *ϕ*(*σn*)*dσ<sup>n</sup>* ··· *dσ*<sup>1</sup> with 0 *ϕ*(*t*) = 0 *<sup>t</sup>* <sup>0</sup> *ϕ*(*σ*)*dσ*, and *n* a positive integer.

Similarly for the case of constant velocity, matrices [*K*1] and [*C*1] contain the stiffness and damping parameters provided by the supports. Therefore, after the integration of the left part of Equation (17) and rearranging terms, we have

$$\begin{aligned} &\int \prescript{(2)}{}{\{\mathcal{K}\}}\_{1} t^{2} \{\delta\} + [\mathcal{C}\_{1}] \left[ \int t^{2} \{\delta\} - 2 \int ^{(2)}{}{t} \{\delta\} \right] = \int \left[ 4Mt - \dot{\phi} \mathcal{C}\_{2} t^{2} \right] \{\delta\} \\ &+ \int \prescript{(2)}{}{\{\mathcal{C}\_{2} \left(\ddot{\phi} t^{2} + 2\dot{\phi} t\right) - 2M - \left(K\_{2} + \ddot{\phi} K\_{3}\right) t^{2}\}} \{\delta\} - [M] t^{2} \{\delta\} \\ &+ \int \prescript{(2)}{}{\{\dot{\phi}^{2} F\_{1}(\phi) + \ddot{\phi} F\_{2}(\phi)\}} t^{2} \end{aligned} \tag{18}$$

It is worth mentioning that Equation (18) can be separated into individual equation systems for each node where the bearings are located. These equations can be written as follows

$$
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zx} \end{bmatrix} \int^{\langle 2 \rangle} t^2 \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} + \begin{Bmatrix} c\_{xx} & c\_{xz} \\ c\_{zx} & c\_{xz} \end{Bmatrix} \left( \int t^2 \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} - 2 \int^{\langle 2 \rangle} t \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} \right) = \begin{Bmatrix} b\_{ni} \\ b\_{ni} \end{Bmatrix} \tag{19}
$$

To solve Equation (19), an equal number of equations and unknows is required. Therefore, Equation (19) is successively integrated three times to obtain the missing equations which are expressed as

$$
\begin{aligned}
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zz} \end{bmatrix} \end{aligned}
\begin{aligned}
\begin{aligned}
\begin{Bmatrix} (3) \\ w\_{i} \end{Bmatrix} \end{aligned} \\
\begin{aligned}
\begin{Bmatrix} \boldsymbol{w}\_{i} \\ \boldsymbol{w}\_{i} \end{Bmatrix} \end{aligned} + \begin{bmatrix} \mathbf{c}\_{xx} & \mathbf{c}\_{xz} \\ \mathbf{c}\_{xx} & \mathbf{c}\_{zz} \end{bmatrix} \left( \int \begin{array}{c} \mathbf{t}^{(2)} \\ \mathbf{w}\_{i} \end{array} \begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \right) - 2 \int \begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned} \\
\begin{aligned}
\begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned} \tag{20}
\end{aligned}
\end{aligned}
\end{aligned}
\begin{aligned}
\begin{aligned}
\begin{Bmatrix} \mathbf{w}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned}$$

$$
\begin{bmatrix} k\_{\text{xx}} & k\_{\text{xz}} \\ k\_{\text{zx}} & k\_{\text{zz}} \end{bmatrix} \int^{(4)} \boldsymbol{t}^2 \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} + \begin{bmatrix} \boldsymbol{c}\_{\text{xx}} & \boldsymbol{c}\_{\text{xz}} \\ \boldsymbol{c}\_{\text{zx}} & \boldsymbol{c}\_{\text{zz}} \end{bmatrix} \left( \int^{(3)} \boldsymbol{t}^2 \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} - 2 \int^{(4)} \boldsymbol{t} \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} \right) = \int^{(2)} \left\{ \begin{array}{c} \boldsymbol{b}\_{\text{ui}} \\ \boldsymbol{b}\_{\text{ui}} \end{array} \right\} \tag{21}
$$

$$
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zz} \end{bmatrix} \int \stackrel{(5)}{\cdot} t^2 \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} + \begin{bmatrix} c\_{xx} & c\_{xz} \\ c\_{zx} & c\_{zz} \end{bmatrix} \left( \int \stackrel{(4)}{\cdot} t^2 \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} - 2 \int ^{(5)} t \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} \right) = \int ^{(3)} \left\{ \begin{array}{cc} b\_{\text{ini}} \\ b\_{\text{ini}} \end{array} \right\} \tag{22}
$$

From Equations (19)–(22), a linear system equation is obtained for each node where the bearings are located. These equations can be expressed as

$$[A\_{\mathfrak{s}}(t)]\{\Theta\_{\mathfrak{s}}\} = \{b\_{\mathfrak{s}}(t)\}\tag{23}$$

where {Θ*s*} <sup>=</sup> {*kxx kxz kzx kzz cxx cxz czx czz* }*<sup>T</sup>* denotes the transposed vector of parameters to be identified and [*As*(*t*)], {*bs*(*t*)} are 8 × 8 and 8 × 1, respectively.

Again, the condition *det*[*As*(*t*)] = 0 must be satisfied to identify the vector {Θ*s*}.

From the solution of Equation (23), a mathematical model for an on-line identifier of stiffness and damping bearing parameters can be obtained as

$$\{\Theta\_{\mathfrak{s}}\} = \left[A\_{\mathfrak{s}}^{-1}\right] \{b\_{\mathfrak{s}}\} \,\,\forall t \in (t\_{0}, t\_{0} + \mathfrak{e}] \tag{24}$$

As can be observed, algebraic identification of stiffness and damping bearing parameters is independent of system initial conditions and only depends on the displacement vector and the type of ramp excitation. It is important to mention that as with the case of constant velocity, to identify the supports parameters, lateral vibration measurements at the node and the nodal slopes are required. Moreover, similar information from the adjacent node is also needed.

## **3. Results**

In Figure 3, a scheme of the rotor-bearing system considered in this work and its discretization is presented.

**Figure 3.** Rotor-bearing system scheme [40].

To obtain the mathematical model for the rotor-bearing system using the FE method, it was discretized into 11 beam-like elements, as is shown in Figure 3. The system includes two inertial disks located at nodes 1 and 12, while supports (bearings) are placed at nodes 4 and 8. The correct nodal location ensures that the simulation replicates the model's real conditions. In addition, two unbalanced masses were considered in two different angular positions located on inertial disks *D*<sup>1</sup> and *D*2.

In Table 1, the mechanical and geometrical properties of the shaft are shown, while the inertial properties of discs and unbalanced masses are presented in Table 2.


**Table 1.** Mechanical and geometrical properties of the rotor-bearing shaft.


**Table 2.** Inertial properties of the disks and unbalance masses.

In Table 3, the stiffness and damping bearing parameters [40] used for numerical simulation are presented.


**Table 3.** Stiffness and damping bearing parameters [40].

On-line algebraic identification of stiffness and damping bearing parameters was determined based on the vibratory response of the rotor-bearing system in the time domain, which was obtained from Equations (3) and (7) by using the Newmark method for numerical integration.

#### *3.1. Algebraic Parameter Identification with Constant System Velocity*

The displacement vector used in the algebraic identification procedure was obtained from Equation (7) by using the Newmark method for numerical integration and taking into account a constant rotational velocity of the rotor-bearing system.

In Figure 4, vibration signals at node 4 (corresponding to bearing 1 location) of the rotor-bearing system of Figure 3 are presented. This response is obtained for an operation rotational velocity Ω = 600 rpm. These signals, the nodal slopes and the corresponding information of the nodes adjacent to the bearing locations are the required data to identify stiffness and damping parameters of the bearings.

**Figure 4.** Vibration signal at node 4 (bearing 1) at 600 rpm: (**a**) X direction; (**b**) Z direction.

Figures 5 and 6 present the obtained results from the numerical simulation for the algebraic identification of stiffness and damping parameters for bearing 1, while Figures 7 and 8

show the results corresponding to bearing 2. It is important to mention that the sample time used in the simulation was 0.1 milliseconds. However, by carrying out numerical simulations with different sample times, it was observed that the shorter the sampling period, the faster the identifier converges.

**Figure 5.** Identified stiffness parameters for bearing 1 at 600 rpm. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 6.** Identified damping parameters for bearing 1 at 600 rpm. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 7.** Identified stiffness parameters for bearing 2 at 600 rpm. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 8.** Identified damping parameters for bearing 2 at 600 rpm. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

As can be observed in Figures 5–8, the identification of both stiffness and damping parameters of the bearings is carried out in less than 0.1 s, and once the parameter reaches the identified value, this remains for the rest of the time period. For a better analysis of the identifier behavior, only results for 0.1 s are presented in Figures 5–8, because it is important to observe the time that the identifier requires to converge to the estimated value.

### *3.2. Algebraic Parameter Identification with Variable System Velocity*

The displacement vector used as input data for the algebraic identification is obtained from Equation (3) by using the Newmark method for numerical integration and taking into account a linear ramp excitation with angular acceleration .. *φ* = 10 *rad*/*s*2. The rotor-bearing system response at node 4 is shown in Figure 9 where the vibratory behavior of the system in the location of bearing 1 can be appreciated.

**Figure 9.** System vibratory response at node 4 with a linear ramp of excitation of 10 *rad*/*s*2.

In Figures 10–13, the behavior of the algebraic identifier for bearing stiffness and damping parameters of both bearings (placed at nodes 4 and 8) is shown as a function of time.

**Figure 10.** Identified stiffness parameters for bearing 1 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 11.** Identified damping parameters for bearing 1 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 12.** Identified stiffness parameters for bearing 2 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 13.** Identified damping parameters for bearing 2 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

As in the case of constant rotor system velocity, stiffness and damping bearing parameters are identified in les tan 0.1 s as can be observed in Figures 10–13. Furthermore, the parameter values remain constant until the rotor-bearing system reaches its nominal operation velocity. For a better analysis of the identifier behavior, the results for 0.1 s are presented in Figures 10–13 because it is important to observe the required time for the identifier convergence. The sample time used to solve Equation (3) using the Newmark method was 0.1 milliseconds. The numerical solution of Equation (3) was used as input data for the proposed algebraic identifier. Moreover, achieving this sample time with diverse data acquisition systems for experimental implementation was verified.

#### **4. Discussion**

Different numerical simulations were carried out in order to determine the robustness of the proposed identifiers under different conditions for the rotor-bearing system velocity. For the constant velocity case, different magnitudes for the rotor system velocity were considered, while for the variable velocity case, different ramps of excitation were explored.

Figure 14 shows results for the algebraic identification of stiffness and damping parameters for bearing 1 at a constant operation velocity of the rotor-bearing system of 50,000 rpm. A rapid identifier convergence to the estimated values can be observed, meaning that an increase in operation velocity does not affect the identifier performance. It is important to mention that, while the results for the identification of damping parameters of bearing 1 and the stiffness and damping parameters of bearing 2 are not presented, these parameters are correctly identified in less than 0.1 s.

The identifier performance for different ramps' excitation was analyzed. The acceleration values considered for numerical simulation were: .. *φ* = 10 *rad*/*s*2, .. *<sup>φ</sup>* <sup>=</sup> <sup>100</sup> *rad*/*s*2, .. *φ* = 1000 *rad*/*s*2, .. *<sup>φ</sup>* <sup>=</sup> <sup>3000</sup> *rad*/*s*<sup>2</sup> and .. *<sup>φ</sup>* <sup>=</sup> <sup>6000</sup> *rad*/*s*2. The result for .. *φ* = 10 *rad*/*s*<sup>2</sup> were reported in the previous section. Due to the similar behavior of the identifier with different acceleration values, only results for .. *φ* = 6000 *rad*/*s*<sup>2</sup> are shown here. The rotor-bearing system response for a ramp of excitation with the mentioned value of acceleration at node 4 (bearing 1 location) is presented in Figure 15. It can be seen that there is a considerable change in the time scale in comparison with Figure 9 because the acceleration is increased 600 times.

**Figure 14.** Identified stiffness parameters for bearing 1. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz* at 50,000 rpm.

**Figure 15.** System vibratory response at node 4 with a linear ramp excitation with acceleration of 6000 *rad*/*s*2.

The algebraic identification performance under the conditions described above is shown in Figures 16–19 where the estimation for the stiffness and damping bearings parameter is visualized. For this simulation, the system response in Figure 15 is used as input data.

**Figure 16.** Identified stiffness parameters for bearing 1 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (d) *kzz*.

**Figure 17.** Identified damping parameters for bearing 1 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 18.** Identified stiffness parameters for bearing 2 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (d) *kzz*.

**Figure 19.** Identified damping parameters for bearing 2 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

From the results presented in Figures 16–19, it can be observed that despite the linear ramp excitation being 600 times faster than the corresponding one in Figure 9, the proposed identifier rapidly converges to the estimated values and remains in these values for the rest of the time period. Note that the algebraic identifier is not affected by the system's acceleration and only depends on the displacement vector at each instance. The robustness of the algebraic identification method against acceleration ramp variations had already been proved by Mendoza-Larios et al. [36] but only for the identification of unbalance parameters in rotor-bearing parameters.

Furthermore, the obtained results for both cases, constant and variable rotor-bearing velocity, show a transient state of the identifiers before the convergence to the estimated values of stiffness and damping bearing parameters. This behavior is due to the sample time used in the numerical simulations of the identifiers for solving the iterated integrals of Equations (15) and (24), which utilize the trapezoidal rule. According to Kharab and Guenther [41], this method presents major calculation errors in comparison with other integration methods. However, it was found that the smaller the sample time the shorter the error in the trapezoidal rule calculation.

#### **5. Conclusions**

In this article the identification problem for stiffness and damping parameters of the supports in rotor-bearing systems was addressed. The system model was obtained by the finite element method and using a finite element type beam, which consider the effects of rotational inertia, gyroscopic moments and shear deformations. The algebraic identification technique was applied to the finite element model to obtain two identifiers for the stiffness and damping parameters attributable to the bearings. The first identifier considers a rotor-bearing system operating at a constant velocity, and the second with a linear ramp of excitation as a system velocity input. The numerical results present the identifier behavior showing a fast convergence and robustness in both operation conditions with different values of constant rotational velocity and ramp of acceleration. The numerical results indicate a fast convergence in the stiffness and damping parameters identification in less than 0.06 s for both considered operation conditions. It is important to mention that the convergence time of the identifier depends mainly on the sample time used in numerical simulations. An important characteristic of the proposed algebraic identifiers is that the unbalance parameters (magnitude and phase) are not needed for their development and implementation because only the vibratory response of the system at the bearings' location and the adjacent nodes is required. As a first approach we have proved the proposed identifiers in rotor-bearing system models with constant rotordynamic coefficients. However, as a future work, the proposed identifiers can be used to numerically and experimentally determine rotordynamic coefficients, which are a function of the system rotational velocity as in the case of pressurized bearings, by adapting the identifier method to estimate non-constant functions. This is possible because in the mathematical model used for the identifiers' development, only the effects (stiffness and damping) of the supports are considered without taking into account the nature of the bearings.

**Author Contributions:** Conceptualization, J.G.M.-L., M.A.-M. and S.J.L.-D.; methodology, J.G.M.-L., E.B. and J.C.-O.; software, J.G.M.-L. and S.J.L.-D.; validation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; formal analysis, J.G.M.-L., M.A.-M. and R.T.-H.; investigation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; resources, J.G.M.-L., E.B., M.A.-M., R.T.-H. and J.C.-O.; data curation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; writing original draft preparation, J.G.M.-L., E.B. and M.A.-M.; writing—review and editing, J.G.M.-L., E.B. and M.A.-M.; supervision, E.B. and M.A.-M.; project administration, J.G.M.-L.; funding acquisition, J.G.M.-L. and M.A.-M. All authors have read and agreed to the published version of the manuscript.

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

**Data Availability Statement:** The datasets generated and supporting the findings in the article are obtainable from the corresponding author upon reasonable request.

**Acknowledgments:** The APC was funded by PRODEP-SEP, Mexico.

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