1 INTRODUCTION
Magnetic bearings offer frictionless high-speed operation, high efficiency, low energy consumption, and active controllability, making them a promising support technology for aero-engines and steam turbines [1]. To further investigate the coupled dynamic behavior of magnetic bearing–rotor systems, researchers have developed various high-precision dynamic models based on their electromagnetic properties, control logic, and dynamic response. These modeling approaches have been integrated with rotor-dynamics theory, enabling systematic investigation of the dynamic behavior of magnetic-bearing–rotor systems [2,3]. Gao et al. [4] investigated the fundamental dynamic characteristics of magnetically suspended rotors by formulating the AMB dynamics on the basis of rotor-dynamics and vibration theory. Wu et al. [5] proposed an analytical dynamic modeling approach for an active magnetic bearing (AMB)–supported rotor system in a magnetic levitation steam turbine. The flexible rotor is discretized using Euler–Bernoulli beam theory, while the disk–shaft connection is represented by an equivalent spring–damper model with stiffness and damping parameters identified via three-dimensional finite element analysis; the AMB control forces are also incorporated into the model. Narsakka et al. [6] extended conventional analytical methodologies to assess rotor stress and fatigue life under drop-down fault conditions in active magnetic bearing systems. The rotor finite element model, incorporating a nonlinear touchdown bearing contact model, was developed to simulate the dynamic process of rotor descent following levitation loss. Despite superior high-speed performance, low friction, and controllability relative to mechanical bearings, magnetic bearings face significant thermal challenges in operation. Uneven heat dissipation among windings causes localized heating and thermal imbalance, while temperature gradients during high-speed rotation can induce rotor deformation, known as thermal bowing, degrading system stability and precision [7]. Regarding the thermal and thermo-mechanical behavior of magnetic-bearing–rotor systems, extensive theoretical and experimental studies have been conducted in recent years. Badykov et al. [8] proposed an improved hybrid axial electromagnetic bearing design that effectively reduces power consumption and heat generation, and performed numerical simulations using two-dimensional magnetostatic and two- and three-dimensional thermal finite element analyses. Takahashi et al. [9] attributed the synchronous vibration instability near critical speeds in a magnetic bearing‑supported, two‑stage overhung centrifugal compressor to the coupling of iron‑loss concentration and thermal bending, clarifying its generation mechanism through theoretical modeling. Xu et al. [10] proposed an integrated multi‑physics simulation framework that couples electromagnetic, thermal, and mechanical fields, enabling quantitative analysis of cross‑scale interactions. This framework provides a systematic approach for investigating thermally induced vibration in magnetically suspended rotors. Jin et al. [11] integrated ANSYS Maxwell with COMSOL to develop a magneto-thermal coupling approach, which was used to analyze the mechanism and parametric influences of iron‑loss concentration and the resulting temperature distribution in magnetic bearings, providing a theoretical foundation for addressing related loss and thermal‑imbalance issues. Zhao et al. [12] developed a simplified thermoelastic rotordynamic model by deriving the thermal potential energy under arbitrary temperature fields and establishing the governing equations of motion via the Euler–Lagrange method. Cui et al. [13] established a dynamic model for an axially loaded thin-walled shaft–disk rotor under non-uniform temperature fields, taking into account temperature-dependent material properties, gyroscopic effects, and centrifugal forces. To further investigate nonlinear dynamics and bidirectional thermo-mechanical coupling in rotor systems, Chang et al. [14] derived the equations of motion using the Lagrangian method while accounting for thermally induced variations in bearing radial clearance. A thermal network model, combining lumped-parameter and full-network approaches, was established for the bearing and coupled with the rotor dynamic equations through the bearing mechanical model. In addition, the same research group proposed a thermo-mechanically coupled model for a dual-rotor system, incorporating multiple nonlinear factors as well as heat generation effects in four bearings, thereby enabling comprehensive analysis of nonlinear characteristics under thermo-mechanical coupling conditions [15]. Yu et al. [16] constructed a spatially distributed heating setup to simulate non‑uniform thermal environments and conducted experiments on the thermally induced bending deformation and dynamic behavior of an aero‑engine high‑pressure rotor. The tests systematically measured rotor temperature profiles, structural deformation, and vibration response. Goldman et al. [17] established a quasi-static rotor heat transfer model based on contact heat generation to derive the complex-variable differential equation for thermal bending, investigating its effects along with imbalance and radial loads on intermittent rub-impact vibration. Larsson [18,19] modeled rotor vibration and journal temperature to analyze unsteady bearing heat transfer, developing a dynamic method for predicting shaft thermal bending vibrations.
Although extensive studies have been conducted on the dynamic modeling of active magnetic bearing (AMB) systems, temperature field analysis, and thermo-mechanical coupling of rotor systems, the dynamic characteristics of rotor systems under the combined effects of AMB heat losses and ambient temperature have not yet been systematically investigated. To address this gap, a thermally coupled rotordynamic model of an AMB-supported rotor system is developed, in which thermally induced loads derived from the temperature field, disk unbalance, gravity, and the equivalent stiffness and damping of the AMB governed by PID control parameters are incorporated. The proposed model is validated by comparison with finite element modal results and with published trends in the equivalent stiffness and damping of AMBs. Based on the developed model, the effects of bias current, rotational speed, and ambient temperature on the system’s dynamic response are systematically analyzed. This study provides a theoretical basis for thermal stability assessment and vibration control of AMB-supported rotor systems under high-temperature conditions.
2 METHODS AND MATERIALS
Given the complexity of dynamic modeling for the multiphysics-coupled system, the following assumptions and simplifications are adopted in this study:
- The AMB rotor is rigidly connected to the shaft–disk system without relative rotation, and is considered an integral part of the shafting system;
- The rotor operates under thermal steady-state conditions, and the temperature field is therefore assumed to be time-invariant;
- Only the thermal bending moments about the X and Y axes are considered, and they are treated as quasi-steady excitations in the vibration analysis.
- The ambient temperature is imposed as an external thermal boundary condition for the AMB-supported rotor system. Thermally induced rotor bending is considered, whereas the temperature-dependent variations in the electromagnetic properties of the AMB materials and the control parameters are neglected.
2.1 Modeling of Shaft and Disk System Dynamics
2.1.1 Shaft Segment
A Timoshenko beam element, shown schematically in Fig. 1, is employed to model individual shaft segments. Each node of the Timoshenko beam element has five degrees of freedom: three translational displacements along the X, Y, and Z axes and two rotational displacements about the X and Y axes. The element therefore has ten generalized coordinates, expressed as follows:[/vc_column_text]
Fig. 1. Timoshenko beam element model
Using Timoshenko beam theory with five degrees of freedom per node, the mass, gyroscopic, and stiffness matrices of the shaft element are denoted by \( \mathbf{M}_{all} \), \( \mathbf{G}_{all} \), and \( \mathbf{K}_{all} \), respectively [20]. Their associated parameters and the formulations of these matrices are provided in the following equations:
where \(\Phi\) denotes the shear factor, \(E\) represents the elastic modulus of the shaft material, \(I\) indicates the cross-sectional second moment of area of the shaft element, \(L\) denotes the length of the shaft element, \(G\) represents the shear modulus of the material, \(A\) represents the cross-sectional area of the shaft element. The shear correction coefficient \(\mu\) depends on the shaft geometry. For the solid circular shaft considered in this study, it is calculated as \(\frac{7 + 6v}{6(1 + v)}\), and \(v\) denotes Poisson’s ratio of the material.
\tag{3}$$
\mathbf{M}_{R}
=
\frac{\rho A L}{(1+\Phi)^2}
\begin{bmatrix}
M_{r1} & 0 & 0 & 0 & M_{r4} & M_{r5} & 0 & 0 & 0 & -M_{r6} \\
0 & M_{r1} & 0 & -M_{r4} & 0 & 0 & M_{r5} & 0 & M_{r6} & 0 \\
0 & 0 & M_{r3} & 0 & 0 & 0 & 0 & M_{r8} & 0 & 0 \\
0 & -M_{r4} & 0 & M_{r2} & 0 & 0 & -M_{r6} & 0 & M_{r7} & 0 \\
M_{r4} & 0 & 0 & 0 & M_{r2} & M_{r6} & 0 & 0 & 0 & M_{r7} \\
M_{r5} & 0 & 0 & 0 & M_{r6} & M_{r1} & 0 & 0 & 0 & -M_{r4} \\
0 & M_{r5} & 0 & -M_{r6} & 0 & 0 & M_{r1} & 0 & M_{r4} & 0 \\
0 & 0 & M_{r8} & 0 & 0 & 0 & 0 & M_{r3} & 0 & 0 \\
0 & M_{r6} & 0 & M_{r7} & 0 & 0 & M_{r4} & 0 & M_{r2} & 0 \\
-M_{r6} & 0 & 0 & 0 & M_{r7} & -M_{r4} & 0 & 0 & 0 & M_{r2}
\end{bmatrix}
\\[5pt]
\begin{aligned}
M_{r1} &= \frac{13}{35}+\frac{7}{10}\Phi+\frac{1}{3}\Phi^2,
&
M_{r2} &= \left(
\frac{1}{105}+\frac{1}{60}\Phi+\frac{1}{120}\Phi^2
\right)L^2,
\\
M_{r3} &= \frac{1}{3}(1+\Phi)^2,
&
M_{r4} &= \left(
\frac{11}{210}+\frac{11}{120}\Phi+\frac{1}{24}\Phi^2
\right)L,
\\
M_{r5} &= \frac{9}{70}+\frac{3}{10}\Phi+\frac{1}{6}\Phi^2,
&
M_{r6} &= \left(
\frac{13}{420}+\frac{3}{40}\Phi+\frac{1}{24}\Phi^2
\right)L,
\\
M_{r7} &= -\left(
\frac{1}{140}+\frac{1}{60}\Phi+\frac{1}{120}\Phi^2
\right)L^2,
&
M_{r8} &= \frac{1}{6}(1+\Phi)^2.
\end{aligned}
\end{gathered}
\tag{4}$$
\mathbf{M}_{T}
=
\frac{\rho I}{L(1+\Phi)^2}
\begin{bmatrix}
M_{t1} & 0 & 0 & 0 & M_{t4} & -M_{t1} & 0 & 0 & 0 & M_{t4} \\
0 & M_{t1} & 0 & -M_{t4} & 0 & 0 & -M_{t1} & 0 & -M_{t4} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & -M_{t4} & 0 & M_{t2} & 0 & 0 & M_{t4} & 0 & M_{t3} & 0 \\
M_{t4} & 0 & 0 & 0 & M_{t2} & -M_{t4} & 0 & 0 & 0 & M_{t3} \\
-M_{t1} & 0 & 0 & 0 & -M_{t4} & M_{t1} & 0 & 0 & 0 & -M_{t4} \\
0 & -M_{t1} & 0 & M_{t4} & 0 & 0 & M_{t1} & 0 & M_{t4} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & -M_{t4} & 0 & M_{t3} & 0 & 0 & M_{t4} & 0 & M_{t2} & 0 \\
M_{t4} & 0 & 0 & 0 & M_{t3} & -M_{t4} & 0 & 0 & 0 & M_{t2}
\end{bmatrix}
\\[5pt]
\begin{aligned}
M_{t1} &= \frac{6}{5},
&
M_{t2} &=
\left(
\frac{2}{15}
+\frac{1}{6}\Phi
+\frac{1}{3}\Phi^2
\right)L^2,
\\
M_{t3} &=
-\left(
\frac{1}{30}
+\frac{1}{6}\Phi
-\frac{1}{6}\Phi^2
\right)L^2,
&
M_{t4} &=
\left(
\frac{1}{10}
-\frac{1}{2}\Phi
\right)L.
\end{aligned}
\end{gathered}
\tag{5}$$
\mathbf{G}_{all}
=
\frac{\rho I}{15L(1+\Phi)^2}
\begin{bmatrix}
0 & -G_{1} & 0 & G_{2} & 0 & 0 & G_{1} & 0 & G_{2} & 0 \\
G_{1} & 0 & 0 & 0 & G_{2} & -G_{1} & 0 & 0 & 0 & G_{2} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-G_{2} & 0 & 0 & 0 & -G_{4} & G_{2} & 0 & 0 & 0 & G_{3} \\
0 & -G_{2} & 0 & G_{4} & 0 & 0 & G_{2} & 0 & -G_{3} & 0 \\
0 & G_{1} & 0 & -G_{2} & 0 & 0 & -G_{1} & 0 & -G_{2} & 0 \\
-G_{1} & 0 & 0 & 0 & -G_{2} & G_{1} & 0 & 0 & 0 & -G_{2} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-G_{2} & 0 & 0 & 0 & G_{3} & G_{2} & 0 & 0 & 0 & -G_{4} \\
0 & -G_{2} & 0 & -G_{3} & 0 & 0 & G_{2} & 0 & G_{4} & 0
\end{bmatrix}
\\[5pt]
\begin{aligned}
G_{1} &= 36,
&
G_{2} &= (3-15\Phi)L,
\\
G_{3} &= (1+5\Phi-5\Phi^2)L^2,
&
G_{4} &= (4+5\Phi+10\Phi^2)L^2.
\end{aligned}
\end{gathered}
\tag{6}$$
\mathbf{K}_{all}
=
\frac{EI}{(1+\Phi)L^3}
\begin{bmatrix}
K_{1} & 0 & 0 & 0 & K_{4} & -K_{1} & 0 & 0 & 0 & K_{4} \\
0 & K_{1} & 0 & -K_{4} & 0 & 0 & -K_{1} & 0 & -K_{4} & 0 \\
0 & 0 & K_{3} & 0 & 0 & 0 & 0 & -K_{3} & 0 & 0 \\
0 & -K_{4} & 0 & K_{2} & 0 & 0 & K_{4} & 0 & K_{5} & 0 \\
K_{4} & 0 & 0 & 0 & K_{2} & -K_{4} & 0 & 0 & 0 & K_{5} \\
-K_{1} & 0 & 0 & 0 & -K_{4} & K_{1} & 0 & 0 & 0 & -K_{4} \\
0 & -K_{1} & 0 & K_{4} & 0 & 0 & K_{1} & 0 & K_{4} & 0 \\
0 & 0 & -K_{3} & 0 & 0 & 0 & 0 & K_{3} & 0 & 0 \\
0 & -K_{4} & 0 & K_{5} & 0 & 0 & K_{4} & 0 & K_{2} & 0 \\
K_{4} & 0 & 0 & 0 & K_{5} & -K_{4} & 0 & 0 & 0 & K_{2}
\end{bmatrix}
\\[5pt]
\begin{aligned}
K_{1} &= 12,
&
K_{2} &= (4+\Phi)L^2,
\\
K_{3} &= \frac{A}{I}(1+\Phi)L^2,
&
K_{4} &= 6L,
\\
K_{5} &= (2-\Phi)L^2.
\end{aligned}
\end{gathered}
\tag{7}$$
In Equations (5-7), \(\rho\) denotes the material density, and the remaining parameters are defined in Eq. (2).
Furthermore, the model incorporates the influence of the shafting’s own weight, which is expressed in the dynamic formulation as:
=
G_w L g
\begin{bmatrix}
0 & -\frac{1}{2} & \frac{L}{12} & 0 & 0 & 0 &
-\frac{1}{2} & -\frac{L}{12} & 0 & 0
\end{bmatrix}^{\mathrm{T}},
\tag{8}$$
where \(G_{w}\) denotes the linear mass density (mass per unit length), \(g\) is the gravitational acceleration, and \(L\) has the same meaning as in Eq. (2).
2.1.2 Disk Element
The disk is modeled as a rigid finite element whose axis of rotation passes through its center of mass. The mass matrix \(\mathbf{M_{\text{disk}}}\) and the gyroscopic matrix \(\mathbf{G_{\text{disk}}}\) generated by the disk element can be represented by the following matrices, where \(m_d\) is the mass of the disk element, \(J_d\) and \(J_p\) denote the diametral and polar mass moments of inertia, respectively.
\mathbf{M}_{\mathrm{disk}}
=
\begin{bmatrix}
m_d & 0 & 0 & 0 & 0 \\
0 & m_d & 0 & 0 & 0 \\
0 & 0 & m_d & 0 & 0 \\
0 & 0 & 0 & J_d & 0 \\
0 & 0 & 0 & 0 & J_d
\end{bmatrix},
\\[10pt]
\mathbf{G}_{\mathrm{disk}}
=
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & J_p \\
0 & 0 & 0 & -J_p & 0
\end{bmatrix}.
\end{gathered}
\tag{9}$$
2.2 AMB Equivalent Stiffness and Damping Model
In this study, an 8-pole active magnetic bearing is selected as the support for the rotor system, and its schematic structure is illustrated in Fig. 2, in which each magnetic pole and corresponding region is clearly labeled.
Fig. 2. Schematic diagram of the active magnetic bearing structure
In this work, the active magnetic bearing (AMB) is modeled in the rotordynamic formulation as equivalent stiffness and damping elements obtained from linearization of the closed-loop system with fixed PID gains. Electromagnetic nonlinearities and temperature-dependent parameter variations are disregarded within the considered operating range. As a result, the equivalent stiffness and damping are not constants but quantities that vary with rotational frequency.
The relevant formulas and dynamic matrices required for the calculations are provided below [21].
The expressions for the displacement stiffness \((k_x)\) and the current stiffness \((k_i)\) are given by:
k_x = -\dfrac{\mu_0 A_b N^2 I_b^2 \cos(\alpha)}{\delta_0^3} \\
k_i = \dfrac{\mu_0 A_b N^2 I_b \cos(\alpha)}{\delta_0^2}
\end{cases},
\tag{10}$$
where \(\mu_0\) is the permeability of free space, \(A_b\) represents the area of the magnetic poles, \(N\) is the number of coil turns, \(I_b\) denotes the bias current, \(\delta_0\) represents the nominal air gap length, and \(\alpha\) corresponds to half of the inter-pole angle.
The formulas for calculating the equivalent stiffness and damping of the active magnetic bearing, considering the influence of PID parameters, are given by:
K_{\mathrm{amb}} = k_x + k_i \, \mathrm{Re}\left\{
P_1 P_2 P_3
\left(k_p + \frac{k_{ic}}{s}
+ \frac{k_d s}{1+T_0 s}\right)
\left(\frac{1}{b_s+1}\right)
\left(\frac{(1-aT_s)^2}{(1+aT_s)^2}\right)
\right\},
\\[3pt]
C_{\mathrm{amb}} = \frac{k_i}{\omega} \, \mathrm{Im}\left\{
P_1 P_2 P_3
\left(k_p + \frac{k_i}{s}
+ \frac{k_d s}{1+T_0 s}\right)
\left(\frac{1}{b_s+1}\right)
\left(\frac{(1-aT_s)^2}{(1+aT_s)^2}\right)
\right\}.
\end{cases}
\tag{11}$$
where \(P_1\), \(P_2\), and \(P_3\) are the gain coefficients of the sensor, power amplifier, and controller, respectively; \(k_p\), \(k_{ic}\), \(k_d\), and \(T_0\) denote the PID controller parameters; \(b_s\) represents the parameter used in the calculation of the power amplifier; \(a\) is the coefficient of the controller model; \(\omega\) denotes the rotational speed; and \(T_s\) denotes the time delay.
In the dynamic modeling framework, the active magnetic bearing is expressed in matrix form as:
\mathbf{K}_{\mathrm{AMB}}=
\begin{bmatrix}
K_{\mathrm{amb}} & 0 & 0 & 0 & 0 \\
0 & K_{\mathrm{amb}} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0
\end{bmatrix},
\\[10pt]
\mathbf{C}_{\mathrm{AMB}}=
\begin{bmatrix}
C_{\mathrm{amb}} & 0 & 0 & 0 & 0 \\
0 & C_{\mathrm{amb}} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0
\end{bmatrix}.
\end{gathered}
\tag{12}$$
2.3 Modeling Additional Loads in AMB-Rotor Systems
In this study, the rotating unbalance forces of the disk and the thermally induced bending loads of the rotor (resulting from the combined effects of ambient temperature and heat losses in the active magnetic bearing coils) are incorporated as additional forces and moments in the AMB–rotor system.
1. The expressions for the unbalance forces \((\mathbf{Q}_{u})\) of the disk in the X and Y directions are given as:
{Q}_{u}{x} = m_d e_d \omega^2 \cos(\omega t+\varphi) \\
{Q}_{u}{y} = m_d e_d \omega^2 \sin(\omega t+\varphi)
\end{cases},
\tag{13}$$
where \(m_d\) is the disk mass, $e_d\) is the disk eccentricity, and $\omega$ is the rotational angular speed.
2. To determine the thermally induced bending moments under realistic thermal conditions, a three-dimensional steady-state thermal analysis is first performed using finite-element software to obtain the temperature field of the AMB-rotor system. The corresponding thermal stress distribution at each nodal cross-section is subsequently obtained through thermal stress analysis. At a given cross-section, the axial thermal stress \(\sigma_z\) acting on an infinitesimal area \(dA_s\) produces an elemental axial force \(dF_z=\sigma_z dA_s\). With \(b_x\) and \(b_y\) denoting the signed moment arms from the area element to the X and Y axes, respectively, the corresponding elemental bending moments are expressed as follows [17]:
dM_x = b_x\,dF_z \\
dM_y = b_y\,dF_z
\end{cases}.
\tag{14}$$
Integrating these elemental moments over the entire cross-sectional area gives the sectional thermal bending moments \(M_x\) and \(M_y\), as expressed in Eq. (15).
M_x = \int_{A_s}\sigma_z b_x\,dA_s \\
M_y = \int_{A_s}\sigma_z b_y\,dA_s
\end{cases}.
\tag{15}$$
The thermal bending moment values at each nodal cross-section in the rotordynamic model are computed using Eq. (16).
M_x’ = M_y’ = 0 \\
M_x'(n) = M_x(n-1)-M_x(n) \\
M_y'(n) = M_y(n-1)-M_y(n)
\end{cases},
\tag{16}$$
where \(n\) denotes the index of the nodal cross-section along the axial direction of the rotor, which is consistent with the node numbering used in the rotor dynamic model. The resulting equivalent nodal moments are assembled into the global thermal-load vector at the rotational degrees of freedom \(\theta_{x}\) and \(\theta_{y}\), thereby forming the thermally induced load vector (\(\mathbf{Q}_{\mathrm{thermal}}\)).
2.4 Thermally Coupled Dynamic Equations of the AMB-Rotor System
Based on the descriptions provided in Subsections 2.1 to 2.3 and the associated assumptions, together with the rotor-dynamic system modeling principles, the governing equations are first expressed in matrix form for each subsystem and then consolidated into a unified differential equation representing the magnetically suspended rotor system subjected to prescribed thermally induced loads. The resulting governing equation is given by Eq. (17) as follows:
&(\mathbf{M}_{\mathrm{all}}+\mathbf{M}_{\mathrm{disk}})\ddot{\mathbf{X}}
+(\mathbf{C}-\omega(\mathbf{G}_{all}+\mathbf{G}_{disk})
+\mathbf{C}_{AMB})\dot{\mathbf{X}}
\\
&\qquad
+(\mathbf{K}_{all}+\mathbf{K}_{AMB})\mathbf{X}
=
\mathbf{Q}_{weight}
+\mathbf{Q}_{thermal}
+\mathbf{Q}_{u}.
\end{split}.
\tag{17}$$
Here, \(\mathbf{C}\) is the proportional damping matrix and is expressed as \(\mathbf{C}=\eta (\mathbf{M}_{\mathrm{all}}+\mathbf{M}_{\mathrm{disk}})+\beta (\mathbf{K}_{all} +\mathbf{K}_{AMB})\). The coefficient \(\eta\) is calculated by \(\eta = 2 \left( \frac{\xi_2}{\omega_2} – \frac{\xi_1}{\omega_1} \right) / \left( \frac{1}{\omega_2^2} – \frac{1}{\omega_1^2} \right)\), and the coefficient \(\beta\) is determined by \(\beta = \frac{2 (\xi_2 \omega_2 – \xi_1 \omega_1)}{\omega_2^2 – \omega_1^2}\), where \(\xi_{1,2}\) denote the first two damping ratios of the system respectively, while \(\omega_{1,2}\) represent the first two natural frequencies of the system. On the other hand, the matrices \(\mathbf{Q}_{\mathrm{weight}}\), \(\mathbf{Q}_{\mathrm{thermal}}\) and \(\mathbf{Q}_{\mathrm{u}}\) represent the total gravitational load of the shaft-disk assembly, the distributed equivalent thermal bending moments at the nodal shaft cross-sections and the disk-induced unbalance forces, respectively. Meanwhile, \(\ddot{\mathbf{X}}\), \(\dot{\mathbf{X}}\), and \(\mathbf{X}\) represent the acceleration, velocity, and displacement vectors, respectively. The time-domain response of Eq. (17) is obtained using the Newmark-\(\beta\) method with the constant-average-acceleration scheme, for which \(\gamma\) = 0.5 and \(\beta\) = 0.25. The time step is determined as 2\(\pi\)/(1024\(\omega\)).
2.5 Finite Element Calculation of Thermally Induced Bending Moments
A three-dimensional finite element model of the AMB-supported rotor system is constructed using the actual structural dimensions and then preprocessed in the finite element software. The three-dimensional analysis consists of a steady-state thermal analysis followed by thermal-stress analysis. The material properties of the shaft, AMB rotor, stator, and coils are defined accordingly. In particular, the temperature dependence of the shaft properties, including Young’s modulus and Poisson’s ratio, is taken into account. Since the present study focuses on the rotor dynamic response to thermal loading, the temperature dependence of the material properties of the remaining AMB components is neglected. Thermal boundary conditions, including coil heat generation, convective heat transfer coefficients, and ambient temperature, are then applied. The model is discretized using an automatic meshing scheme, with an element size of 4.5 mm for the shaft and 4 mm for the remaining components, as shown in Fig. 3.
Fig. 3. Finite element mesh model of the active magnetic bearing rotor system
After the steady-state thermal analysis is completed, the three-dimensional temperature field distributions under different ambient temperatures are obtained, as shown in Fig. 4. Based on the obtained temperature field, the corresponding thermal stress distribution at each rotor nodal cross-section is derived through thermal-stress analysis. Furthermore, based on the thermal stress results derived from the finite element temperature field, the equivalent thermal bending moments between adjacent nodal cross-sections are calculated and extracted according to Eqs. (14-16). These equivalent thermal bending moments are subsequently assembled into the thermal load vector \(\mathbf{Q}_{\mathrm{thermal}}\) and included in the governing rotordynamic equation.
Fig. 4. Steady-state temperature distributions of the AMB-supported rotor system under different ambient temperatures; A) 100 °C, B) 300 °C, C) 500 °C, and D) 700 °C)
Table 1. Parameters of the Thermally Coupled AMB–Rotor System
| Parameter | Symbol | Value |
|---|---|---|
| Number of rotor nodes | Nnode | 16 |
| Length of each shaft segment (m) | L1 | 1.5×10-2 |
| L2 | 8×10-3 | |
| L3 | 6×10-2 | |
| L4 | 7×10-2 | |
| L5 | 2×10-2 | |
| L6 | 2.7×10-2 | |
| Outer diameter of each shaft section (m) | D1(L1) | 2×10-2 |
| D2(L2) | 2.5×10-2 | |
| D3(L3,L4) | 3.5×10-2 | |
| D4(L5) | 3×10-2 | |
| D5(L6) | 5×10-2 | |
| Young’s modulus of the shaft (20 °C) (GPa) | E | 200 |
| Material density of shaft and AMB rotor (kg/m3) | ρ | 7850 |
| Poisson’s ratio of the shaft and AMB rotor (20 °C) | ν | 0.3 |
| Disk mass (kg) | md | 4.8 |
| Disk diametral mass moment of inertia (kg·m2) | Jd | 1.25×10-2 |
| Disk polar moment of inertia (kg·m2) | Jp | 2.5×10-2 |
| Eccentricity of the disk (m) | ed | 5×10-5 |
| Permeability of free space | μ0 | 1.3×10-6 |
| Magnetic-pole area (m2) | Ab | 2.6×10-4 |
| Number of turns per coil | N | 160 |
| Half-angle between the pole legs of the AMB (°) | α | 22.5 |
| The nominal air gap length of AMB (m) | δ0 | 3×10-4 |
| Sensor gain coefficient | P1 | 1.6×104 |
| Power-amplifier gain coefficient | P2 | 1 |
| Controller coefficient | P3 | 4 |
| PID controller parameters | kp | 0.6 |
| kic | 30 | |
| kd | 2×10-3 | |
| T0 | 1×10-4 | |
| Gravitational acceleration (m/s2) | g | 10 |
| Thermal loss dissipation of AMB at A1, A2 (W/m3) | WA1, WA2 | 4.1×106 |
| Volumetric heat-generation rate A3–A8 (W/m3) | WA3–WA8 | 1.65×105 |
| Air convection heat transfer coefficient (W/m2·°C) | Airc | 15 |
| Overall heat-transfer coefficient of the coil (W/m2·°C) | Coilc | 13.3 |
| Overall heat-transfer coefficient of the rotor (W/m2·°C) | Rotorc | 14 |
| Overall heat-transfer coefficient of the stator (W/m2·°C) | Statorc | 16.7 |
| Density of the coil material (kg/m3) | ρcoil | 8960 |
| Density of the stator material (kg/m3) | ρstator | 7650 |
| Thermal conductivity of the coil material (W/m·°C) | CoilT | 400 |
| Thermal conductivity of stator and rotor (W/m·°C) | RotorT, StatorT | 35 |
| Specific heat capacity of the coil material (J/kg·°C) | Coils | 385 |
| Specific heat capacity of the stator and rotor materials (J/kg·°C) | Rotors, Stators | 450 |
2.6 Model Parameters and Analysis Procedure
This section outlines the model parameters, thermal boundary conditions, and research procedure. The system schematic and overall research flowchart are shown in Figs. 5 and 6, respectively.
Fig. 5. Schematic of the AMB–rotor system; A) three-dimensional configuration, and B) rotor-node discretization
Fig. 6. Research and analysis workflow
The model parameters used in this study are summarized in Table 1. The overall analysis procedure is illustrated in the workflow shown in Fig. 6. Additionally, this study validates the rotor dynamics model and the equivalent stiffness-damping model of the active magnetic bearings, solves the dynamic model using the Newmark-\(\beta\) method, and analyzes the thermally coupled dynamic characteristics of the AMB-supported rotor system. The validation results and parametric analyses are presented in the following section.
3 VALIDATION
To assess the basic accuracy of the developed rotordynamic model, a modal validation is first performed. Since the construction of the proportional damping matrix requires the first two modal frequencies, the modal frequencies predicted by the present model are compared with those obtained from the ANSYS finite element model over a rotational speed range of 0 rpm to 2000 rpm. The comparison results are shown in Fig. 7. The first two modal frequencies predicted by the present model show good agreement with the finite element results in both magnitude and trend, and the maximum relative error of the modal frequencies is 8.8 %. These results support the basic reliability of the developed rotordynamic model and provide a reasonable basis for the subsequent dynamic analyses.
To assess the validity of the equivalent stiffness and damping model used in this study, their variations with rotational frequency were calculated and compared with those reported in [22], as shown in Fig. 8. Since the bearing structural parameters and control parameters used in this study are not identical to those in [22], this comparison is mainly intended for qualitative trend validation.
Fig. 7. Comparison of Campbell diagrams obtained from the present and ANSYS models: FW forward, BW backward
The overall trends predicted by the present model are consistent with the published results. The variation in the equivalent stiffness and damping can be explained as follows. The equivalent stiffness is jointly determined by the displacement stiffness term and the real part of the feedback current stiffness term. Since the displacement stiffness term is negative and the feedback current stiffness term varies with rotational frequency, the equivalent stiffness first increases and then decreases as the rotational frequency increases, and may become negative at higher frequencies. The equivalent damping is mainly governed by the imaginary part of the feedback control term and the rotational frequency. In the low-frequency range, the integral term of the PID controller may lead to negative equivalent damping. As the rotational frequency increases, the equivalent damping rapidly changes from negative to positive. With a further increase in rotational frequency, the system time-delay effect becomes more pronounced, leading to a gradual decrease in the equivalent damping, which eventually approaches zero.
Fig. 8. Equivalent stiffness and damping of the AMB as functions of rotational frequency; A) present model, and B) the model from the literature [22])
Fig. 9. Relative deviations in the average shaft temperature for different mesh sizes (4.5 mm baseline)
Furthermore, a mesh-sensitivity analysis is performed to assess the dependence of the finite element temperature predictions on the shaft mesh size. A representative high-temperature condition with an ambient temperature of 100 °C is considered, while all other boundary conditions and the mesh settings of the remaining components are kept unchanged. The shaft mesh sizes are set to 2.5, 3.5, 4.5, 5.5, and 6.5 mm, and the corresponding average shaft temperature is calculated for each mesh. Using the result obtained with the 4.5 mm mesh as the baseline, the relative deviation of each result from the baseline value is evaluated, as shown in Fig. 9. When the shaft mesh is further refined from 4.5 mm to 3.5 and 2.5 mm, the relative deviations in the average shaft temperature are only 0.25 % and 0.33 %, respectively. These results indicate that the predicted average shaft temperature is insensitive to further mesh refinement below 4.5 mm. Therefore, a shaft mesh size of 4.5 mm is adopted in the subsequent finite element analyses.
4 DISCUSSION
This section examines the effects of bias current, rotational speed, and ambient temperature on the dynamic response of the AMB-supported rotor system. The study focuses on three key aspects: first, the influence of different bias currents on the dynamic behavior of the system under fixed rotational speed and temperature field conditions; second, the effect of varying rotational speeds on the dynamic response of the system under fixed bias current and temperature field; and third, the impact of different ambient temperatures on the dynamic characteristics of the system under fixed bias current and rotational speed. The rotor disk node (Node 8) is selected as the key observation point. Its dynamic behavior under various operating conditions is analyzed using displacement–time history plots, shaft orbit, and frequency spectrum response diagrams.
4.1 Effect of Bias Current on System Dynamics
To investigate the influence of different bias currents on the dynamic characteristics of the system, a baseline operating condition is set with an ambient temperature of 100 °C and a rotational speed of 838 rad/s.
Fig. 10. Displacement responses in X and Y directions at different bias currents; A) X-direction, and B) Y-direction
Fig. 11. Shaft orbit plot at different bias currents
Fig. 12. Frequency spectra in the X and Y directions at different bias currents; A) X-direction, and B) Y-direction
Figures 10–12 indicate that increasing the bias current leads to a moderate reduction in the displacement amplitude and a contraction of the shaft orbit, while the response remains dominated by the fundamental frequency component ($1f$). Neither additional harmonic components nor changes in the motion period are observed. This indicates that the bias current mainly regulates the response magnitude through the support characteristics of the AMB, rather than altering the frequency composition or motion pattern of the rotor system.
4.2 Effect of Rotational Speed on System Dynamics
To investigate the influence of rotational speed on the dynamic characteristics of the system, the bias current is fixed at 3 A and the ambient temperature is set to 100 °C as the baseline operating condition.
Fig. 13. Displacement responses in X and Y directions at different speeds; A) X-direction, and B) Y-direction
Fig. 14. Shaft orbit plot at different speeds
Fig. 15. Frequency spectra in the X and Y directions at different speeds; A) X-direction, and B) Y-direction
Figures 13–15 show that rotational speed primarily governs the synchronous vibration component. As the rotational speed increases, both the amplitude of the fundamental (\(1f\)) response and the size of the shaft orbit increase. In particular, the response increases more markedly between 1047 rad/s and 1257 rad/s, indicating enhanced sensitivity near the upper end of the investigated speed range. This behavior results from the combined effects of the speed-dependent unbalance excitation and the gyroscopic matrix term in the governing equations, which jointly influence the vibration response at the system nodes. Nevertheless, the system retains a single-periodic motion throughout the investigated speed range. As the rotational speed increases, the fundamental-frequency peaks in both the X and Y directions shift to higher frequencies and continue to track the corresponding rotational frequency.
Fig. 16. Displacement responses in X and Y directions at different ambient temperatures; A) X-direction, and B) Y-direction
4.3 Effect of Ambient Temperature on System Dynamics
To investigate the influence of different ambient temperatures on the dynamic characteristics of the system, the rotational speed is fixed at 838 rad/s and the bias current is set to 3 A.
Fig. 17. Shaft orbit plot at different ambient temperatures
Fig. 18. Frequency spectra in the X and Y directions at different ambient temperatures; A) X-direction, and B) Y-direction
Figures 16–18 show that the predominant effect of increasing ambient temperature is a pronounced migration of the shaft-orbit center, rather than the emergence of additional frequency components. As the ambient temperature increases, the mean displacement offsets in both directions increase markedly, and the shaft-orbit center shifts toward the upper-right region. Meanwhile, the system response remains dominated by the synchronous fundamental frequency and retains a single-periodic motion pattern. These observations indicate that the thermally induced bending moments act primarily as quasi-steady loads. As the ambient temperature increases, the overall temperature of the rotor system rises, thereby increasing the bending moments generated by the nonuniform temperature field and causing pronounced changes in the dynamic response.
5 CONCLUSIONS
This study focuses on the thermally induced dynamic behavior of an active magnetic bearing (AMB)-supported rotor system operating under high-temperature conditions. A one-way sequential thermal–mechanical–rotordynamic analysis framework is established by considering the combined effects of AMB heat generation and ambient thermal loading. The model incorporates thermally induced loads derived from the temperature field, disk mass unbalance, gravity, and the equivalent support stiffness and damping of the AMB governed by PID control parameters. Based on the developed model, the dynamic response characteristics of the system under varying bias currents, rotational speeds, and ambient temperatures are systematically investigated. The main conclusions are summarized as follows:
1. Increasing the bias current alters the effective support characteristics of the active magnetic bearing and reduces the displacement amplitude and gradually contracts the shaft orbit. Nevertheless, the spectral content of the response and its single-periodic motion pattern remain essentially unchanged, indicating that the bias current mainly regulates the response magnitude without altering the spectral composition or response periodicity.
2. Rotational speed plays a dominant role in governing the synchronous vibration response. As the speed increases, both the response amplitude and the spatial extent of the shaft orbit increase progressively. The more pronounced growth observed between 1047 rad/s and 1257 rad/s further indicates enhanced response sensitivity in the upper portion of the investigated speed range.
3. Ambient temperature predominantly influences the displacement response and the location of the shaft-orbit center, while its effect on the spectral composition remains comparatively limited. The pronounced shaft-orbit center migration is mainly associated with the thermal bending moments induced by the temperature field. Throughout the investigated temperature range, the rotor maintains a single-periodic motion state, with the response consistently dominated by the fundamental frequency.
References
[1] Bleuler, H. A survey of magnetic levitation and magnetic bearing types. JSME Int J Ser III 35 335–342 (1992) DOI:10.1299/jsmec1988.35.335.
[2] Huang, Z., Li, C., Zhou, Z., Liu, B., Zhang, Y., Yang, M. et al. Magnetic bearing: Structure, model, and control strategy. Int J Adv Manuf Technol 131 3287–3333 (2024) DOI:10.1007/s00170-023-12389-8.
[3] Liu, C., Zhan, J., Yang, Y., Liu, Z. Review of research status and development of flexible rotor–magnetic bearing. Proc CSEE 40 4602–4614, 4739 (2020).
[4] Gao, R., Luo, G., Yan, C. Dynamic modeling and analysis of active magnetic bearings. Appl Mech Mater 494–495 685–688 (2014) DOI:10.4028/www.scientific.net/AMM.494-495.685.
[5] Wu, C., Su, Z., Wang, D., Jiang, H. Dynamic modeling method for active magnetic bearings–rotor system of steam turbines. J Mech Sci Technol 37 1665–1673 (2023) DOI:10.1007/s12206-023-0308-x.
[6] Ranjan, G., Narsakka, J., Choudhury, T., Nevaranta, N., Sopanen, J. Method for evaluating stresses and fatigue life of a coupled multi-shaft rotor supported by active magnetic bearings in dropdown events. Nonlinear Dyn 113 17705–17724 (2025) DOI:10.1007/s11071-025-11109-z.
[7] de Jongh, F. The synchronous rotor instability phenomenon—morton effect. Proc. 47th Turbomachinery Symp. Turbomachinery Laboratory, Texas A&M Engineering Experiment Station (2018), Houston, USA.
[8] Badykov, R., Liu, Z., Yurtaev, A., Diligenskiy, D., Zhou, W., Molokanov, A. Magnetostatic and thermal analysis of classical and hybrid axial active magnetic bearings. Actuators 15 30 (2026) DOI:10.3390/act15010030.
[9] Takahashi, N., Kaneko, S. Thermal instability in a magnetically levitated doubly overhung rotor. J Sound Vib 332 1188–1203 (2013) DOI:10.1016/j.jsv.2012.09.001.
[10] Xu, Y., Li, Z., Xing, D. Multi-physics coupling dynamics simulation of thermally induced vibration of magnetically suspended rotor in small and micro nuclear reactors. Energies 18 2433 (2025) DOI:10.3390/en18102433.
[11] Jin, C., Su, H., Dong, Y., Zhou, J., Xu, Y., Yan, X. Iron loss analysis of magnetic bearing system considering magnetic–thermal coupling. Int J Appl Electromagn Mech 70 1–20 (2022) DOI:10.3233/JAE-220011.
[12] Zhao, Y., Zhou, J., Guo, M., Xu, Y. A thermal flexible rotor dynamic modelling for rapid prediction of thermo-elastic coupling vibration characteristics in non-uniform temperature fields. Appl Math Model 138 115751 (2025) DOI:10.1016/j.apm.2024.115751.
[13] Cui, Y.S., Wang, Y.Q. Free vibrations of axially loaded thin-walled shaft-disk rotors subjected to non-uniform temperature field. Thin-Walled Struct 196 111461 (2024) DOI:10.1016/j.tws.2023.111461.
[14] Chang, Z., Hou, L., Chen, Y. Nonlinear dynamics and thermal bidirectional coupling characteristics of a rotor–ball bearing system. Appl Math Model 119 513–533 (2023) DOI:10.1016/j.apm.2023.03.009.
[15] Chang, Z., Hou, L., Masarati, P., Lin, R., Li, Z., Chen, Y. Modeling and nonlinear analysis of a coupled thermo-mechanical dual-rotor system. Nonlinear Dyn 112 17811–17842 (2024) DOI:10.1007/s11071-024-09969-y.
[16] Yu, X., Liu, Z., He, P., Qi, J. Experimental research on thermal bow deformation and dynamics of aero-engine rotor. J Aerosp Power 40 20230188 (2025) DOI:10.13224/j.cnki.jasp.20230188.
[17] Goldman, P., Muszynska, A. Rotor-to-stator, rub-related, thermal/mechanical effects in rotating machinery. Chaos Solitons Fractals 5 1579–1601 (1995) DOI:10.1016/0960-0779(94)00165-M.
[18] Larsson, B. Journal asymmetric heating—part I: Nonstationary bow. J Tribol 121 157–163 (1999) DOI:10.1115/1.2833811.
[19] Larsson, B. Journal asymmetric heating—part II: Alteration of rotor dynamic properties. J Tribol 121 164–168 (1999) DOI:10.1115/1.2833798.
[20] Du, B., Xiang, H., Ge, Y., Zhu, L. Derivation and application of 3D-beam’s element stiffness and mass matrix with shear effect. J Chongqing Jiaotong Univ Nat Sci 27 502–507 (2008).
[21] Jin, C., Xu, Y., Zhou, J., Cheng, C. Active magnetic bearings stiffness and damping identification from frequency characteristics of control system. Shock Vib 2016 1067506 (2016) DOI:10.1155/2016/1067506.
[22] Xu, Y. Identification of Supporting Characteristics of Flexible Rotor–AMB System. PhD dissertation, Nanjing University of Aeronautics and Astronautics, Nanjing, China (2018).