## 2 Methods  
### Harmonic Balance Method  
Since traditional time integral method suffers from expensive cost when solving nonlinear equations, the HBM is applied in this paper, which mainly concentrates on the steady-state regimes of vibration response. In time domain, the equation of motion for a rotor rubbing with rigid stators can be expressed by:  
\[\mathbf{M}\ddot{\mathbf{x}}+(\mathbf{C}_{o}+\omega\mathbf{C}_{c})\dot{\mathbf{x }}+\Big{(}\mathbf{K}_{o}-\omega^{2}\mathbf{K}_{s}\Big{)}\mathbf{x}+\mathbf{f}_ {nl}(\mathbf{x},\lambda,t)+\mathbf{f}_{ub}(t)=0 \tag{1}\]  
where \(\mathbf{M}\), \(\mathbf{C}_{o}\) and \(\mathbf{K}_{o}\) are, respectively, the linear mass, damping and stiffness matrices of the system; \(\omega\) is the rotation speed of the system; \(\mathbf{x}\) is the vector of displacements, and, \(\ddot{\mathbf{x}}\), \(\ddot{\mathbf{x}}\) are the corresponding derivatives of displacements with respect to time \(t\); \(\omega\mathbf{C}_{c}\) is the Coriolis matrix; \(\omega^{2}\mathbf{K}_{s}\) is the matrix describing stiffness softening effect; \(\mathbf{f}_{ub}\) is the vector of unbalanced force; \(\mathbf{f}_{nl}\) is the nonlinear force existing in the rubbing interfaces between rotor and stator, which is dependent on the displacement vector \(\mathbf{x}\) and \(\lambda\), the parameter of contact interface used for sensitivity analysis. According to the HBM, the time-variant displacement vector \(\mathbf{x}\) in Equation (1) is transformed into frequency domain represented as a restricted number of Fourier series, i.e.,:  
\[\mathbf{x}=\mathbf{Q}_{0}+\sum_{j=1}^{n}\Big{[}\mathbf{Q}_{cj}\cos\big{(}m_{j }\tau\big{)}+\mathbf{Q}_{sj}\text{sin}(m_{j}\tau)\Big{]} \tag{2}\]  
where \(\mathbf{Q}_{0}\) is the vector of constant component in Fourier expansion; \(\mathbf{Q}_{cj}\) and \(\mathbf{Q}_{sj}\) are vectors of harmonic coefficients corresponding to cosine and sine terms respectively; \(m_{j}\) is the harmonic index; \(\tau\) is the non-dimensional time; \(n\) is the total number of harmonic terms retained, it should be selected carefully so that moderate number of harmonics are kept in harmonic expansion to reflect realistic vibration. The terms of nonlinear force and unbalanced force in Equation (1) are also expanded in frequency domain by sequentially multiplied by cosine terms \(\cos m_{j}\tau\) and sine terms \(\text{sin}m_{j}\tau\) for all harmonic indexes and integrated over the whole vibration period T, where T = 2\(\tau\)/\(\omega\). Since all time-variant terms have been expanded in the frequency domain, the multi-harmonic representation of Equation (1) can be obtained by equating harmonic coefficients, i.e.,:  
\[\mathbf{R}(\mathbf{Q})=\mathbf{Z}(\omega)\mathbf{Q}+\mathbf{F}_{n}+\mathbf{F} _{u}=0 \tag{3}\]  
where \(\mathbf{Q}=\Big{\{}\mathbf{Q}_{0}^{T},\mathbf{Q}_{c1}^{T},\mathbf{Q}_{s1}^{T},\cdots,\mathbf{Q}_{cn}^{T},\mathbf{Q}_{sn}^{T}\Big{\}}^{T}\) is the vector of harmonic components containing all coefficients for displacements described in Equation (2); \(\mathbf{F}_{n}\) and \(\mathbf{F}_{u}\) are vectors of multi-harmonic coefficients for rubbing forces and unbalanced forces; \(\mathbf{Z}(\omega)\) is the block-diagonal dynamic stiffness matrix and it has the following form:  
\[\mathbf{Z}(\omega)=\text{diag}(\mathbf{K},\mathbf{\Lambda}_{1},\cdots,\mathbf{ \Lambda}_{n}) \tag{4}\]  
\[\mathbf{\Lambda}_{i}=\left[ \begin{array}{cc}\mathbf{K}-(m_{i}\omega)^{2} \mathbf{M}&m_{i}\omega\mathbf{C}\\ -m_{i}\omega\mathbf{C}&\mathbf{K}-(m_{i}\omega)^{2}\mathbf{M}\end{array} \right] \tag{5}\]  
In Equation (5), \(\mathbf{K}\) and \(\mathbf{C}\) are simplified representation forms of generalized stiffness and damping matrices of system DOFs, with \(\mathbf{K}=\mathbf{K}_{o}-\omega^{2}\mathbf{K}_{s}\) and \(\mathbf{C}=\mathbf{C}_{o}+\omega\mathbf{C}_{s}\) respectively. The left-hand side of Equation (3) is defined as the residual terms which are used to describe if the sought vector converges to the actual solution.