## 4 Methods Used to Solve the Equation of Motion
### Harmonic Balance Method  
The HBM uses truncated Fourier series for solving the nonlinear differential equations by approximating the steady-state solution [36; 37] with a finite number of harmonic terms with amplitude coefficients and frequencies. The HBM involves approximating [36] the vectors of unknown kinematic parameters \(\mathbf{q},\ \dot{\mathbf{q}},\ \ddot{\mathbf{q}}\), the vector of nonlinear hydraulic forces \(\mathbf{f}_{\rm D}(\mathbf{q},\dot{\mathbf{q}})\), and the vector of external forces \(\mathbf{f}_{\rm A}\) in the equation of motion (1) with the terms of Fourier series:  
\[\mathbf{q}(t)=\sum_{n=1}^{N}\ddot{\mathbf{q}}_{n}\mathrm{e}^{\mathrm{i}(n \omega/v)t},\mathbf{f}_{\rm D}(t)=\sum_{n=1}^{N}\ddot{\mathbf{f}}_{\rm Dr} \mathrm{e}^{\mathrm{i}(n\omega/v)t},\mathbf{f}_{\rm A}(t)=\sum_{n=1}^{N} \ddot{\mathbf{f}}_{\rm An}\mathrm{e}^{\mathrm{i}(n\omega/v)t}, \tag{18}\]  
where \(\ddot{\mathbf{q}}_{n},\ \ddot{\mathbf{f}}_{\rm Dr}\), and \(\ddot{\mathbf{f}}_{An}\) are Fourier coefficients of the amplitude vector of the displacements, hydraulic, and external forces. The number of harmonic terms is denoted \(N\), the imaginary unit is \(\mathrm{i}\), and the optional integer parameter \(v\) allows to accommodate possible subharmonic multiples of the rotor rotational speed \(\omega\). By substituting relations (18) into the equation of motion (1), a set of nonlinear algebraical equations is obtained:  
\[\mathbf{\tilde{K}}\ \mathbf{\tilde{q}}_{n}=\mathbf{\tilde{f}}_{\rm Dn}\Big{(} \mathbf{\tilde{q}}_{n}\Big{)}+\mathbf{\tilde{f}}_{\rm A}, \tag{19}\]  
where the matrix of dynamical stiffness \(\mathbf{\tilde{K}}\) have the form  
\[\mathbf{\tilde{K}}=\bigg{[}-\Big{(}\frac{n\omega}{v}\Big{)}^{2}\mathbf{M}+ \mathrm{i}\frac{n\omega}{v}\mathbf{C}+\mathbf{K}\bigg{]}. \tag{20}\]  
The nonlinear set of algebraical equations (19) is solved for the unknown Fourier coefficients by the Newton-Raphson like method. In each iteration step, the Alternating Frequency-Time Scheme [37] is used to approximate the displacements in the time domain and evaluate the nonlinear force. The stable and unstable vibration response is determined by means of the continuity technique. The computational procedure uses a predictor-corrector method called the arc length method. In the corrector step, the nonlinear algebraical Equation (19) is parameterized with an additional equation [38]  
\[(\mathbf{q}-\mathbf{q}_{0})^{\rm T}(\mathbf{q}-\mathbf{q}_{0})+(\omega- \omega_{0})^{2}=s^{2}, \tag{21}\]  
where \(s\) is the arc-length parameter, \(\mathbf{q}\) and \(\omega\) are values in the current step, and \(\mathbf{q}_{0}\) and \(\omega_{0}\) denote values from the previous continuation step. The predictor step uses a secant predictor, where guess values for the following step are extrapolated from the two previous steps. The set of nonlinear equations (19) and the parametrization condition, Equation (21), often have elements with different orders of magnitude and appropriate scaling [37] should be applied. The vibration stability of the approximated solution in relations (18) is assessed by the application of Floquet theory on the perturbed equations of motion [39; 40]. The vector of nonlinear hydraulic forces \(\mathbf{f}_{\mathrm{D}}\big{(}\mathbf{q},\dot{\mathbf{q}}\big{)}\) is expanded into the Taylor series in the neighbourhood of the phase trajectory with omitting the higher terms:  
\[\mathbf{f}_{\mathrm{D}}\big{(}\mathbf{q}+\Delta\mathbf{q},\dot{\mathbf{q}}+ \Delta\dot{\mathbf{q}}\big{)}=\mathbf{f}_{\mathrm{D}}\big{(}\mathbf{q},\dot{ \mathbf{q}}\big{)}+\mathbf{J}_{\mathrm{K}}\Delta\mathbf{q}+\mathbf{J}_{ \mathrm{C}}\Delta\dot{\mathbf{q}} \tag{22}\]  
(23)  
where \(\mathbf{J}_{\mathrm{K}}\), \(\mathbf{J}_{\mathrm{C}}\) are Jacobi matrices of partial derivations and \(\Delta\mathbf{q}\), \(\Delta\dot{\mathbf{q}}\), are the vectors of deviations of displacements and velocities, respectively.