Variation patterns of modal characteristics in a CAES multi-shaft integrally geared compressor rotor system under variable operating conditions

Abstract

The modal characteristics of the multi-rotor coupling system in an integrally geared compressor (IGC) for Compressed Air Energy Storage (CAES) have a critical influence on system dynamics and reliability under variable operating conditions. Taking into account the coupling effects of multi-gear meshing, a mathematical model of the IGC rotor system was established. Firstly, the modal characteristics of each independent shaft train were calculated before and after the gear meshing. Analysis revealed that gear engagement alters the natural frequencies and mode shapes of individual shaft trains and introduces new gear-related frequency components. Subsequently, an orthogonal experiment was designed for the control parameters (power, rotational speed, inlet oil temperature), followed by range analysis and analysis of variance (ANOVA). This process revealed the effects of each control parameter on the various modal characteristics and the derived natural frequencies of the rotor system. Rotational speed exhibited the highest sensitivity, followed by power and then inlet oil temperature. Multivariate linear regression prediction models were then established for each dependent variable, yielding the regression coefficient values for each control parameter. Finally, validation using field vibration data confirmed the accuracy of the method and models, providing a theoretical basis for optimizing CAES control strategies.

Keywords

variable operating conditions integrally geared compressor coupled rotor system control parameters modal characteristics CAES

Introduction

With the continuous increase in the proportion of renewable energy installed capacity worldwide, wind and photovoltaic (PV) power generation are gradually becoming the mainstay of new power systems. This trend, however, brings numerous new challenges to power systems, including power balance, energy accommodation, and grid security and stability control.^1,2 To address the inherent randomness, volatility, and intermittency of renewable power generation, energy storage technology has emerged as a crucial supporting technology for building new-type power systems and enhancing power ancillary service capabilities.

CAES is an emerging long-duration, large-scale energy storage technology. It offers significant advantages such as large storage capacity, short construction period, long operational lifespan, and environmental compatibility, and has been developing rapidly in major electricity-consuming countries in recent years. CAES holds broad application prospects in areas including peak shaving and valley filling, peak regulation, renewable energy accommodation, and ancillary services. It is of great importance for promoting the development of new-type power systems in the context of “Carbon Peak and Carbon Neutrality”.^3–5

CAES technology is based on a modified gas-turbine cycle in which compression and expansion are separate processes. In the charging phase, excess electricity is used to compress air, which is then stored; the heat of compression is concurrently stored in a Thermal Energy Storage (TES) system. In the discharging phase, the stored air is released and heated before being expanded through a turbine to generate electricity.^6,7 The two commercial CAES facilities—the Huntorf plant in Germany (commissioned in 1969) and the McIntosh plant in the United States (commissioned in 1991)—employ underground storage with capacities of 532,000 m³ and 270,000 m³, respectively. The potential of CAES to compensate for fluctuations in renewable energy was first mentioned as early as 1976 and again in 1981, although it did not receive significant attention at the time. However, this situation has changed dramatically with the widespread integration of intermittent renewable energy sources, such as wind and solar PV, into power supply systems across many countries.^8–10 In recent years, CAES research has primarily focused on thermodynamic cycles, design and analysis of key components, and system integration and validation. The main components of a CAES system include the motor, compressor, heat exchanger, storage vessel, expander, and generator.^11–13

The turbocompressor is one of the core components in a CAES system, and its performance directly determines the overall efficiency and economy of the system. The pressure ratio required in a CAES system is typically as high as 40–80 or above, making multistage compression a common configuration. Multistage centrifugal compressors are mainly divided into single-shaft and multi-shaft types. Among these, the multi-shaft geared centrifugal compressor offers advantages such as high efficiency, compact structure, and ease of inter-stage heat recovery, making it more suitable for large-scale CAES systems.¹⁴ Efficient energy conversion is central to energy storage technology. Therefore, enhancing the operational stability of the compressor under variable working conditions is one of the key technical challenges for CAES systems.

In terms of key technologies for IGCs, scholars have conducted multidimensional research. Wang et al.¹⁵ considered the gear tooth surface friction effect, calculated tooth profile deformation caused by contact temperature variation, and investigated the influence patterns of rotational speed, torque load, and lubricant viscosity on system response. Li et al.¹⁶ proposed a control strategy for the grid-connection process of an Advanced Adiabatic Compressed Air Energy Storage (AA-CAES) system, achieving efficient and stable operation of the energy storage system by regulating parameters such as the power and speed of the turbocompressor and expander. Song et al.¹⁷ investigated the mechanical and electrical run-out errors at key journal locations of the IGC rotor. They demonstrated that a combined approach involving demagnetization, cold extrusion, and micro-peening can effectively reduce mechanical-electrical run-out errors. To mitigate inducer stall and extend the low-flow operating range, Pelton et al.¹⁸ proposed an open impeller design with a passive recirculation casing treatment for wide-range centrifugal compressor stages in supercritical CO₂ (sCO₂) power cycles. Modekurti et al.¹⁹ proposed an integrated model for post-combustion CO₂ compression using an integrally geared multistage centrifugal compressor. Using dimensionless performance curves, they simulated its off-design performance and analyzed the response of key variables to disturbances. Zhang et al.²⁰ developed a dynamic model of a gear-rotor-bearing system, including bearing dynamics, axial force, and torque. They analyzed the dynamic response using a five-shaft transmission system as a case study. Yang et al.²¹ introduced a general method for calculating the maximum shaft radius in the unbalanced response orbit of geared rotor systems, providing mathematical proof and comparing it with classical modal synthesis. Zhang et al.²² established a finite element model of a centrifugal compressor rotor using a 6-DOF rotating beam element. They investigated how changes in bearing stiffness under load affect critical speeds and resonance. Wang et al.²³ studied the effect of CAES variable conditions on a three-gear rotor system. Their finite element model, incorporating unbalance and torque excitation, was validated experimentally. Yang et al.²⁴ investigated the internal loss mechanisms in a supercritical CO₂ centrifugal compressor, comparing the effects of three working fluids—sCO₂, ideal CO₂ (ICO₂), and ideal air—on compressor performance and loss distribution. Lin et al.²⁵ conducted a numerical study on the full mainstream channel coupled with the impeller backside cavity (IBC) in a centrifugal compressor for CAES. They analyzed the internal flow field within the IBC and the compressor coupling characteristics under variable operating conditions. Wang et al.²⁶ proposed a self-recirculation casing treatment implemented at the diffuser to extend the stability margin. This method defines the recirculation angle and establishes physical parameters to characterize improvements in stall margin (η_stall) and stable flow range, effectively mitigating compressor stall.

As indicated by the above literature, research on IGCs has primarily focused on the internal flow characteristics under single operating conditions or the dynamic response of individual gears, whereas few studies have addressed the modal characteristics of the coupled multi-rotor system. Regarding the parameters incorporated into mathematical models, the effects of mesh stiffness and damping of multiple gear pairs are rarely considered. Variable operating conditions are a typical feature of multi-shaft IGCs used in CAES systems. As the external load and rotational speed continuously vary, the meshing states of the gear pairs, as well as the stiffness and damping parameters of the bearings, undergo changes, leading to dynamic variations in the modal characteristics of the entire rotor system.

Due to the multi-shaft configuration of IGCs, where individual rotors often operate above their first or second critical speeds and multiple gear meshing zones are present, their dynamic behavior and modal characteristics are exceptionally complex. Consequently, it is essential to study the rotor system, which comprises impellers, gears, shafts, and bearings, as an integrated whole. This necessity is even more pronounced in the context of CAES applications under variable operating conditions, where the modal characteristics and their variation patterns in IGC rotor systems become particularly intricate and volatile, posing significant constraints on the unit’s operational stability and reliability.

Mathematical model

Multi-shaft rotor system configuration

The process flow of the CAES system is shown in Figure 1, in which the core device is the multi-shaft integrally geared compressor (IGC). The configuration of its rotor system is illustrated in Figure 2. It consists of a bull gear component (BGC), a first- and second-stage high-speed shaft component (H1), a third-stage high-speed shaft component (H2), and a fourth- and fifth-stage high-speed shaft component (H3). The bull gear simultaneously drives the three high-speed shafts, which in turn rotate the corresponding impeller stages to achieve multistage compression of air, thereby delivering the optimum flow rate and pressure required for energy storage.

Figure 1.

Flow diagram of the CAES system.

Figure 2.

Layout of the three high-speed shaft five-stage compression rotor system.

Equations of motion

By treating the bearing housings or foundation as rigid, the rotor system, which comprises bearings, gears, transmission shafts, impellers, and other components, can be considered an independent system. Within this system, the individual rotors are coupled through gear-pair meshing, which consequently influences the modal characteristics of the entire rotor system. In the dynamic modeling of a gear meshing unit, the primary considerations are the mesh stiffness excitation and displacement excitation, which interact with each other.²⁷ The following assumptions are made: (1) A plane cross-section that is perpendicular to the neutral axis before deformation remains plane after deformation, but is not necessarily perpendicular to the deformed neutral axis. (2) The material is linearly elastic and isotropic, the Poisson effect is neglected, and shear deformation and rotatory inertia effects are incorporated.²⁸ Utilizing the Timoshenko beam element,²⁹ the developed dynamic model of the rotor system is depicted in Figure 3, where K and C represent the support stiffness and damping of the bearings or gears, respectively, and e and 2b denote the transmission error and backlash of the gear pair, respectively.

Figure 3.

Dynamic model of the rotor system.

The time-varying meshing stiffness of the gear pair is expressed by the following equation:

k (t) = k_{0} + \sum_{i = 1}^{\infty} k_{i} \cos (i ω t + φ_{i})

(1)

where k₀ is the mean meshing stiffness, k_i denotes the Fourier expansion coefficients, ω represents the gear meshing frequency, and φ_i stands for the phase angle.

To simplify the calculation, only the first-order term of the time-varying meshing stiffness expansion is retained. Equation (1) can then be expressed as follows:

k (t) = k_{0} [1 + ε \sin (ω t + φ)]

(2)

where ε is the stiffness fluctuation coefficient, and φ is the initial phase angle.

Similarly, by performing a Fourier expansion of the gear-pair transmission error function and retaining only the first-order term, an approximate expression for the gear pair transmission error function can be obtained:

e (t) = e_{s} \sin (ω t + φ)

(3)

where e_s is the amplitude of transmission error fluctuation.

Displacement function along the gear meshing line:

P (X_{i j}) = {BX}_{i j} - e (t)

(4)

where B is a constant vector related to the gear meshing parameters, and X is the generalized displacement vector (see Equation (8)).

The relative deformation along the gear meshing line is expressed as follows:

g (x) = {\begin{cases} p (x) - b p (x) > b \\ 0 | p (x) | \leq b \\ p (x) + b p (x) < b \end{cases}

(5)

where b is half of the gear backlash value.

Based on the time-varying meshing stiffness, transmission error, and displacement function along the meshing line, a multi-degree-of-freedom nonlinear mathematical model that includes the meshing state of the gear pair can be derived using the lumped mass method and Newton’s law. The governing differential equations for each high-speed gear are given as follows:

\begin{array}{l} m_{i} \frac{d^{2} x_{i}}{d t} + c_{i j} \frac{d x_{i}}{d t} + k_{i j} (t) g_{i j} (t) \cos β_{i j} \cos ψ_{i j} = 0 \\ m_{i} \frac{d^{2} y_{i}}{d t} + c_{i j} \frac{d y_{i}}{d t} + k_{i j} (t) g_{i j} (t) \cos β_{i j} \sin ψ_{i j} = 0 \\ m_{i} \frac{d^{2} z_{i}}{d t} + c_{i j} \frac{d z_{i}}{d t} + k_{i j} (t) g_{i j} (t) \sin β_{i j} = 0 \\ {I_{d}}_{i} \frac{d^{2} θ_{x i}}{d t} + r_{i} c_{i j} \frac{d θ_{x i}}{d t} + r_{i} k_{i j} (t) g_{i j} (t) \sin β_{i j} \cos ψ_{i j} = 0 \\ {I_{d}}_{i} \frac{d^{2} θ_{y i}}{d t} + r_{i} c_{i j} \frac{d θ_{y i}}{d t} + r_{i} k_{i j} (t) g_{i j} (t) \sin β_{i j} \sin ψ_{i j} = 0 \\ {I_{p}}_{i} \frac{d^{2} θ_{z i}}{d t} + r_{i} c_{i j} \frac{d θ_{y i}}{d t} + r_{i} k_{i j} (t) g_{i j} (t) \cos β_{i j} = T_{i} \end{array}

(6)

The bull gear engages simultaneously in three meshing zones. Considering the vector sum of interactions, the comprehensive equations of motion are given as follows:

\begin{array}{l} m_{j} \frac{d^{2} x_{j}}{d t} + \sum c_{i j} \frac{d x_{j}}{d t} + \sum k_{i j} (t) g_{i j} (t) \cos β_{i j} \cos ψ_{i j} = 0 \\ m_{j} \frac{d^{2} y_{j}}{d t} + \sum c_{i j} \frac{d y_{j}}{d t} + \sum k_{i j} (t) g_{i j} (t) \cos β_{i j} \sin ψ_{i j} = 0 \\ m_{j} \frac{d^{2} z_{j}}{d t} + \sum c_{i j} \frac{d y_{j}}{d t} + \sum k_{i j} (t) g_{i j} (t) \sin β_{i j} = 0 \\ I_{d j} \frac{d^{2} θ_{x j}}{d t} + \sum r_{j} c_{i j} \frac{d θ_{x j}}{d t} + \sum r_{j} k_{i j} (t) g_{i j} (t) \sin β_{i j} \cos ψ_{i j} = 0 \\ I_{d j} \frac{d^{2} θ_{y j}}{d t} + \sum r_{j} c_{i j} \frac{d θ_{y j}}{d t} + \sum r_{j} k_{i j} (t) g_{i j} (t) \sin β_{i j} \sin ψ_{i j} = 0 \\ I_{p j} \frac{d^{2} θ_{z j}}{d t} + \sum r_{j} c_{i j} \frac{d θ_{z j}}{d t} + \sum r_{j} k_{i j} (t) g_{i j} (t) \cos β_{i j} = T_{j} \end{array}

(7)

Let:

\begin{array}{l} X_{i j} = {[x_{i}, y_{i}, z_{i}, θ_{x i}, θ_{y i}, θ_{z i}, x_{j}, y_{j}, z_{j}, θ_{x j}, θ_{y j}, θ_{z j}]}^{T} \\ M_{i j} = d i a g ([m_{i}, m_{i}, m_{i}, I_{d i}, I_{d i}, I_{p i}, m_{j}, m_{j}, m_{j}, I_{d j}, I_{d j}, I_{p j}]) \\ B = [\cos β_{i j} \cos ψ_{i j}, \cos β_{i j} \sin ψ_{i j}, \sin β_{i j}, r_{i} \sin β_{i j} \cos ψ_{i j}, r_{i} \sin β_{i j} \sin ψ_{i j}, r_{i} \cos β_{i j}, \\ - \sum \cos β_{i j} \cos ψ_{i j}, - \sum \cos β_{i j} \sin ψ_{i j}, - \sum \sin β_{i j}, - \sum r_{i} \sin β_{i j} \cos ψ_{i j}, - \sum r_{i} \sin β_{i j} \sin ψ_{i j}, - \sum r_{i} \cos β_{i j}] \\ C_{i j} = c_{i j} B_{i j}^{T} \cdot B_{i j} \\ K_{i j} = k_{i j} B_{i j}^{T} \cdot B_{i j} \\ F_{i j} = [0, 0, 0, 0, 0, T_{i}, 0, 0, 0, 0, 0, T_{j}] \end{array}

(8)

Taking the generalized form of g(x) as G_ij (X_ij), equations (6) and (7) can be simplified to the following form:

M_{i j} \frac{d^{2} X_{i j}}{d t} + C_{i j} \frac{d X_{i j}}{d t} + K_{i j} G_{i j} (X_{i j}) = F_{i j}

(9)

The primary methods for solving the nonlinear dynamic equation (9) include the transfer matrix method, the finite element method, and the modal synthesis method. By applying appropriate boundary conditions, the eigenvalues of the equation can be obtained.

Variation patterns of modal characteristics

Case parameters

Taking a multi-shaft IGC for energy storage as an example, a variable-frequency drive motor is employed to accommodate the variable operating conditions of the air storage system. The unit is designed with a power of 20 MW and an input speed of 1800 r/min. The support bearings for the three high-speed shafts are of the five-pad tilting-pad type, while the drive shaft uses cylindrical sleeve bearings. Other design parameters are listed in Table 1.

Table 1.

Design parameters of the rotor system.

Shaft components	Tooth number	Bearing bore diameter (mm)	Bearing width(mm)	Bearing span(mm)
BGC	239	300	150	620
H1	38	150	150	610
H2	30	105	100	610
H3	21	90	90	530

Inherent characteristics under coupled and decoupled modes

Under the design operating conditions, the inherent characteristics of each independent shaft system were analyzed, excluding gear meshing coupling effects. Within a frequency range extending to 1.5 times the maximum design operating frequency of the highest-speed shaft system (approximately 500 Hz), the calculated natural frequencies of each shaft system are listed in Table 2. The results reveal that the 1st and 2nd, 3rd and 4th, as well as the 5th and 6th natural frequencies are relatively close to each other, and their corresponding mode shapes are similar (Figure 4), exhibiting either linear tilting or bending modes. This is attributed to the difference in stiffness values between the X and Y directions of the support bearings. Within the examined frequency range, H1 and H3 exhibit six natural frequencies, while BGC and H2 show four natural frequencies. As depicted in Figure 4, the displacement amplitudes of the axial nodes of rotors H1, H2, and H3 are contingent upon the distinct natural frequencies. For rotor H1, the nodes attaining the maximum amplitude within the initial four mode shapes are positioned at the two extremities of the shaft, whereas those in the subsequent two mode shapes are situated at the mid-span region. Rotor H3 demonstrates an identical evolutionary pattern. For rotor H2, the nodes exhibiting the peak amplitude in the first two mode shapes are localized at the right-hand end of the shaft (i.e., the overhung impeller terminus), while those in the latter two mode shapes are located at the shaft’s mid-span.

Table 2.

The first six natural frequencies of uncoupled rotors.

Shaft components	1^st(Hz)	2^nd(Hz)	3^rd (Hz)	4^th (Hz)	5^th (Hz)	6^th (Hz)
BGC	94.3	106.5	371.8	373.0	/	/
H1	111.7	115.9	152.6	157.2	418.0	437.4
H2	176.6	182.9	323.7	351.5	/	/
H3	211.3	235.7	260.4	285.4	408.7	479.9

Figure 4.

The mode shapes of rotor: (a) H1, (b) H2, and (c) H3.

Considering the gear meshing coupling effects, an inherent characteristic analysis of the rotor system was carried out. The resulting first 30 natural frequencies of each shaft are presented in Table 3. Due to the gear meshing effect, the composition of the system’s natural frequencies becomes more complex. While some natural frequencies are close to those of the individual uncoupled rotors, their corresponding mode shapes exhibit distinct characteristics. The natural frequencies of the coupled rotor system manifest in three primary types.

(1) Exclusive to a Single Shaft: Certain modes are dominated by a single shaft, such as Modes 6, 9, and 12 in Table 3, where the respective shaft accounts for 100% of the kinetic energy. Comparing the first four mode shapes and phases of rotor H2 before and after meshing, it can be observed from Figure 5(a) that the variation trends of the first two mode shapes after meshing are similar. Specifically, the maximum nodal amplitude exhibits a substantial reduction, whereas the phase value undergoes only a minor change. The third mode shape shows considerable variations in both amplitude and phase. For the fourth mode shape, the amplitude change is relatively small, but the phase value experiences a notable shift (Figure 5(b)).

(2) Shared by Multiple Shafts via Strong Meshing Coupling: In these modes, multiple shafts are strongly coupled through gear teeth and share the kinetic energy. Examples include Modes 3, 5, 8, and 11. Mode 11 is particularly representative, where the BGC, H1, and H2 shafts share this natural frequency with kinetic energy distributions of 30.5%, 30.6%, and 38.9%, respectively. The corresponding mode shapes (in the X and Y directions) are shown in Figure 6. Rotor H2 exhibits the greatest vibration amplitude and the most substantial phase alteration, which aligns with its predominant energy contribution. In contrast, due to their larger inertial masses and smaller kinetic energy fractions, both BGC and rotor H1 show relatively attenuated peak amplitudes and phase deviations.

(3) Weakly Coupled, Dominated by a Single Shaft: These modes, such as Modes 1, 2, 4, and 7, exhibit weak coupling effects. One shaft plays a dominant role, possessing an absolute kinetic energy advantage (≥95% of the total kinetic energy). It is noted that while Modes 1 & 2 and Modes 10 & 11 have identical calculated natural frequency values, they correspond to distinct eigenvectors, reflected in their different mode shapes and kinetic energy distributions.

Table 3.

First 30 natural frequencies and kinetic energy distribution of the rotor system with gear meshing coupling.

Mode order	Frequencies (Hz)	Belonging to	Kinetic energy distribution (%)	Mode order	Frequencies (Hz)	Belonging to	Kinetic energy distribution (%)
1	76.4	BGC	98.5	16	271	H3	100
2	76.4	BGC	97.0	17	287.8	H2	99.9
3	85.3	BGC,H1	84.5,15.0	18	330.4	H1	98.4
4	95.5	H1	98.6	19	332.4	H2	98.8
5	100.8	BGC,H1	79.4,17.7	20	380.5	H1	99.9
6	106	H1	100	21	396.9	H1	100
7	120.5	M	96.4	22	429.4	BGC	99.7
8	126.4	BGC,H1	8.6,91.2	23	429.5	BGC	99.8
9	134.5	H1	100	24	430.3	H1	99.9
10	155.2	H2	98	25	469.8	H3	100
11	155.2	BGC,H1,H2	30.5,30.6,38.9	26	474.9	H3	99.6
12	171.3	H2	100	27	548.2	H1	99.3
13	206	H3	100	28	550.8	H2	99.8
14	227.5	H3	100	29	586.9	H1	100
15	252	H3	100	30	657.4	H2	100

Figure 5.

The mode shapes and phases of rotor H2 at distinct orders prior to and subsequent to coupling: (a) the 1st and 2nd orders, and (b) the 3rd and 4th orders.

As shown in Tables 2 and 3 and Figures 5 and 6, incorporating gear-pair meshing coupling significantly increases the modal complexity of the multi-shaft geared rotor system. This coupling not only alters the natural frequencies, mode shapes, and phases of the individual shafts but also introduces numerous system-level natural frequencies exhibiting distinct strongly and weakly coupled characteristics.

Figure 6.

Mode shapes and phases of individual shafts under shared Mode 11.

Orthogonal experimental design

In CAES, the output flow and pressure of a multi-shaft IGC must be regulated to adapt to variations in the pressure and temperature of the terminal gas storage environment. During the regulation of a multi-shaft IGC, changes in input power and operating speed also alter the inherent dynamic characteristics of its rotor system. Owing to the increased number of natural frequencies resulting from gear-mesh coupling, many of which lie within 1.5 times the operating frequency of individual shafts, both the values and the number of these frequencies critically constrain the stability of the geared rotor system. To investigate this influence, the key operating parameters power, rotational speed, and inlet oil temperature were selected for a focused study on their effects on the inherent characteristics of the rotor system. Numerical analyses were performed over the following ranges: 0.5–1.1 times the rated power, 0.5–1.1 times the rated speed, and an inlet oil temperature of 36–46 °C. To maintain computational efficiency without compromising the validity of the results, an orthogonal experimental design was adopted. The design considered three factors at seven levels each (see Table 4) and employed an L₄₉ (7³) orthogonal array, requiring 49 test runs. Interactions between factors were neglected in this design.

Table 4.

Orthogonal experimental parameters.

Factor/Level	1	2	3	4	5	6	7
Power Ratio (PR)	0.5	0.6	0.7	0.8	0.9	1.0	1.1
Speed Ratio (SR)	0.5	0.6	0.7	0.8	0.9	1.0	1.1
Inlet Oil Temperature (IOT)(°C)	34	36	38	40	42	44	46

Influence of parameters on modal characteristics

Through the orthogonal experiment, natural frequency data were collected for the BGC, H1, H2, and H3 shafts under each test condition, corresponding to both tilting and bending mode characteristics. Data for shared-mode frequencies and their kinetic energy ratios were also obtained, with detailed results presented in Table 5. The table lists distinctive exclusive natural frequencies (orders 6, 9, 12, 14, 16, and 21) as well as shared natural frequencies. The results indicate that under variable operating conditions, both the exclusive natural frequencies of individual shafts and the shared natural frequencies of the multi-shaft system undergo considerable variation. Concurrently, the kinetic energy distribution exhibits significant changes, reflecting increased complexity in the inherent dynamic characteristics. In addition to two-shaft shared frequencies, cases of three-shaft shared frequencies emerged, such as the coupled natural frequencies in test runs No. 24, 36, and 49, which are shared by the H1, BGC, and H2 shafts.

Table 5.

Orthogonal experimental data.

No.	PR	SR	IOT	6^th	9^th	12^th	14^th	16^th	21^th	Shared	H1	BGC	H2
1	0.8	0.8	40	106.9	135.2	174.7	228.9	272.6	399.0	157.9	42.0	55.8	/
2	1.1	0.8	36	112.0	140.6	181.5	238.6	285.7	429.4	166.4	33.9	65.9	/
3	0.6	0.8	38	106.1	129.8	162.1	203.0	248.9	381.6	150.2	58.7	41.2	/
4	1.1	1.1	34	105.7	134.2	171.4	226.9	270.3	396.2	154.3	53.6	46.1	/
5	0.9	0.9	38	105.4	134.0	171.9	228.1	271.7	393.5	154.2	48.5	51.1	/
6	0.9	1	42	104.1	132.9	168.3	223.6	267.0	391.1	153.8	51.6	48.2	/
7	0.7	0.5	34	113.1	142.1	183.4	240.4	288.9	429.6	172.2	25.5	73.6	/
8	0.7	1.1	44	96.9	126.4	155.4	205.8	249.4	374.3	147.9	63.5	36.5	/
9	0.6	1	46	94.8	125.5	175.3	203.0	248.9	372.0	147.6	62.8	37.2	/
10	1.1	0.7	46	114.4	145.6	186.0	242.5	293.2	429.7	179.8	19.6	78.5	/
11	0.5	0.6	36	103.4	132.3	167.0	221.9	265.3	388.3	154.6	46.1	53.7	/
12	1.1	0.6	42	119.8	149.4	189.1	245.0	299.0	450.6	186.7	14.8	78.9	/
13	0.7	1	40	97.9	144.6	159.0	211.0	256.7	378.1	144.6	74.7	23.3	/
14	0.6	0.9	42	97.1	127.5	157.6	208.9	252.5	376.2	148.9	60.1	39.8	/
15	0.5	0.7	40	99.4	129.2	161.0	213.6	258.0	379.7	151.2	52.7	47.2	/
16	0.5	0.9	34	100.9	123.8	164.8	198.4	244.9	369.8	145.9	69.5	30.5	/
17	0.9	0.6	40	113.9	143.8	189.6	241.7	291.5	437.4	176.0	22.1	77.1	/
18	0.5	1	38	99.8	121.8	146.6	192.9	234.2	362.2	145.1	71.4	28.6	/
19	0.5	1.1	42	98.9	120.1	144.5	188.0	230.8	353.9	144.5	72.7	27.3	/
20	0.7	0.8	46	104.4	133.2	168.7	224.0	267.4	391.2	156.5	41.7	58.0	/
21	0.9	0.8	34	104.4	133.2	168.7	224.0	267.4	391.2	159.1	44.0	55.8	/
22	1	0.5	44	129.6	153.5	197.8	246.5	302.9	473.4	197.8	10.6	87.2	/
23	1.1	0.9	40	110.2	138.5	178.2	235.3	280.7	420.1	162.5	38.2	61.6	/
24	1	1	36	106.0	134.5	171.3	227.5	271.0	396.9	155.2	30.6	30.5	38.9
25	0.6	0.7	34	103.4	132.3	167.1	222.1	265.6	388.9	153.3	51.7	48.2	/
26	1	1.1	40	104.0	132.8	168.1	223.4	266.8	391.1	153.3	53.4	46.5	/
27	0.8	1	34	102.2	130.3	163.4	217.2	261.1	384.1	150.4	60.2	39.7	/
28	0.7	0.7	42	107.4	135.7	173.4	229.9	273.6	400.5	160.2	36.9	62.7	/
29	0.8	0.9	44	103.0	132.0	168.5	223.8	267.2	387.7	154.2	46.7	53.0	/
30	0.6	0.5	40	111.4	140.0	180.2	237.2	283.5	429.4	170.3	24.6	74.9	/
31	1	0.9	46	109.0	137.2	176.1	232.9	277.4	407.1	162.0	35.9	63.5	/
32	0.5	0.5	46	108.7	138.3	175.6	232.2	276.5	405.1	167.8	23.6	75.8	/
33	0.8	0.5	42	118.6	147.9	187.2	243.5	295.3	431.3	184.5	15.4	80.2	/
34	0.9	0.7	44	111.8	140.4	181.0	238.1	284.9	429.4	169.4	27.0	72.6	/
35	0.6	1.1	36	100.7	123.3	149.7	197.5	241.0	369.8	145.5	71.3	28.7	/
36	1.1	1	44	108.4	136.7	175.1	231.9	276.0	404.9	159.6	12.9	17.6	69.5
37	0.8	0.6	46	112.8	143.0	182.8	239.7	287.7	429.5	175.5	20.5	78.7	/
38	0.5	0.8	44	95.7	126.6	155.6	205.9	247.1	373.6	148.8	58.1	41.8	/
39	1	0.7	38	112.9	141.8	183.1	240.1	288.3	429.5	171.2	26.9	71.4	2.0
40	0.9	1.1	46	101.9	131.9	164.9	219.1	262.8	385.7	148.3	69.5	28.9	/
41	0.7	0.9	36	100.5	130.1	162.9	216.4	260.4	382.9	150.6	58.2	41.8	/
42	1.1	0.5	38	128.5	154.0	197.3	247.2	304.8	481.6	197.3	11.1	86.0	/
43	1	0.8	42	111.0	139.4	179.5	236.6	282.6	429.0	165.7	31.1	64.5	/
44	0.8	1.1	38	101.5	128.6	160.0	212.5	257.0	379.7	149.1	62.6	37.4	/
45	0.9	0.5	36	118.2	148.8	188.9	244.9	298.6	447.9	185.1	16.4	80.6	/
46	0.6	0.6	44	108.0	136.3	174.4	230.8	274.8	402.5	163.0	31.0	68.6	/
47	0.7	0.6	38	110.3	138.6	178.4	235.4	280.8	421.3	155.5	43.1	55.9	/
48	1	0.6	34	114.7	144.5	186.6	243.2	294.6	430.0	175.4	24.1	73.1	/
49	0.8	0.7	36	109.5	137.7	176.9	233.9	278.6	411.6	162.0	32.6	54.4	13.1

Influence of power on modal characteristics

The influence of power variation on the natural frequencies of the rotor system is shown in Figure 7. The results indicate that the natural frequencies of each order increase with rising power, exhibiting a positive correlation; the derived coupled (shared) natural frequencies follow the same trend. Regarding the kinetic energy ratio under shared frequencies, that of H1 shows a decreasing trend, while BGC exhibits an increasing trend. All characteristic variations demonstrate weakly nonlinear behavior.

Figure 7.

Influence of power ratio variation on modal characteristics.

Influence of rotational speed on modal characteristics

The influence of rotational speed variation on the natural frequencies is presented in Figure 8. It can be seen that the natural frequencies of all orders decrease with increasing speed ratio, showing a negative correlation, a trend also followed by the derived coupled (shared) natural frequencies. Regarding the kinetic energy ratio under shared frequencies, that of H1 increases with the speed ratio, whereas BGC first decreases and then increases. All characteristic variations exhibit weakly nonlinear behavior. At a speed ratio of 1, three-shaft shared natural frequencies appear in test runs 24 and 36. The corresponding kinetic energy distributions are 30.6 % (H1), 30.5 % (BGC), and 38.9 % (H2) for the first case, and 12.9 % (H1), 17.6 % (BGC), and 69.5 % (H2) for the second. In these instances, the kinetic energy share of H2 increases while those of H1 and BGC decrease, a pattern consistent with the trend observed in Figure 8 at the speed ratio of 1.

Figure 8.

Influence of speed ratio variation on modal characteristics.

Influence of inlet oil temperature on modal characteristics

The influence of inlet oil temperature variation on the natural frequencies is shown in Figure 9. The data demonstrate that as the oil temperature changes, some natural frequencies (e.g., orders 6, 9, 12, and the coupled frequencies) follow a pattern of increasing, then decreasing, then increasing again, and finally decreasing. In contrast, other frequencies (e.g., orders 14 and 16) exhibit the opposite trend. For the kinetic energy ratios of the shafts under shared natural frequencies, both H1 and BGC show a pattern of decreasing, increasing, decreasing, and then increasing again. All observed trends demonstrate considerable variability.

Figure 9.

Influence of inlet oil temperature variation on modal characteristics.

Sensitivity analysis

To determine the degree of influence of power, rotational speed, and inlet oil temperature on the inherent characteristics of the rotor system, range analysis was performed on the results of the orthogonal experiment. Table 6 presents the range data, where K_i represents the average value of the i-th level of a given factor and R denotes the range value. The data reveal that for indicators such as the 6^th order frequency, the shared frequencies, and the kinetic energy ratio of the H1 shaft, the range values follow the order SR > PR > IOT. The same trend is observed for other frequency indicators. Based on these results, it can be concluded that the parameters, ranked in descending order of their influence on the natural frequencies, are: SR > PR > IOT.

Table 6.

Range analysis.

Factors	Indicators	K₁	K₂	K₃	K₄	K₅	K₆	K₇	R
PR	6^th	101.0	103.1	104.4	107.8	108.5	112.5	114.1	13.1
	Shared	151.1	154.1	155.4	161.9	163.7	168.7	172.4	21.3
	H1	56.3	51.5	49.1	40.0	39.9	30.4	26.3	30.0
SR	6^th	118.3	111.8	108.4	105.8	103.7	101.9	101.4	16.9
	Shared	182.1	169.5	163.9	157.8	154.0	150.9	149.0	33.1
	H1	18.2	28.8	35.3	44.2	51.0	52.0	63.8	45.6
IOT	6^th	106.3	107.2	109.2	106.2	108.1	107.6	106.6	2.9
	Shared	158.7	159.9	160.4	159.4	163.5	163.0	162.5	4.8
	H1	46.9	41.3	46.0	44.0	40.4	35.7	39.1	11.2

Variance analysis

To examine whether significant correlations exist between the various factors and the natural frequencies, a variance analysis was conducted. The resulting analysis of variance (ANOVA) table is presented in Table 7. In the table, the degrees of freedom for each factor are 6, and the error term is 30, where an asterisk (*) indicates a significant correlation. Taking α = 0.05, the critical value is F₀.₀₅ (6,30) = 2.42. All calculated F-values for the PR and the SR exceed this threshold, indicating that both factors show significant correlation with the natural frequencies. The IOT exhibits a significant effect on the kinetic energy ratio of H1, but no significant effect on the 6^th order frequency or the shared frequency itself. Based on the magnitude of the F-values, the factors can be ranked in descending order of significance as: SR> PR > IOT, which is consistent with the results of the range analysis.

Table 7.

Analysis of variance table.

Factors	Indicators	Sum of squares	Degrees of freedom	F-value	p-value	Significance
PR	6^th	992.261	6	25.625	0	*
	Shared	2609.668	6	28.737	0	*
	H1	5140.278	6	20.271	0	*
SR	6^th	1556.461	6	40.195	0	*
	Shared	5830.905	6	64.209	0	*
	H1	10135.301	6	39.97	0	*
IOT	6^th	49.21	6	1.271	0.30
	Shared	152.605	6	1.68	0.16
	H1	672.41	6	2.652	0.035	*
Errors	6^th	193.612	30
	Shared	454.059	30
	H1	1267.863	30

Multiple linear regression prediction model

Drawing from the parameter variation patterns illustrated in Figures 7 and 9 and the orthogonal experimental data provided in Table 5, a multiple linear regression model was developed to predict the natural frequencies of different orders, the coupled (shared) frequency, and the kinetic energy ratios within the studied parameter ranges.

Define y_i as the i-th natural frequency of the rotor system, where i ≤ 30. In addition, let y₃₁, y₃₂, and y₃₃ denote the shared modal frequency, the kinetic energy ratio of the H1 and the BGC shaft, respectively. The independent variables are defined as x₁ for PR, x₂ for SR, and x₃ for IOT.

For the i-th dependent variable, the following relationship is formulated:

y_{i} = a_{0 i} + a_{1 i} x_{1 i} + a_{2 i} x_{2 i} + a_{3 i} x_{3 i} + e_{i} i \leq 33

(10)

Let:

\begin{array}{l} Y = {(y_{1}, y_{2}, \dots, y_{33})}^{T} \\ A = {(a_{01}, a_{11}, a_{21}, a_{31}; a_{02}, a_{12}, a_{22}, a_{32}; \dots, a_{0, 33}, a_{1, 33}, a_{2, 33}, a_{3, 33})}^{T} \\ E = {(e_{1}, e_{2}, \dots, e_{33})}^{T} \\ X = {(1, x_{1 i}, x_{2 i}, x_{3 i})}^{T} \end{array}

(11)

Hence, the multiple linear regression model takes the following form:

Y_{33 \times 1} = A_{33 \times 4} X_{4 \times 1} + E_{33 \times 1}

(12)

On the basis of the orthogonal experimental data, regression analysis was performed on Equation (12) to obtain the regression coefficients. Partial results are listed in Table 8. A comparison between the orthogonal experimental results and the fitted results for the 6^th-order frequency and the shared frequency is presented in Figure 10. The small discrepancies observed verify the accuracy of the regression prediction model. In practical engineering applications, the established prediction model can be directly used to guide the dynamic control strategies of multi-shaft IGCs under variable operating conditions.

Table 8.

Regression coefficients.

y _i /a _i	a _0i	a _1i	a _2i	a _3i	R²
6^th	110.67	22.31	-26.92	0.01	0.86
9^th	133.09	24.13	-27.79	0.14	0.87
12^th	168.06	37.49	-43.65	0.23	0.91
14^th	218.59	50.79	-52.74	0.22	0.93
16^th	267.28	59.30	-64.40	0.21	0.96
21^th	404.99	87.31	-102.83	0.29	0.92
Shared	159.21	36.13	-52.34	0.37	0.89
H1	54.36	-50.50	71.07	-0.72	0.89
BGC	67.06	33.05	-78.38	0.60	0.86

Figure 10.

Comparison between fitted and experimental results.

Engineering case validation

To acquire vibration data from the rotor systems of the multi-shaft IGC, operational data were collected for the rotor system with the parameters listed in Table 1. The schematic diagram of the established test procedure is shown in Figure 11. The vibration acceleration measurement points are located at the H1 bearing housing in the X and Y directions (orthogonal, i.e., 90° apart), and the sampling frequency range is set to at least four times the maximum operating rotational speed. The shaft vibration measurement points are located near the bearing positions of H1 and H2 on the motor side, and the sensor type is a non-contact eddy current sensor.

Figure 11.

Schematic diagram of data acquisition.

The captured frequency-domain and time-domain vibration characteristics allow for targeted analysis. In this case study, the operating conditions during data acquisition were: SR= 1, PR = 1, and IOT = 36 °C. The frequency-domain and time-domain vibration data obtained under these conditions for the X and Y directions are shown in Figures 12 and 13, respectively.

Figure 12.

Frequency-domain plots.

Figure 13.

Time -domain plots.

Figure 12 reveals that the frequency-domain data contain peaks corresponding to the rotational speeds of the individual rotors, as well as frequency components arising from system-level natural frequencies induced by gear meshing. The measured value (158 Hz) shows an error of 1.8 % relative to the computed value (155.2 Hz) and an error of 1.1 % relative to the predicted value (156.3 Hz), which validates the accuracy of both the computational model and the regression prediction model. Although the vibration amplitude associated with the coupled frequency is relatively low, it still contributes a noticeable share to the overall response. This demonstrates that in practical engineering applications, the coupling effect introduced by gear meshing cannot be neglected.

Under operating conditions with a PR of 0.7 and an IOT of 44 °C, the multi-shaft IGC was regulated by reducing the SR from 1.0 to 0.7. During this process, the monitoring system captured the shaft vibration data for H1 and H2. The curves of vibration amplitude versus speed ratio are shown in Figure 14. It can be observed that the critical speed ratios for both H1 and H2 are 0.72, corresponding to rotational frequencies of 135.9 Hz and 172.1 Hz, respectively. Each resonant frequency corresponds to the 9^th and 12^th natural frequencies of the rotor system. As derived from the multiple linear regression model equation, the predicted frequencies for this operating condition are 136.1 Hz and 173.0 Hz, with errors within 1 Hz. This result further validates the accuracy of the prediction model. Based on this model, the inherent characteristics of the IGC rotor system under CAES variable operating conditions can be predicted, thereby providing a theoretical foundation for system control strategies.

Figure 14.

Vibration curves of the H1 and H2.

Conclusion

Variable operating conditions pose a significant challenge to the stability of multi-shaft IGC rotor systems used in CAES. This study investigates the evolution of natural frequencies and mode shapes of the rotor system under both single and variable operating conditions, accounting for multi-gear meshing coupling effects. The sensitivity of inherent characteristics to key control parameters is analyzed, and a multiple linear regression model is developed and validated against engineering data. The main findings are summarized as follows.

(1) Before and after multi-gear-pair meshing coupling, the natural frequencies and mode shapes of each rotor undergo noticeable changes, while several new coupled shared natural frequencies are generated, making the modal behavior of the rotor system more complex.

(2) Under variable operating conditions, rotational speed and power have a significant influence on the modal characteristics of the rotor system, mainly reflected in variations in frequency values and kinetic energy distribution, following a weakly nonlinear trend. The parameters are ranked in descending order of sensitivity as: rotational speed, power, and inlet oil temperature.

(3) Vibration data collected from the engineering case validate the accuracy of the developed multiple linear regression prediction model, which provides important guidance for the design of multi-shaft IGCs used in CAES under variable operating conditions.

Footnotes

ORCID iD

Jizhe Mao

Author contributions

Jizhe Mao was responsible for data curation and drafting the initial manuscript. Yongchao Yuan conducted the formal analysis and applied software tools. Jianning Gong provided critical review and editing of the manuscript. All authors approved the final version.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Key R&D Program Project of Shandong Province, China (Grant Nos. 2023CXPT108 and 2024CXPT038).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data Availability Statement

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.*

Appendix

References

Qiu

Yuan

Han

, et al. Performance investigation of a novel power-cooling-heating-compressed air multigeneration system based on compressed air-thermal combined energy storage. Energy Conversion And Management 2026; 351: 121027. https://doi.org/10.1016/j.enconman.2025.121027

Tang

Huang

, et al. Comparative study of analytical and numerical solutions for the thermodynamic processes of compressed air energy storage reservoirs in abandoned mines. Renewable Energy 2026; 260: 125245. https://doi.org/10.1016/j.renene.2026.125245

Sun

Zhang

, et al. Thermal-flow coupling mechanism of multi-nozzle air injection for compressed air energy storage in abandoned mine roadway. Energy 2026; 342: 139586. https://doi.org/10.1016/j.energy.2025.139586

Fang

Hou

. The impact of energy supply side on the diffusion of low carbon transformation on energy demand side under low-carbon policies in China. Energy 2024; 307: 132817. https://doi.org/10.1016/j.energy.2024.132817

Zhang

Gao

Zhou

, et al. Advanced Compressed Air Energy Storage Systems: Fundamentals and Applications. Engineering 2024; 34: 246–269. https://doi.org/10.1016/j.eng.2023.12.008

Barbour

Pottie

DLF

. In: Cabeza

(ed). Adiabatic Compressed Air Energy Storage Systems. Encyclopedia of Energy Storage. Elsevier, 2022, pp. 188–203.

Braff

Mueller

Trancik

. Value of storage technologies for wind and solar energy. Nature Climate Change 2016; 6(10): 964–969. https://doi.org/10.1038/nclimate3045

Budt

Wolf

Span

, et al. A review on compressed air energy storage: Basic principles, past milestones and recent developments. Applied Energy 2016; 170: 250–268. https://doi.org/10.1016/j.apenergy.2016.02.108

Bazdar

Sameti

Nasiri

, et al. Compressed air energy storage in integrated energy systems: A review. Renewable & Sustainable Energy Reviews 2022; 167: 112701. https://doi.org/10.1016/j.rser.2022.112701

10.

Agyekum

Abdullah

Rashid

, et al. Compressed air energy storage (CAES) systems: technological progress, challenges, and future prospects in renewable energy grids. Energy Conversion and Management: X 2025; 28: 101402. https://doi.org/10.1016/j.ecmx.2025.101402

11.

Miao

Luo

, et al. Dynamic modelling and technoeconomic analysis of adiabatic compressed air energy storage for emergency back-up power in supporting microgrid. Appl Energy 2020; 261: 114448. https://doi.org/10.1016/j.apenergy.2019.114448

12.

Guo

Zhang

, et al. Dynamic characteristics and control of supercritical compressed air energy storage systems. Appl Energy 2021; 283: 116294. https://doi.org/10.1016/j.apenergy.2020.116294

13.

Zhang

Chen

, et al. Distributed generation with energy storage systems: a case study. Appl Energy 2017; 204: 1251–1263. https://doi.org/10.1016/j.apenergy.2017.05.063

14.

Liang

Zuo

Zhou

, et al. Design of a Centrifugal Compressor with Low Solidity Vaned Diffuser (LSVD) for Large-Scale Compressed Air Energy Storage (CAES). J Therm Sci 2020; 29(2): 423–434. https://doi.org/10.1007/s11630-019-1204-7

15.

Wang

, et al. Dynamic characteristics of the gear-rotor system in compressed air energy storage considering friction effects. Mech Sci 2021; 12(1): 677–688. https://doi.org/10.5194/ms-12-677-2021

16.

Yang

Sun

, et al. Dynamic characteristics and operation strategy of the discharge process in compressed air energy storage systems for applications in power systems. Int J Energy Res 2020; 44(8): 6363–6382. https://doi.org/10.1002/er.5362

17.

Song

Zhou

Zhang

, et al. A comprehensive method for controlling the electrical runout of the rotor for a geared centrifugal compressor. Review of Scientific Instruments 2021; 92(6): 065104. https://doi.org/10.1063/5.0048285

18.

Pelton

Jung

Allison

, et al. Design of a Wide-Range Centrifugal Compressor Stage for Supercritical CO2 Power Cycles. Journal of Engineering for Gas Turbines and Power 2018; 140(9): 092602. https://doi.org/10.1115/1.4039835

19.

Modekurti

Eslick

Omell

, et al. Design, dynamic modeling, and control of a multistage CO2 compression system. International Journal of Greenhouse Gas Control 2017; 62: 31–45. https://doi.org/10.1016/j.ijggc.2017.03.009

20.

Zhang

Jiang

Gao

. Dynamic Analysis of Integrally Geared Compressors with Varying Workloads. Shock and Vibration 2016; 2016: 1–13. https://doi.org/10.1155/2016/2594635

21.

Yang

Wang

, et al. A general method to predict unbalance responses of geared rotor systems. Journal of Sound and Vibration 2016; 381: 246–263. https://doi.org/10.1016/j.jsv.2016.06.031

22.

Zhang

Zhai

Han

, et al. Dynamics of a geared parallel-rotor system subjected to changing oil-bearing stiffness due to external loads. Finite Elements in Analysis and Design 2015; 106: 32–40. https://doi.org/10.1016/j.finel.2015.07.008

23.

Wang

Zhu

, et al. Vibration characteristics of triple-gear-rotor system in compressed air energy storage under variable torque load. Science Progress 2021; 104(1): 0036850420987058. https://doi.org/10.1177/0036850420987058

24.

Yang

Cai

Zhuge

, et al. Influence of Real Gas Properties on Aerodynamic Losses in a Supercritical CO2 Centrifugal Compressor. J Therm Sci 2024; 33(6): 2032–2046. https://doi.org/10.1007/s11630-024-2027-8

25.

Lin

Zuo

Sun

, et al. Flow characteristics of impeller backside cavity and its effects on the centrifugal compressor for compressed air energy storage. Journal of Energy Storage 2022; 49: 104024. https://doi.org/10.1016/j.est.2022.104024

26.

Wang

, et al. Investigation of methods to enhance the stable operating range of high-pressure ratio centrifugal compressors adopted in the compressed air energy storage system. Journal of Energy Storage 2025; 121: 116539.

27.

Liu

Pang

, et al. Dynamic modeling and vibration analysis of a flexible gear transmission system. Mechanical Systems and Signal Processing 2023; 197: 110367. https://doi.org/10.1016/j.ymssp.2023.110367

28.

Tang

, et al. Three-dimensional vibration investigation of the thin web gear pair based on Timoshenko beam. Thin-Walled Structures 2023; 184: 110507. https://doi.org/10.1016/j.tws.2022.110507

29.

Babu

Tandon

Pandey

. Nonlinear Vibration Analysis of an Elastic Rotor Supported on Angular Contact Ball Bearings Considering Six Degrees of Freedom and Waviness on Balls and Races. Journal of Vibration and Acoustics 2014; 136(4): 044503, 044503. 29. https://doi.org/10.1115/1.4027712