## 3 Computational Methods and Validation  
### Harmonic Balance Method with Alternating Frequency/Time-Domain Technique  
Since the system response contains an excitation frequency related to the rotational speed, the periodic solution of the rotor system and nonlinear oil film forces can be represented by H-order truncated Fourier series:  
\[ \begin{cases}\mathbf{q}=\mathbf{a}_{0}+\sum\limits_{j=1}^{H}\mathbf{a}_{j}\text{cos}(j \Omega t)+\mathbf{b}_{j}\text{sin}(j\Omega t)\\ \mathbf{F}_{sfd}=\mathbf{c}_{0}+\sum\limits_{j=1}^{H}\mathbf{c}_{j}\text{cos}(j\Omega t)+ \mathbf{d}_{j}\text{sin}(j\Omega t)\end{cases}  \tag{13}\]  
where \(\mathbf{a}_{0}\) and \(\mathbf{c}_{0}\) are constant terms, and \(\mathbf{a}_{j}\), \(\mathbf{b}_{j}\), \(\mathbf{c}_{j}\), \(\mathbf{d}_{j}\) are the coefficients to be determined for each truncated Fourier term. Substituting Equation (13) into Equation (1) and eliminating the cosine and sine harmonic terms yields the system residual equation:  
\[\mathbf{R}(\mathbf{X})=\mathbf{P}\mathbf{X}+\mathbf{Q}_{nl}-\mathbf{Q}_{ex} \tag{14}\]  
where the coefficient matrix \(\mathbf{P}\) is expressed as follows:  
\[ \begin{split}\mathbf{P}=\text{diag}\Big{(}\mathbf{P}_{0}\quad\mathbf{P}_{1} \quad\cdots\quad\mathbf{P}_{j-1}\quad\mathbf{P}_{j}\quad\mathbf{P}_{j+1}\quad\cdots\quad\bm {P}_{H}\Big{)}\\ \mathbf{P}_{j}=\Big{[}\mathbf{K}_{s}+\mathbf{K}_{brg}+\mathbf{K}_{base}-(j\Omega)^{2}\mathbf{M}_{ s}\Big{]}-\text{i}\Big{[}j\Omega\Big{(}\mathbf{C}_{brg}+\mathbf{C}_{base}-\Omega\mathbf{G}_{s} \Big{)}\Big{]}\end{split}  \tag{15}\]where \(\text{diag}()\) denotes the diagonal arrangement of the enclosed matrices to form a block diagonal matrix, and i is the imaginary unit. The excitation force term \(\mathbf{Q}_{ex}\) is given as follows:  
\[\mathbf{Q}_{ex}=[\mathbf{Q}_{0};\mathbf{Q}_{1};\mathbf{0};\cdots;\mathbf{0}] \tag{16}\]  
where \(\mathbf{Q}_{0}=F_{base}\), \(\mathbf{Q}_{1}=F_{mn,c}+\text{i}F_{mn,s}\), and \(\mathbf{0}\) represents the zero vector. The analytical expression of nonlinear squeeze film forces is difficult to expand through trigonometric function relations. Therefore, the alternating frequency time-domain (AFT) method is required to obtain discrete time-domain information of nonlinear forces. The discrete time-domain information of system response is obtained through inverse discrete Fourier transform (IDFT):  
\[\mathbf{x}(n)=\text{Re}\left(\sum\limits_{j=0}^{H}\mathbf{X}_{j}\mathbf{e}^{\text{i}(2\pi jn )/N_{f}}\right) \tag{17}\]  
where \(\mathbf{X}_{j}=\mathbf{a}_{j}-\text{i}\mathbf{b}_{j}\); \(N_{f}\) is the number of selected time discretization points, which must be greater than 2 times of the total harmonic number \(H\); and \(\mathbf{x}(n)\) represents the time-domain information of generalized displacement Fourier coefficients at the \(n\)-th time point, \(n=0,1,2,...,N_{f}-1\). Similarly, the expression for \(\dot{\mathbf{x}}(n)\) can be obtained. Then, the time-domain information of nonlinear force Fourier coefficients \(\mathbf{\widetilde{f}}_{nl}\big{(}\dot{\mathbf{x}}(n),\mathbf{x}(n)\big{)}\) is determined through the explicit expression presented in Equation (12). The frequency-domain information of nonlinear force Fourier coefficients can be obtained using discrete Fourier transform (DFT):  
\[\mathbf{Q}_{nl}^{j}=\mathbf{c}_{j}-\text{i}\mathbf{d}_{j}=\left\{ \begin{array}{l}\frac{ 1}{N_{f}}\sum\limits_{n=0}^{N_{f}-1}\text{Re}(\mathbf{\widetilde{f}}_{nl}(\dot\mathbf{{x}}(n),x(n))\mathbf{e}^{\text{i}(-2\pi jn)/N_{f}}),j=0\\ \frac{2}{N_{f}}\sum\limits_{n=0}^{N_{f}-1}\mathbf{\widetilde{f}}_{nl}(\dot{\mathbf{x }}(n),x(n))\mathbf{e}^{\text{i}(-2\pi jn)/N_{f}},j\
eq 0\end{array} \right. \tag{18}\]  
Then, the nonlinear force term \(\mathbf{Q}_{nl}\) in Equation (14) can be obtained by \(\mathbf{Q}_{nl}=\left[\mathbf{Q}_{nl}^{0};\ \ \mathbf{Q}_{nl}^{1};\ \ \cdots\ \ \mathbf{Q}_{nl}^{H-1};\ \ \mathbf{Q}_{nl}^{H}\right]\). Nonlinear systems (such as the Duffing oscillator [29]) exhibit solution trajectories characterized by bifurcation phenomena and turning points.