## 4 Solution of Equations Using the Newmark \(\beta\) Method  
Due to the significant impact of nonlinearity in the study of SFD on rotor response, the Newmark method was employed to solve the control equations established previously, which include the SFD oil film force and the nonlinear restoring force of the supports. After discretizing the Euler-Bernoulli beam vibration equation using the finite element method and similar approaches, the vibration equation takes the following form:  
\[\mathbf{M\ddot{u}}(t)+\mathbf{C\ddot{u}}(t)+\mathbf{f}(t)=\mathbf{F}(t) \tag{61}\] \begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{**Order Number**} & **Critical Speed from Custom Code** & **Critical Speed from** & \multirow{2}{*}{**Difference (\%)**} \\  & **Solution (rpm)** & **ANSYS Solution (rpm)** & \\ \hline First Order & 2071.14 (34.52 Hz) & 2077.8 (34.63 Hz) & 0.32 \\ \hline Second Order & 3929.41 (65.50 Hz) & 3973.8 (66.23 Hz) & 1.11 \\ \hline Third Order & 4587.68 (76.46 Hz) & 4621.2 (77.02 Hz) & 0.72 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Comparison of critical speeds for a beam. Figure 3: Comparisons between custom code computations and ANSYS simulations. Figure 2: Rotor-support system model. It is commonly assumed that the mass matrix \(\mathbf{M}\) and the damping matrix \(\mathbf{C}\) are constant, while \(\mathbf{f}(t)\) represents the nonlinear restoring force. At time \(t+\Delta t\), the dynamic equilibrium equation must also be satisfied:  
\[\mathbf{M}\ddot{\mathbf{u}}(t+\Delta t)+\mathbf{C}\ddot{\mathbf{u}}(t+\Delta t )+\mathbf{f}(t+\Delta t)=\mathbf{F}(t+\Delta t) \tag{62}\]  
Expanding the restoring force \(\mathbf{f}(t+\Delta t)\) in a Taylor series around time \(t\), and neglecting higher-order terms, yields  
\[\mathbf{f}(t+\Delta t)\doteq\mathbf{f}(t)+\mathbf{K}(t)\Delta u \tag{63}\]  
where \(\mathbf{K}(t)\) is the tangent stiffness matrix at time \(t\). It can be observed that Equation (63) utilizes the tangent stiffness matrix at time \(t\), which is an approximation; hence, it introduces an error. To refine the solution, the Newton-Raphson method is employed for iteration, progressively converging towards the true solution. Substituting Equation (63) into (62) results in  
\[\mathbf{M}\ddot{\mathbf{u}}(t+\Delta t)+\mathbf{C}\ddot{\mathbf{u}}(t+\Delta t )+\mathbf{K}(t)\Delta\mathbf{u}=\mathbf{F}(t+\Delta t)-\mathbf{f}(t) \tag{64}\]  
In Equation (64), there are 3n unknowns \(\ddot{\mathbf{u}}(t+\Delta t)\), \(\dot{\mathbf{u}}(t+\Delta t)\) and \(\Delta\mathbf{u}\), yet there are only three equations, necessitating the addition of 2n constraint conditions. The Newmark Beta method is proposed to construct these constraint conditions, which will be elaborated later. Once the constraint conditions are incorporated, \(\ddot{\mathbf{u}}(t+\Delta t)\), \(\dot{\mathbf{u}}(t+\Delta t)\) and \(\Delta\mathbf{u}\) can be determined, thus providing an approximate value for the displacement at time \(t+\Delta t\). \[\mathbf{u}(t+\Delta t)=\mathbf{u}(t)+\Delta\mathbf{u} \tag{65}\]  
Assuming that a precise solution has been obtained at time \(t\), the iterative process at time \(t+\Delta t\) in accordance with the Newton-Raphson method adopts the state at time \(t\) as the initial condition, that is,  
\[\mathbf{u}^{(0)}(t+\Delta t)=\mathbf{u}(t)\] \[\mathbf{f}^{(0)}(t+\Delta t)=\mathbf{f}(t) \tag{66}\] \[\mathbf{K}^{(0)}(t)=\mathbf{K}(t)\]  
For the \(k=1\), \(2\),...iteration steps, it follows that  
Figure 3:0  
Figure 3:1  
During the first iteration step \(k=1\), Equations (67) and (68) simplify to Equations (64) and (65), respectively. In subsequent iterations, the most recent estimates of displacement are utilized to recalculate the restoring force \(\mathbf{f}^{(k-1)}(t+\Delta t)\) and the tangent stiffness \(\mathbf{K}(t)\). If the balance in Equation (62) is not achieved, an increment in displacement is required, continuing until equilibrium in Equation (62) is reached.