Speed-Dependent Eigenmodes for Efficient Simulation of Transverse Rotor Vibration

Abstract

Accurate, computationally efficient simulations enable engineers to design high-performing, cost-efficient, lightweight machines that can leverage models of predictive controls and digital twin predictive maintenance schedules. This study demonstrates a new speed-dependent eigenmode method for accurately and efficiently simulating shaft transverse vibrations. The method involves first independently computing shaft eigenmodes over a range of operating speeds, then correlating the eigenmodes across the different speeds during compilation, and finally adjusting modal properties gradually in accordance with a lookup method during simulation. The new method offers several distinct advantages over the traditional static eigenmodes and Craig-Bampton methods. The new method maintains accuracy over a large range of shaft rotation speeds whereas the static eigenmodes method does not. The new method typically requires fewer modal degrees of freedom than the Craig-Bampton method. Whereas the Craig-Bampton method is limited to modeling changes at the boundaries, the new method is suitable for modeling changing body properties as well as boundary-based changes. For this paper, a fluid-bearing-supported 10 MW direct-drive wind turbine drive shaft is tested virtually in a simulation model developed in Simscape™ Driveline™. Using the simulation statistics, this study compares the accuracy and computational efficiency of the speed-dependent eigenmode method to the traditional finite lumped element, static eigenmode, and Craig–Bampton methods. This paper shows that the new method simulates the chosen system 5 times faster than the traditional lumped mass method and 2.4 times faster than the Craig-Bampton method.

1. Introduction

Rotors play crucial energy transfer roles in machines used for energy generation, transportation, and other industrial applications. Since the rotor and its loads are never perfectly symmetrical about the rotor axis, the rotational motion is always accompanied by some transverse vibration. This vibration can contribute to lower quality power output, machine downtime due to failures, and a reduced machine lifetime. If the vibration is excessive, catastrophic failure occurs [1].

1.1. The Role of Rotor Dynamics Modeling

Modeling rotor dynamics throughout the machine design process and deployment can improve machine performance. Designers need to keep rotating machinery vibration within an acceptable tolerance [2]. Modeling the vibration response of rotors with different structural and control parameters can help optimize the machine performance and cost. During the design verification stage for complex machines, simulations of many different load cases help build confidence that a design will operate successfully throughout its lifetime. Simulating the machine response in many scenarios is lower cost and less time-consuming than what is possible using a hardware prototype. Simulation often plays a key role before creating a hardware prototype [3].

The ability to simulate physical systems in real time is a key stage of model-based design for efficiently and robustly developing control systems, software, and hardware. If the plant model can run in real time, then it can be used in place of a hardware prototype to test a processor and production code. Simulation permits testing in conditions that would damage equipment or personnel and before a hardware prototype exists. Moreover, simulation reduces costs and enables non-stop testing. In addition to the development process, real-time simulation can be used in the final product, such as in human-in-the-loop simulators for training, model-predictive controllers, and digital twins for predictive maintenance [4].

1.2. Rotor Dynamics Modeling Challenges

Traditional rotor dynamics models require tradeoffs between simulation fidelity and computational cost, as illustrated in Figure 1. Rotor dynamics behavior is complex. It involves elastic deformation with multiple resonant frequencies; asymmetric, speed-dependent bearing stiffness and damping; and disk gyroscopic effects. These effects can lead to instabilities that make modeling them both challenging and important [5].

Model order reduction methods for rotors enable detailed and accurate simulations at a lower computational cost, which can accelerate and enhance the design process or permit real-time simulation. Rotor model order reduction methods are often used to design shaft structures and controllers, analyze them, and optimize their sensitivity. Key goals for any model reduction method include linear independence, completeness, high accuracy, low computational expense, simplicity, and automation. There are several highly common model reduction methods for rotors, each of which achieve some but not all of these goals [5].

Guyan reduction (or static reduction) partitions the system into primary and secondary degrees of freedom (DOFs). External loads excite the primary DOFs that, in turn, excite the secondary DOFs. This method is accurate only when the stiffness properties dominate the inertia properties [5,6].

Real modal analysis uses eigenanalysis to reduce and decouple the rotor equations of motion. Complex modal analysis that retains damping in the mode computation is not commonly used because it is overly complex and frequency-dependent. Typically, only modes less than twice the excitation frequency are retained. This method can reduce the computational cost of up to tens of thousands of DOFs by a factor of 10–100. However, for rotor systems that include fluid film bearings and other speed-dependent shaft properties, this traditional eigenanalysis can be inaccurate for some shaft rotation speeds, so it often cannot be used for transient system analysis [5].

Component mode synthesis is one of the most common model order reduction methods for rotor dynamics because it can accurately model speed-dependent effects. As with Guyan reduction, component mode synthesis partitions a system into substructures and then further partitions the substructures into primary and secondary DOFs. The primary DOFs interface with the other substructures. Static modes relate the primary DOFs to the secondary DOFs. By using assumed states for the primary DOFs, eigenanalysis reduces the secondary DOFs to inner modal DOFs. To compute the inner modes, some variations of component mode synthesis assume that the primary DOFs are free whereas the popular Craig–Bampton method assumes that they are stationary [6,7].

Component mode synthesis is well-suited for modeling the damping and time-varying properties of rotor supports because the inner mode properties are independent of the supports acting on the interface DOFs. However, in component mode synthesis the interface modes and inner modes are not independent and this DOF dependency complicates conceptual analysis of the system. Furthermore, mixing dynamics that vary both quickly and slowly results in a numerically stiff system that limits the computational efficiency of a model, and in component mode synthesis, bearing support DOFs tend to be stiffer than the modal DOFs [5].

This paper describes a new speed-dependent eigenmode method that maintains accuracy when shaft properties are speed-dependent during varying rotation speed in a simulation. The method involves first choosing several shaft operating speeds and independently computing eigenmodes of the shaft-bearing system at each speed. Then, the method correlates the eigenmodes across the different speeds to generate lookup tables of modal properties that gradually vary with the shaft speed. Finally, during simulation, a lookup method gradually adjusts modal properties as the shaft speed changes. We hypothesize that, for systems with certain characteristics, the new method will have better numerical performance than the Craig–Bampton method because incorporating support properties into the modes eliminates interface DOFs and decreases system numerical stiffness. We expect that the new method will have a simulation performance penalty compared to the traditional static eigenmode method because of the cost of lookup table computations. However, compared to the static eigenmode method, the new method maintains accuracy during transient simulations because it updates parameters that depend on transient variables.

1.3. Paper Organization

The main goals of this work are to describe the new speed-dependent eigenmode method for modeling transverse vibrations and to compare its performance to that of the non-reduced lumped mass model, static eigenmode reduction method, and Craig–Bampton reduction method in a simulation with variable shaft speed. Section 2 summarizes the four modeling approaches. Section 3 describes a wind turbine driveline simulation model implemented in Simscape™ Driveline™ that is used for the model performance comparison [8]. Section 4 describes and discusses the results of the case study. Section 5 summarizes key conclusions and future steps.

2. Modeling Methods

This section describes four methods for modeling the shaft bending vibrations: the lumped mass model, the Craig–Bampton method, the static eigenmodes method, and the new speed-dependent eigenmodes method. The last three methods are model order reduction methods performed on the lumped mass model. The lumped mass model considers shaft bending caused by a force unbalance that excites the shaft at a multiple (pole) of the rotation. The shaft diameter varies along the shaft axis. The shaft is supported by fluid film bearings that exhibit different linearized stiffness and damping properties as the shaft speed changes.

In this study, all four methods are implemented in a Simscape™ model in the Simulink simulation environment [9,10]. The lumped mass model, static eigenmode method, and speed-dependent eigenmode method are implemented in the Simscape Driveline Flexible Shaft block [11]. The Craig–Bampton method is implemented using a custom Simscape component. A link to the models and scripts used to simulate the system is provided in [12].

2.1. Lumped Mass Model

This section summarizes the well-known lumped mass model for shaft bending [2,11]. The model divides the shaft into flexible linear elements with mass concentrated at nodes on either end. Each node has four DOFs: translational and rotational motions in both the x and y directions perpendicular to the shaft axis, z. That is, the DOFs in each node are:

$${\overrightarrow{x}}_i=\left[\begin{array}{c}{x}_i\\ y_i\\ {\theta }_i\\ {\varphi }_i\end{array}\right].$$

(1)

The lumped mass equation of motion is:

$$M\ddot{\overrightarrow{x}}+\left(B+G_{Disk}\Omega \right)\dot{\overrightarrow{x} }+\left(K+G_{Disk}\dot{\Omega }\right)\overrightarrow{x}=\overrightarrow{f},$$

(2)

where $\Omega$ is the shaft speed.

The mass matrix of the ith flexible element, which affects the ith and (i + 1)th nodes, is:

$$M_i=\left[\begin{array}{cccccccc}\frac{m}{2}& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{m}{2}& 0& 0& 0& 0& 0& 0\\ 0& 0& I_d& 0& 0& 0& 0& 0\\ 0& 0& 0& I_d& 0& 0& 0& 0\\ 0& 0& 0& 0& \frac{m}{2}& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{m}{2}& 0& 0\\ 0& 0& 0& 0& 0& 0& I_d& 0\\ 0& 0& 0& 0& 0& 0& 0& I_d\end{array}\right],$$

(3)

where m is the mass of the shaft element and $I_d$ is the diametric mass moment of inertia of half the shaft element about the x and y axes.

The mass matrix of a rigid disk attached to the ith node is modeled as:

$$M_{Disk,i}=\left[\begin{array}{cccc}{m}_{Disk,i}& 0& 0& 0\\ 0& m_{Disk,i}& 0& 0\\ 0& 0& I_{d,Disk,i}& 0\\ 0& 0& 0& I_{d,Disk,i}\end{array}\right],$$

(4)

where $m_{Disk,i}$ is the disk mass and $I_{d,Disk,i}$ is the disk diametric mass moment of inertia about the x and y axes.

The polar mass moment of inertia of a rigid disk attached to the ith node is modeled as:

$$G_{Disk,i}=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& I_{p,Disk,i}\\ 0& 0& I_{p,Disk,i}& 0\end{array}\right].$$

(5)

The stiffness matrix of the ith flexible element is modeled as an Euler-Bernoulli beam:

$$K_i=\frac{2EI}{l^3} \left[\begin{array}{cccccccc}6& 0& 0& 3l& -6& 0& 0& 3l\\ 0& 6& -3l& 0& 0& -6& -3l& 0\\ 0& -3l& 2l^2& 0& 0& 3l& l^2& 0\\ 3l& 0& 0& 2l^2& -3l& 0& 0& l^2\\ -6& 0& 0& -3l& 6& 0& 0& -3l\\ 0& -6& 3l& 0& 0& 6& 3l& 0\\ 0& -3l& l^2& 0& 0& 3l& 2l^2& 0\\ 3l& 0& 0& l^2& -3l& 0& 0& 2l^2\end{array}\right],$$

(6)

where EI is shaft rigidity and l is the element length along the shaft axis. The described model neglects stress stiffening and spin softening due to centrifugal forces, which are negligible because the shaft spins at low speeds. We may incorporate these effects into the models as part of future work.

Fluid film bearings that support the shaft at the ith node are modeled as linearized stiffness and damping coefficients that vary with the shaft speed, $\Omega$. The linearized coefficients are valid for small shaft amplitude motions about the equilibrium position. Details on the linearization of fluid film bearings are provided in [13]. A bearing at the ith node is characterized by a stiffness matrix, K_Support,i, and damping matrix, B_Support,i:

$$K_{Support,i}(\Omega )=\left[\begin{array}{cccc}{K}_{xx}(\Omega )& K_{xy}(\Omega )& 0& 0\\ K_{yx}(\Omega )& K_{yy}(\Omega )& 0& 0\\ 0& 0& K_{\Theta \Theta (\Omega )}& 0\\ 0& 0& 0& K_{\Phi \Phi }(\Omega )\end{array}\right],$$

(7)

$$B_{Support,i}(\Omega )=\left[\begin{array}{cccc}{B}_{xx}(\Omega )& B_{xy}(\Omega )& 0& 0\\ B_{yx}(\Omega )& B_{yy}(\Omega )& 0& 0\\ 0& 0& B_{\Theta \Theta (\Omega )}& 0\\ 0& 0& 0& B_{\Phi \Phi }(\Omega )\end{array}\right].$$

(8)

Additional damping is modeled proportionally to the system mass and stiffness through coefficients α and β, respectively, so that the total system damping is:

$$B=\alpha M+\beta K(\Omega )+B_{Support}(\Omega ).$$

(9)

Damping proportional to stiffness typically represents structural hysteretic damping [14]. In the case study of Section 3, the main contribution to damping in the rotor is generally due to the supports, while the small proportional damping contributes to model numeric stability. The model neglects the internal rotational damping that is typically caused by material stress–strain hysteresis and sliding friction between rotor components because the shaft speed remains well below the whirling mode frequencies [14]. We may add internal rotational damping effects as part of future work.

A static unbalanced magnetic force that acts on the ith node to excite bending in the x and y directions is:

$${\overrightarrow{f}}_i=\left[\begin{array}{c}{f}_i\cos \left(p\omega t+{\xi }_{Offset}\right)\\ f_i\sin \left(p\omega t+{\xi }_{Offset}\right)\end{array}\right],$$

(10)

where $f_i$ is the magnitude of the unbalanced force of a permanent magnetic generator, $p$ is the number magnetic pole pairs, and ${\xi }_{Offset}$ is a phase offset [15,16].

2.2. Static Eigenmode Reduced Model

Standard eigenmode model order reduction reduces the number of simulated DOFs for the lumped mass model of Section 2.1. For rotating shafts, the standard eigenmode method chooses a nominal shaft speed, ${\Omega }_{Nominal}$, for calculating the modal properties. The method neglects the effects of modal properties changing with the shaft speed. Therefore, if the system properties change significantly as the shaft rotation speed changes, then the simulation results are inaccurate at rotation speeds that differ from the nominal speed. The properties of some fluid film bearings change significantly depending on the shaft rotation speed.

For a selected operating speed, ${\Omega }_{Nominal}$, the undamped eigenmodes are computed according to the extended modal reduction method for non-symmetric mass and stiffness matrices [17]. We use a traditional eigenproblem solver implemented in the MATLAB^® eigs function [18,19]. The eigs function takes the form:

$$\left[R_{static}, D_{static}\right]=eigs\left(K\left({\Omega }_{Nominal}\right), M, numModes\right).$$

(11a)

$$\left[L_{static}, D_{static}\right]=eigs\left(K^T\left({\Omega }_{Nominal}\right), M^T, numModes\right)$$

(11b)

The MATLAB eigs function outputs the right eigenvectors in the columns of $R_{static}$, the left eigenvectors in the columns of $L_{static}$, and the corresponding eigenvalues in the diagonal matrix, $D_{static}={\omega }_{static}^2$. The integer input $numModes$ limits the results to the $numModes$ eigenmodes with the lowest frequencies.

Following the traditional eigenmode method, the governing equation for shaft bending, Equation (2), can be reduced from 4 × numNodes DOFs to $numModes$ DOFs. In the case study of Section 3, the eigenmode model order reduction is applied to the lumped mass model with 50 nodes (200 DOFs per Equation (1)), so $K$ and $M$ are matrices with sizes of 200-by-200. The eigenvector matrices $R_{static}$ and $L_{static}$ have sizes of 200-by-$numModes$, where the number of retained modes $numModes$ ranges from 8 to 64, as shown in

Table A1. Performance metrics of the different modeling methods for computing the maximum horizontal force (%ErrX) and vertical force (%ErrY) in the second bearing. # DOF is the number of degrees of freedom in the models: four times the number of lumped elements for the lumped mass model and number of modes for the other three models. T_Exec is the execution time in seconds for a 100 s simulation time. %ErrX and %ErrY are the errors of the predicted forces relative to the lumped mass model simulation with 200 DOFs (50 nodes), the reference result. The bold highlighted results are the simulation for each method with the highest computational efficiency that converges within 2% of the reference result. For further reference, the execution time is 0.08 s when shaft flexibility is disabled.

Lumped Mass Model				Static Eigenmodes Model				Speed-Dependent Eigenmodes Model				Craig–Bampton Model
# DOF	TExec	% ErrX	%ErrY	# DOF	TExec	% ErrX	% ErrY	# DOF	TExec	% ErrX	%ErrY	# DOF	TExec	% ErrX	%ErrY
40	45.6	0.4	−3.0	8	4.4	−20.6	−8.2	8	6.2	15.8	−6.9	20	32.7	−0.0003	−0.13
80	69.2	0.1	−0.5	12	5.0	−14.7	−7.3	12	8.6	9.1	−0.04	24	35.9	−0.04	−0.07
120	130.2	0.05	−0.1	16	5.4	−15.7	−8.2	16	10.2	3.9	−0.7	28	39.4	0.01	−0.21
160	216.2	0.002	−0.2	20	6.1	−16.2	−9.1	20	13.7	1.36	−0.3	32	51.8	−0.005	0.08
200	341.2	0	0	24	7.0	−16.4	−9.5	24	19.1	0.41	−0.2	36	55.1	−0.0006	0.06
				32	8.9	−16.4	−9.8	32	28.1	0.20	−0.4	44	63.9	0.01	0.02
				40	11.3	−16.4	−9.8	40	40.0	0.17	−0.5	52	74.3	−0.006	0.04
				48	14.4	−16.4	−9.8	48	53.1	0.10	−0.06	60	87.5	0.007	0.07
				56	17.9	−16.4	−9.8	56	63.3	0.08	−0.05	68	126.9	0.01	0.01
				64	22.5	−16.4	−9.8	64	85.2	0.06	−0.05

. The reduced equation of motion has the form:

$$M_{Modal}{\ddot{\overrightarrow{x}}}_{Modal}+\left(B_{Modal}+G_{Modal}\Omega \right){\dot{\overrightarrow{x}}}_{Modal}+\left(K_{Modal}+G_{Modal}\dot{\Omega }\right){\overrightarrow{x}}_{Modal}={\overrightarrow{f}}_{Modal},$$

(12)

where

$${\overrightarrow{x}}_{Modal}=R_{static}\overrightarrow{x},$$

(13)

$$M_{Modal}=L_{static}^T M{R}_{static},$$

(14)

$$K_{Modal}=L_{static}^{T}K\left({\Omega }_{Nominal}\right)R_{static},$$

(15)

$$G_{Modal}=L_{static}^T{G}_{Disk}{R}_{static}$$

(16)

$$B_{Modal}=L_{static}^{T}B\left({\Omega }_{Nominal}\right)R_{static},$$

(17)

$${\overrightarrow{f}}_{Modal}=L_{static}^T\overrightarrow{f}.$$

(18)

2.3. Speed-Dependent Eigenmode Reduced Model

In this paper, we consider an extension of the standard eigenmode model order reduction method to increase simulation accuracy while incurring a moderate computational cost for rotational shafts with speed-dependent properties. This extended method involves choosing several reference rotational velocities, ${\Omega }_{Reference}$, for calculating the modal properties. The method correlates modes across neighboring reference rotational velocities by comparing and maximizing correlated mode shape similarity. Then, the method uses the mode shapes to generate modal property lookup tables. As the mode shapes gradually change, the modal properties gradually change. During simulation, as the shaft speed changes, the lookup table interpolation adjusts the modal properties based on the shaft speed.

The speed-dependent eigenmode method is similar to the traditional eigenmode method in that it starts by computing the undamped modal properties. However, the speed-dependent eigenmode method independently computes the modal properties at multiple reference shaft speeds, ${\Omega }_{Reference,i}$, rather than at a single speed. The case study of Section 3 uses reference speeds of [1, 2.5, 5, 7.5, 10] RPM to compute the eigenmodes:

$$\left[R_{SD,i}, D_{SD,i}\right]=eigs\left(K\left({\Omega }_{Reference,i}\right), M, numModes\right).$$

(19a)

$$\left[L_{SD,i}, D_{SD,i}\right]=eigs\left(K^T\left({\Omega }_{Reference,i}\right), M^T, numModes\right)$$

(19b)

The MATLAB eigs function outputs the right eigenvectors in each column of $R_{SD,i}$, the left eigenvectors in each column of $L_{SD,i}$, and the corresponding eigenvalues in $D_{SD,i}={\omega }_{SD,i}^2$.

The method computes the modal properties at each shaft speed i according to:

$${\overrightarrow{x}}_{Modal,i}=R_{SD,i}\overrightarrow{x},$$

(20)

$$M_{Modal,i}=L_{SD,i}^{T}M{R}_{SD,i},$$

(21)

$$K_{Modal,i}=L_{SD,i}^{T}K\left({\Omega }_{Reference,i}\right)R_{SD,i},$$

(22)

$$G_{Modal,i}=L_{SD,i}^T{G}_{Disk}{R}_{SD,i},$$

(23)

$$B_{Modal,i}=L_{SD,i}^{T}B\left({\Omega }_{Reference,i}\right)R_{SD,i},$$

(24)

$${\overrightarrow{f}}_{Modal,i}=L_{SD,i}^T\overrightarrow{f}.$$

(25)

Simscape Driveline uses an automated approach to sort modes at the different reference speeds. The core concept of the sorting method is a modal assurance criterion similar to the ones described in Mogenier et al. [20]. This approach is robust to mode frequency crossing behavior, which is when the modal frequencies of two distinct modes swap order at neighboring reference shaft speeds. Section 3.2 describes a Campbell diagram generated using the speed-dependent eigenmode method and shows that the algorithm captured multiple frequencies crossing each other.

Once the modes at different shaft speeds are correlated with each other and sorted, the corresponding modal properties are stored in three-dimensional tables, where each two-dimensional sheet within the table corresponds to the modal properties at each reference shaft speed. At each time step of the simulation, the 2D modal property matrices are updated based on lookup table linear interpolation of the 3D tables. Therefore, as shaft speed varies during a simulation, the modal properties remain accurate.

The final speed-dependent modal governing equation is:

$${\ddot{\overrightarrow{x}}}_{Modal}+M_{Modal}^{-1}\left(B_{Modal}+G_{Modal}\Omega \right){\dot{\overrightarrow{x}}}_{Modal}+M_{Modal}^{-1}\left(K_{Modal}+G_{Modal}\dot{\Omega }\right){\overrightarrow{x}}_{Modal}=M_{Modal}^{-1}{\overrightarrow{f}}_{Modal},$$

(26)

where the terms are normalized by the modal mass matrix to reduce the number of lookup table computations.

During simulation, the four modal properties, $M_{Modal}^{-1}{B}_{Modal}$, $M_{Modal}^{-1}{G}_{Modal}$, $M_{Modal}^{-1}{K}_{Modal}$, and $M_{Modal}^{-1}{\overrightarrow{f}}_{Modal}$ are computed according to the aforementioned lookup table method:

$$C\left(t\right)=tablelookup\left({\Omega }_{Reference}, C_{Reference}, \Omega \right),$$

(27)

where $C$ represents the four modal properties listed above and $\Omega$ is the transient shaft speed during simulation.

2.4. Craig–Bampton Reduced Model

In this paper, we implement the standard Craig–Bampton reduction method as a performance benchmark for the speed-dependent eigenmode method of Section 2.3 [6,7]. As with the eigenmode methods, we apply the Craig–Bampton method to the shaft lumped mass model described in Section 2.1.

Following the traditional Craig–Bampton method, the partitioning matrix, $T_{part}$, partitions the DOFs into primary nodes, which interface with the supports and external forcing, and secondary nodes. Each primary node contains four DOFs as described in Equation (1). Mathematically, the partitioning process is:

$${\overrightarrow{x}}_{part}=T_{part}\overrightarrow{x},$$

(28)

where

$${\overrightarrow{x}}_{part}=\left[\begin{array}{c}{\overrightarrow{x}}_{Primary}\\ {\overrightarrow{x}}_{Secondary}\end{array}\right].$$

(29)

The partitioned equation of motion matrices are:

$$M_{part}=T_{part}MT{\prime }_{part},$$

(30)

$$K_{part}=T_{part}KT{\prime }_{part},$$

(31)

$$G_{part}=T_{part}GT{\prime }_{part},$$

(32)

$$B_{part}=T_{part}BT{\prime }_{part},$$

(33)

$${\overrightarrow{f}}_{part}=T_{part}\overrightarrow{f}.$$

(34)

The partitioned stiffness matrix has quadrants corresponding to the primary and secondary DOFs:

$$K_{part}=\left[\begin{array}{cc}{K}_{pp}& K_{ps}\\ K_{sp}& K_{ss}\end{array}\right].$$

(35)

The static reduction matrix is obtained from the partitioned stiffness matrix:

$$S_{CB}=\left[\begin{array}{c}{I}_{pp}\\ -inv\left(K_{ss}\right)K_{sp}\end{array}\right].$$

(36)

The columns of the static reduction matrix are the Craig–Bampton static modes. In each static mode, one primary (interface) DOF undergoes a unit displacement. The condensation matrix, $-inv\left(K_{ss}\right)K_{sp}$, relates the secondary (inner) DOFs to each primary DOF.

In addition to the static modes, the Craig–Bampton method computes the eigenmodes of the secondary (inner) DOFs. This study uses the MATLAB eigs function for the eigenmode computation according to:

$$\left[V_s, ~\right]=eigs\left(K_{ss}, M_{ss}, numInnerModes\right),$$

(37)

where $numInnerModes$ is an integer limit to the number of retained inner modes. The inner modes have stationary boundaries that correspond to the primary DOFs. The inner modes are the columns of the dynamic transformation matrix:

$$D_{CB}=\left[\begin{array}{c}{0}_{pp}\\ V_s\end{array}\right]$$

(38)

The final Craig–Bampton transformation matrix is:

$$T_{CB}=\left[S_{CB} D_{CB}\right]$$

(39)

The final Craig–Bampton modal properties are:

$$M_{CB}=T_{CB}{M}_{part}T{\prime }_{CB},$$

(40)

$$K_{CB}=T_{CB}{K}_{part}T{\prime }_{CB},$$

(41)

$$G_{CB}=T_{CB}{G}_{part}T{\prime }_{CB},$$

(42)

$$B_{CB}=T_{CB}{B}_{part}T{\prime }_{CB},$$

(43)

$${\overrightarrow{f}}_{CB}=T_{CB}\overrightarrow{f}.$$

(44)

The final Craig–Bampton governing equation is:

$${\ddot{\overrightarrow{x}}}_{CB}+M_{CB}^{-1}\left(B_{CB}+G_{CB}\Omega \right){\dot{\overrightarrow{x}}}_{CB}+M_{CB}^{-1}\left(K_{CB}+G_{CB}\dot{\Omega }\right){\overrightarrow{x}}_{CB}=M_{CB}^{-1}{\overrightarrow{f}}_{CB}.$$

(45)

The system forcing accounts for both the static unbalanced magnetic force given by Equation (10) and the support forcing. That is, $K_{CB}$ and $B_{CB}$ do not account for the support loads. The support loads acting at the primary nodes i are:

$${\overrightarrow{f}}_i=-K_{Support,i}(\Omega ){\overrightarrow{x}}_{Support,i}-B_{Support,i}(\Omega ) {\dot{\overrightarrow{x}}}_{Support,i},$$

(46)

where lookup table interpolation computes the support stiffness and damping at the varying shaft speed, according to:

$$K_{Support,i}\left(t\right)=tablelookup\left({\Omega }_{Reference}, K_{Support,i}{}_{Reference}, \Omega \right).$$

(47)

The Craig–Bampton transformation reduces the original 4 × numNodes DOFs into 4P + S DOFs, where P is the number of interface (primary) nodes and S is the number of retained inner (secondary) modes. The number of inner modes required for simulation convergence generally depends on the system stiffness, excitation frequency, and location of transverse excitation. Systems with weaker shaft stiffnesses, higher excitation frequencies, and transverse excitations located at mode shape displacement peaks typically require more inner modes.

3. Case Study: Wind Turbine Driveline Model

To benchmark the different modeling methods, a fluid-bearing-supported 10 MW direct-drive wind turbine drive shaft with transverse vibrations is tested virtually in a simulation model developed in Simscape Driveline. The wind industry is generally moving towards using direct-drive generators in large offshore wind turbines because direct-drive generators likely offer better availability and efficiency and have reduced operation and maintenance requirements compared to geared wind turbines [21].

In wind turbines, rotor vibration contributes to main bearing failure, the rate of which can be as high as 30% over a 20-year wind turbine lifetime. Reliability of the main bearings is often considered one of the most important mechanical reliability challenges in wind turbines, after gearbox reliability [22].

Wind turbine driveshafts are subject to complex sources of transverse vibration. Although rotors spin at low-rated speeds ranging from about 8–15 rpm, there are many sources of higher harmonic transverse vibration excitation. As blades pass over the tower, their frequency causes a third harmonic excitation. Gearbox and generator unbalanced magnetic forces can generate excitations ranging from approximately 80–2000 times the rotor speed [23]. In this paper, the source of transverse vibration is a linearized unbalanced magnetic force corresponding to a generator with 200 poles [24,25].

Fluid film bearings offer several advantages over ball or roller bearings. Most importantly, when designed properly, they can have very low friction and nearly infinite life. Tolerance-related manufacturing issues due to increasing wind turbine size negatively impact the reliability of ball and roller bearing systems. Reliability tends to also be negatively affected by larger drive shaft diameters as the required number of rolling elements increases. While repairs for traditional main bearings often require some disassembly of the rotor, it is possible to design plain journal bearings or other fluid film bearings to avoid disassembly for repairs. Additionally, when designed properly, fluid film bearings can dampen vibration [26,27].

From a modeling standpoint, the speed-dependent properties of fluid film bearings make them more challenging to model. While this paper considers the modeling challenges of rotors supported by plain journal bearings, ball and roller bearings with lubrication face similar modeling challenges because they also have speed-dependent stiffness and damping properties [28].

3.1. System Concept

Figure 2 shows a schematic of the wind turbine and

Table 1. 10 MW direct drive wind turbine parameters. These parameters are based on [24,31,32].

Parameter	Value
Wind turbine electric power rating	10 MW
Wind turbine power efficiency	94%
Wind turbine radius	99 m
$\mathrm{Mass} \mathrm{of} \mathrm{blades}, m_{Blades}$	228 tons
$\mathrm{Diametric} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{blades}, I_{D,Blades}$	3.4 × 107 kg m2
$\mathrm{Polar} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{blades}, I_{P,Blades}$	6.8 × 107 kg·m2
$\mathrm{Mass} \mathrm{of} \mathrm{hub}, m_{Hub}$	106 tons
$\mathrm{Diametric} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{hub}, I_{D,Hub}$	9.29 × 106 kg·m2
$\mathrm{Polar} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{hub}, I_{P,Hub}$	1.68 × 106 kg·m2
$\mathrm{Generator} \mathrm{axal} \mathrm{location}, z_{Generator}$	6 m
$\mathrm{Mass} \mathrm{of} \mathrm{generator} \mathrm{rotor}, m_{Generator}$	137 tons
$\mathrm{Diametric} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{generator} \mathrm{rotor}, I_{D,Generator}$	8.5 × 106 kg·m2
$\mathrm{Polar} \mathrm{moment} \mathrm{of} \mathrm{inertia} \mathrm{of} \mathrm{generator} \mathrm{rotor}, I_{P,Generator}$	5.8 × 106 kg·m2
Rated shaft speed	8.6 RPM
Number of poles of generator, p	200
Generator maximum static unbalance magnitude, f	20 kN
Shaft length	10 m
Shaft maximum outer diameter	1.5 m
Shaft minimum outer diameter	1.25 m
Shaft maximum radial thickness	0.625 m
Shaft minimum radial thickness	0.375 m
Shaft density	7800 kg/m3
Shaft elastic modulus	200 GPa
System damping proportional to mass, α	10
System damping proportional to stiffness, β	1 × 10−9
$\mathrm{Bearing} \mathrm{axial} \mathrm{locations}, z_{Support}$	[0, 5, 8, 10] m
Bearing fluid density	1077 kg/m3
Bearing fluid viscosity	3.2 × 10−4 m2/s
Bearing axial length	30 cm
Bearing radius	0.6875 m
Bearing clearance	0.4 mm

lists the parameters. The wind turbine rotor model consists of three 99 m rigid blades mounted on a hub. We assume that the blades are rigid since the goal of the paper is to test the computational performance of the different shaft models while accounting for the large inertia effects of the rotor. Future work may involve modeling flexible blades in addition to the flexible shaft. The flexible shaft has a length of 10 m. Its front end has a diameter of 1.5 m that linearly decreases to 1.25 m at its rear end. Four fluid film bearings support the shaft at axial positions of 0, 5, 8, and 10 m from the front end. This model assumes that the bearing casings are stationary. A generator with 200 poles is located at the shaft midpoint. Unbalanced magnetic pull due to rotor sag causes a transverse force with a maximum amplitude of 20 kN.

Figure 3 shows the wind turbine simulation model developed in Simulink and Simscape Driveline. This model is available at [12]. Simscape blocks model the torque generated by wind acting on the rotor, the rotor and hub inertia, the shaft with transverse vibrations, the generator torque, and a disc brake. The blade torque block computes the wind turbine torque using a lookup table of power coefficient versus tip-speed ratio and blade pitch. Blade pitch is modeled quasi-statically, responding instantaneously to the pitch controller command.

The controllers in the model are based on the well-known National Renewable Energy Laboratory 5 MW Reference offshore wind turbine controller scaled up to 10 MW [29,30]. Generator torque and blade pitch controllers are modeled using Simulink blocks while a simple brake controller is modeled using Simscape blocks.

The generator controller uses shaft speed and wind speed as inputs. It has four main control regimes: a start-up region when no torque is applied, an optimized control region that tries to maintain an optimal tip-speed ratio for maximum power, a rated power region when the generator provides the rated power, and a cut-out region where the generator turns off during wind speeds that exceed 25 m/s.

When the wind turbine is in the rated power region, a blade pitch controller adjusts the collective blade pitch to maintain the rated turbine speed. The blade pitch controller uses scheduled proportional and integral gains to reduce the shaft speed error. When the wind speed exceeds 25 m/s, the blades are feathered to 90 degrees and the disc brake is engaged to avoid overloading the generator.

Figure 4 shows the mechanical power, collective blade pitch, and shaft rotation speed in response to the simple wind speed scenario used in the case study of Section 3.3.

3.2. Rotordynamics Characteristics

Figure 5 shows the bearing properties versus rotation speed. These stiffness and damping properties represent a bearing with a radius of 0.6875 m, clearance of 0.4 mm, and Sommerfeld numbers over the relatively low range of 2 × 10⁻³ to 2 × 10⁻². Low Sommerfeld numbers are typical for bearings with low speed and high radial load. When generating these bearing properties, we decrease the bearing damping compared to a typical journal bearing because of the current limitation that Simscape Driveline computes only undamped eigenmodes. Decreasing the support damping permits closer convergence between the lumped mass model and the speed-dependent eigenmode model. While computation of the mode shapes assumes undamped modes (Equations (11) and (19)), the modal properties account for support damping (Equations (17) and (24)). Future work will consider extending the model to account for complex mode shapes. As is typical for fluid film bearings, the vertical stiffness, subject to the shaft gravitational load, has the higher stiffness. For more detailed information, see [13].

Figure 6 shows the Campbell diagram for the 20 lowest eigenmodes of the system generated using the speed-dependent eigenmode reduced model. To avoid static unbalance instabilities, none of the whirling frequencies match the shaft speed. Several of the modes pass through the 200th harmonic of the shaft rotation speed, which corresponds to the generator magnetic unbalance. These include modes 1 and 2 at a shaft speed of 0.5 RPM, mode 4 at a shaft speed of 5.75 RPM, and mode 5 at a shaft speed of 6.5 RPM. Therefore, to avoid excess vibration in this system, it is best for the controller to quickly pass through these shaft rotation speeds. The frequencies of several modes cross each other when the shaft speed passes through 2 RPM. To produce an accurate response prediction, the Simscape Driveline algorithm for sorting modes during compile time that is described in Section 2.3 can capture this behavior.

Figure 7 shows the shapes of the eight lowest frequency eigenmodes at the rated shaft speed of 8.6 RPM. The mode shapes do not significantly change with rotation speed even though some mode frequencies do because the bearings are relatively stiff. Modes 1, 3, 5, and 6 have large amplitude in the X direction and small amplitudes in the Y direction, while modes 2, 4, 7, and 8 have large amplitudes in the Y direction and small amplitudes in the X direction. These modes are relatively planar because the bearing cross-coupling stiffnesses are relatively small.

3.3. Case Study Scenario

The case study uses the wind speed scenario shown in Figure 4, The scenario consists of the wind ramping up from 0 m/s to 20 m/s over 30 s, holding steady at 20 m/s for 30 s, ramping back down to 0 m/s over 30 s, and then remaining at 0 m/s for the final 10 s. In response to this wind speed, the wind turbine shaft accelerates from rest until it reaches a peak speed of 10.9 RPM and then gradually decelerates to maintain a constant power.

To fully exercise the lookup-table-based algorithm, we choose this scenario, wherein the simulation passes through a large range of shaft speeds. To ensure that the measured simulation execution times are long enough that compilation start-up effects are not significant, and that data can still be gathered relatively quickly, we choose a simulated time of 100 s. A more typical wind turbine simulation might span an hour or more of simulated time and have more wind speed and shaft speed fluctuations. Additionally, as mentioned in the introduction, a typical wind turbine simulation will account for different sources of shaft bending excitation with different harmonics. This simplified scenario permits the simulation performance to be dominated by the vibration model performance.

3.4. Numeric Methodology

The simulation data are from an Intel^® Xeon^® W-2133 CPU @ 3.60GHz x64-based processor. For each modeling method and chosen number of DOFs, a MATLAB script runs the simulation multiple times and averages the execution time.

The Simulink model uses the ode23t solver with an absolute tolerance of 1 × 10⁻³. This model uses variable nominal value scaling to improve the simulation speed for each modeling method. Nominal values scale the variables seen by the solver. When the solver sees variables with scaled magnitudes on the order of 1, the simulation can have improved robustness and execution speed. We use the Simscape Results Explorer and Simscape Variable Scaling Analyzer to determine the magnitudes of each shaft DOF and its time-derivative. For each modeling method, we set all of the shaft DOFs to have nominal values equal to the peak value of the DOF with the largest observed response.

4. Results and Discussion

This section compares the performance of the four different methods for simulating the wind turbine shaft transverse vibrations. These results correspond to the wind speed and shaft speed scenario shown in Figure 4. Figure 8 plots the percent error in the maximum horizontal and vertical forces acting on the second bearing versus the run/sim time ratio. The run/sim time ratio is the ratio of the execution time to the simulated time of 100 s. Figure 8 is based on simulations with transient shaft speeds. Figure 9 plots a snapshot of the transient time series for each simulation method. The performance trends for the chosen variables are representative of the performance trends for other variables.

Along each curve in Figure 8, the markers represent simulations with different numbers of DOFs. Each DOF corresponds to four times the number of nodes for the lumped mass method and the number of modes for the order reduction methods. Each DOF is associated with a displacement variable and a velocity variable. Simulations with fewer DOFs have shorter execution times and smaller run/sim time ratios. Table A1 in the Appendix A provides details on the number of DOFs in each simulation.

The circle on each curve indicates the best performance for each simulation method. For each method, we define the result with the best performance as the most computationally efficient simulation that converges within 2 percent of the reference result. We use the results of the lumped mass model with 200 DOFs (50 nodes) as the reference result. The reference results are 2858 N for the peak horizontal force and 3175 N for the peak vertical force. The bottom subplots in Figure 8 summarize the simulations with the best performance for each method. Table A1 in Appendix A provides more details on the results shown in Figure 8.

The speed-dependent eigenmode method has the best performance in terms of accuracy and computational efficiency. Considering the best performance for each simulation method, the speed-dependent eigenmode method simulates 5-times faster than the lumped mass method and 2.4-times faster than the Craig–Bampton method.

While the static eigenmode method is the most computationally efficient simulation, it does not converge to an accurate result, regardless of how many eigenmodes are added. The static eigenmode properties are based on a nominal shaft speed of 4.3 RPM, so the response of the model is inaccurate at the instants of peak bearing force, which occur at shaft speeds greater than 7 RPM. We expect accuracy errors for the static eigenmode method because the simulation has varying shaft speeds and the model properties are sensitive to the shaft speed. Therefore, these static eigenmode results serve as a lower-bound reference on the run/sim time ratio of the other methods rather than providing a meaningful % error comparison.

The Craig–Bampton method generally requires more DOFs than the speed-dependent eigenmode method because the Craig–Bampton method requires four static modes for each external load and support interface in addition to the inner modes that capture flexible dynamics. Since the system in this case study has four support nodes and one forced node, the Craig–Bampton system requires 20 static modes (at least 20 DOFs). The eigenmode approaches capture external load, support, and flexibility dynamics all within each mode. The Craig–Bampton static modes, which are coupled with the inner modes, tend to have different stiffnesses than the inner modes, and this makes it a more challenging problem for the numeric solver. An advantage of the Craig–Bampton approach for this system is that it requires few inner modes to converge to the accurate solution. More Craig–Bampton inner modes are required when a shaft is more flexible.

The speed-dependent eigenmode method requires a sufficient number of modes to converge. Some of the modes have natural frequencies that are much higher than the excitation frequency, but the small excitation of these modes still has a non-negligible effect on the force since higher modes with higher bending tend to impart higher forces on the bearings.

For all modeling methods, the vertical force requires more DOFs than the horizontal force to converge to the reference value. This is related to how the system has higher stiffness in the vertical direction.

Previous practical rotating shaft model order reduction methods have focused on component mode synthesis approaches such as the Craig–Bampton method. The hypothesis of this study was that the additional cost of lookup table interpolation would be more than offset by the computational efficiency of the new method because it generally requires fewer DOFs (none of which are coupled) than the Craig–Bampton method. The results of the wind turbine driveline case study confirm this hypothesis. As mentioned in the introduction, detailed and computationally efficient rotor dynamics models can accelerate and strengthen conceptual design exploration, parameter optimization, and system design verification and validation. Detailed, computationally efficient models can also enable real-time simulation for controller hardware testing and predictive maintenance, which can reduce the risk of failures and improve machine performance. These computational benefits become more significant when shaft models are based on higher-cost finite element analysis approaches, and simulations require longer simulated scenarios.

In addition to more computationally efficient simulation, the uncoupled DOFs of the speed-dependent eigenmode method make it easier to fundamentally understand the system behavior and make design decisions.

The speed-dependent eigenmode approach has a few limitations. The method considers the entire shaft-bearing-disk system when computing the eigenmodes. This approach is more restrictive than the Craig–Bampton method for interfacing with other machine subcomponents. Additionally, the speed-dependent eigenmode method in Simscape Driveline currently restricts the bearing casing to be stationary. A limitation of this study is that the eigenmode methods consider only undamped vibration mode shapes. In order to ensure that the eigenmode results could sufficiently converge with the reference lumped mass method, we included only light damping in the system. Future extensions could loosen these restrictions to permit additional external excitations and parameters of the system.

In the future, this approach could be expanded to complex modes and additional gradually changing variables and systems. Unlike the Craig–Bampton approach, the speed-dependent eigenmode approach is well-suited for capturing the behavior of changing body characteristics in addition to changing boundary characteristics.

5. Conclusions

In this paper, we show that a speed-dependent eigenmode method is capable of accurately simulating the transient transverse vibrations of a wind turbine driveline with a run to simulation time ratio of 0.14 for excitation frequencies that range from 0 Hz to 28 Hz. The speed-dependent eigenmode method simulates 2.4 times faster than the Craig–Bampton method. A limitation of this Craig–Bampton method performance is that it requires four DOFs for every interface. A limitation of the speed-dependent eigenmode method performance is that it requires lookup table computations for the varying modal properties during every simulation time step. Accordingly, we observe that the speed-dependent eigenmode method computational performance is generally better than that of the Craig–Bampton method when characterizing the system behavior requires fewer eigenmodes and more Craig–Bampton interface nodes for an accurate simulation. For future work in this space, we will consider extending the method to model complex modes and more detailed finite element analysis models in the Simscape Driveline simulation model. The authors encourage readers to adopt the model at [12] for their own studies and contact J.K. with any questions.