## 3 Computational Methods and Validation
### Stability Analysis Method  
Floquet theory provides a framework for analyzing periodic solutions in linear and nonlinear systems. The stability of these periodic solutions is determined through the  
Figure 4: Diagram of the arc-length continuation-based HB-AFT method. eigenvalues of the Floquet monodromy matrix, also known as Floquet multipliers. Based on this theory, Equation (1) can be rewritten as follows:  
\[\mathbf{M}_{s}\ddot{\mathbf{q}}_{0}+(\mathbf{C}_{s}-\Omega\mathbf{G}_{s}+\mathbf{C}_{base})\dot{\mathbf{ q}}_{0}+(\mathbf{K}_{s}+\mathbf{K}_{base})\mathbf{q}_{0}=\mathbf{F}_{mu}+\mathbf{F}_{base}+\mathbf{F}_{sfd} \tag{28}\]  
where \(\mathbf{q}_{0}\) represents the periodic solution with period \(T\) and \(N\) degrees of freedom, obtainable through the HB-AFT method. To investigate the stability of steady-state periodic solutions, a small perturbation \(\Delta\mathbf{q}(t)\) is superimposed on the periodic solution \(\mathbf{q}_{0}(t)\), yielding the perturbed periodic solution:  
\[\mathbf{q}(t)=\mathbf{q}_{0}(t)+\Delta\mathbf{q}(t) \tag{29}\]  
Substituting Equation (29) into Equation (28) and applying Taylor expansion around \(\mathbf{q}_{0}\) with linear approximation yields the perturbation differential equation:  
\[\mathbf{M}_{s}\Delta\ddot{\mathbf{q}}+(\mathbf{C}_{s}-\Omega\mathbf{G}+\mathbf{C}_{base}-\mathbf{C}_{sfd })\Delta\dot{\mathbf{q}}+(\mathbf{K}_{s}+\mathbf{K}_{base}-\mathbf{K}_{sfd})\Delta\mathbf{q}=0 \tag{30}\]  
where \(\mathbf{K}_{sfd}=\frac{\partial\mathbf{F}_{sfd}}{\partial\mathbf{q}}\Big{|}_{(\mathbf{q}_{0} \dot{\mathbf{q}}_{0})}\), \(\mathbf{C}_{sfd}=\frac{\partial\mathbf{F}_{sfd}}{\partial\dot{\mathbf{q}}}\Big{|}_{(\mathbf{q }_{0}\dot{\mathbf{q}}_{0})}\). The stability of the periodic solution \(\mathbf{q}_{0}\) depends on \(\Delta\mathbf{q}\). Based on Floquet theory [31], one period \(T\) is divided into \(M\) equal time intervals \(\Delta t\), assuming the monodromy matrix remains constant within each interval \(\Delta t\). The equation can be solved through numerical integration, requiring a sufficiently large value of \(M\). To balance computational efficiency and accuracy, \(M=400\) is selected. For the j-th time interval, the equation becomes:  
\[\mathbf{M}_{s}\Delta\ddot{\mathbf{q}}_{j}+\mathbf{C}_{j}\Delta\dot{\mathbf{q}}_{j}+\mathbf{K}_{j} \Delta\mathbf{q}_{j}=0 \tag{31}\]  
where \(\mathbf{C}_{j}=\mathbf{C}_{s}-\Omega\mathbf{G}_{s}+\mathbf{C}_{base}-\mathbf{C}_{sfd},\mathbf{K}_{j} =\mathbf{K}_{s}+\mathbf{K}_{base}-\mathbf{K}_{sfd}\). Primary numerical methods for approximating the monodromy matrix include the Runge-Kutta method [4], Newmark method [32], and Hsu exponential expansion method [33]. Notably, for computing the monodromy matrix of multi-degree-of-freedom finite element rotor systems, the Newmark method demonstrates superior computational efficiency [6]. Therefore, the Newmark numerical integration method is employed to solve for Floquet multipliers. Under the Newmark scheme, the expressions for generalized velocity and acceleration of the perturbation differential equation are as follows:  
\[ \begin{cases}\Delta\dot{\mathbf{q}}_{j+1}=-\Delta\dot{\mathbf{q}}_{j}+\frac{2}{\Delta t }(\Delta\mathbf{q}_{j+1}-\Delta\mathbf{q}_{j})-\Delta t\Delta\ddot{\mathbf{q}}_{j}\\ \Delta\ddot{\mathbf{q}}_{j+1}=-\Delta\ddot{\mathbf{q}}_{j}+\frac{4}{\Delta t^{2}}( \Delta\mathbf{q}_{j+1}-\Delta\mathbf{q}_{j}-\Delta t\Delta\dot{\mathbf{q}}_{j})\end{cases}  \tag{32}\]  
Substituting Equation (32) into Equation (31) yields:  
\[\mathbf{D}_{j}\Delta\mathbf{q}_{j+1}=\mathbf{M}_{s}\Bigg{(}\Delta\ddot{\mathbf{q}}_{j}+\frac {4\Delta\dot{\mathbf{q}}_{j}}{\Delta t}+\frac{4\Delta\mathbf{q}_{j}}{\Delta t^{2}} \Bigg{)}+\mathbf{C}_{j+1}\Bigg{(}\Delta t\Delta\ddot{\mathbf{q}}_{j}+\Delta\dot{\mathbf{q }}_{j}+\frac{2\Delta\mathbf{q}_{j}}{\Delta t}\Bigg{)} \tag{33}\]  
where \(\mathbf{D}_{j}=\frac{4\mathbf{M}_{s}}{\Delta t^{2}}+\frac{2\mathbf{C}_{j+1}}{\Delta t}+\bm {K}_{j+1}\).