Modeling and Robust Control of Magnetic Bearing-Rotor System Considering Interface Contact
-
摘要:
磁悬浮轴承-转子系统中，螺栓连接装配形成的界面接触会在转子悬浮状态下激发弯曲模态振动，同时振动频率随转速明显变化. 为实现全转速范围内对弯曲模态振动的主动控制，提出一种考虑频率不确定性的H∞鲁棒控制器设计方法. 首先，通过建立考虑界面接触的动力学模型进行仿真预测，获得振动频率的变化范围；其次，通过频响拟合的方式对转子传递函数进行重构，并将仿真得到的振动频率变化范围以加性不确定性方式引入到重构传递函数中，得到考虑模态频率不确定性的被控对象模型；最后，基于该模型设计兼顾参数摄动和外力扰动的鲁棒性、闭环系统稳定性和防止控制电压饱和等多功能的H∞ 鲁棒控制器. 数值仿真结果表明：该控制器在模态频率处具有宽频带阻的频响特性，能够抑制磁悬浮轴承转子系统的弯曲模态振动；采用该方法设计的H∞ 鲁棒控制器后，转子弯曲模态振动幅值最少减小90%.
Abstract:In magnetic bearing-rotor systems, the bending mode vibration may be excited by the interface contact formed by bolt joints during rotor levitation, and the vibration frequency varies with rotation speed. To actively control bending mode vibration at any speed, the design method of a robust H∞ controller considering frequency uncertainty was proposed. Firstly, the dynamic model considering interface contact was established for numerical simulation, and the vibration frequency variation was obtained. Then, the rotor transfer function was reconstructed by frequency response fitting, and the variation range of vibration frequency obtained by simulation was introduced into the reconstructed transfer function by means of additive uncertainty. As a result, a controlled object model considering mode frequency uncertainty was obtained. Finally, based on the model, the robust H∞ controller was designed by taking the robustness to parameter perturbations and external disturbance, closed-loop system stability, control voltage saturation, and other functions into account. The numerical simulation results show that the controller has the frequency response characteristic of wide band resistance at the mode frequency, which is able to suppress the bending mode vibration of the magnetic bearing-rotor system. After the robust H∞ controller designed by this method is used, the bending mode vibration amplitude of the rotor is reduced by more than 90%.
-
Key words:
- active magnetic bearing /
- interface contact /
- mode vibration /
- robust control
-
平面电机式磁悬浮转台具有多自由度平移、旋转运动能力[1]，且克服了传统机械传动平台存在的接触摩擦、回程误差等限制因素[2]，在微细加工、复杂微结构成型等精密制造应用场景中具有极大优势. 然而，由于采用Halbach永磁阵列，磁悬浮转台气隙磁场的空间周期性使得输出磁力大小和方向随运动次级的位置变化而变化，且在不同运动自由度之间存在非线性的磁力耦合. 此外，磁悬浮系统固有的开环不稳定性使得工作台易受扰动影响，对运动控制算法的稳定性和抗扰性要求较高.
为实现磁悬浮工作台对指定轨迹的快速精准跟踪，研究人员基于PID（proportion integration differentiation）控制器对其进行控制器设计[3-4]. 但PID控制难以兼顾动态性能和运动精度，如增大微分增益会导致在提升响应速度的同时放大测量噪声，对进一步提高跟踪精度具有局限性. 此外，建模和测量的不精确以及外界扰动导致系统存在不确定性，易于引起控制系统品质恶化，从而导致运动精度变低. 滑模控制在系统存在参数摄动和不确定性扰动时表现优异[5]，魏静波等[6]结合终端滑模干扰观测器和全局快速终端滑模控制优化了磁悬浮球系统的控制效果. 此外，采用分数幂函数代替原有的符号函数，利用非线性部分辅助线性部分削弱切换项带来的抖振问题[7-8]. 近年来，具有全局记忆性的分数阶微积分被引入滑模控制策略[9]，相较于传统的整数阶滑模控制更具优越性. Sun等[10]针对直线电机提出一种离散时间分数阶终端滑模控制策略，提升了其运动跟踪性能；Kuang等[11]利用交叉耦合二阶分数阶滑模控制方法提高了直线电机的同步控制和跟踪精度，对磁悬浮转台系统运动控制方法设计具有借鉴意义. 此外，相比于直接针对系统的非线性模型进行控制器设计[12]，采用动态解耦方法将六自由度运动模型转换为6个二阶线性模型可简化控制器设计过程[13-14]. 为充分考虑系统的耦合效应，采用干扰观测器进行扰动补偿[15].
为抑制系统非线性、强耦合性和不确定性对磁悬浮转台运动精度的影响，首先，采用非线性干扰观测器对包含模型误差、测量误差、六自由度间耦合项和其他外界干扰在内的系统集总扰动进行补偿；然后，设计了一种离散时间分数阶滑模控制策略，分数幂函数削弱抖振，分数阶微积分对跟踪误差的全局记忆性来提升系统的动态性能和轨迹跟踪性能；最后，通过实验验证了所提方法的优越性和实用性.
1. 磁悬浮转台动力学分析
图1中的磁悬浮转台具有中心对称结构，主体部分由8个矩形线圈绕组（定子，Coil 1～Coil 8）和1个搭载圆周型Halbach永磁阵列的圆台（动子）构成，适合执行平面运动和绕竖直方向的大范围旋转运动. 以线圈阵列上表面的几何中心为坐标原点$ O $，通过设计控制器调节线圈通电电流，可使线圈长边区域的电流与Halbach永磁阵列产生的空间磁场相互作用，产生可控的磁力和磁力矩，从而驱动磁悬浮转台动子在$ x、y、{\textit{z}}、\alpha 、\beta 、\gamma $方向上进行六自由度运动.
1.1 电磁力模型
圆周型Halbach永磁阵列由12个磁极组成，其中每个磁极包含4块永磁体. 在动子坐标系{m}中，以Halbach永磁阵列下表面的几何中心为原点$ {}^{\text{m}}O $，将该Halbach永磁阵列沿其切向方向展开，则其下方空间的磁通密度分布如图2所示，图中：$ {}^{\text{m}}r $、$ {}^{\text{m}}\theta $、$ {}^{\text{m}}{\textit{z}} $分别为坐标系{m}的3个坐标变量，$ {}^{\text{m}}\theta \in [0,2{\text{π )}} $；wm、hm分别为永磁体宽和高；ϕm为相邻永磁体间夹角； r为永磁阵列中心半径；$ \tau $为永磁阵列周期角.
忽略高次谐波和边端效应，Halbach永磁阵列下方空间的磁通密度$ {}^{\text{m}}{\boldsymbol{B}} $可表示为[16]
mB=[mBrmBθmBz]=[0−4Brπsin(πwmτ)(1−e−wm/λc)sin(πmθ2ϕm)emz/λc4Brπsin(πwmτ)(1−e−wm/λc)cos(πmθ2ϕm)emz/λc]， (1) 式中：$ {}^{\text{m}}{B_r} $、$ {}^{\text{m}}{B_\theta } $、$ {}^{\text{m}}{B_{\textit{z}}} $分别为$ {}^{\text{m}}B $在$ {}^{\text{m}}r $、$ {}^{\text{m}}\theta $、$ {}^{\text{m}}{\textit{z}} $方向上的分量，${B_{\rm{r}}}$为永磁体剩磁；$ {\lambda _{\text{c}}} = \tau /2{\text{π }} $，为特征波长.
根据洛伦兹积分公式和牛顿第三定律，动子受到的驱动力为通电线圈所受电磁力的反作用力，即
{Fi=−R∫Vcoil iJi×mBdν,Ti=−R∫Vcoil iri×Ji×mBdν, (2) 式中： $ {{\boldsymbol{F}}_i},{{\boldsymbol{T}}_i} \in {\mathbb{R}^3} $，分别为第$i$相线圈提供的磁力和磁力矩，i = 1，2，…，8，$ {\boldsymbol{R}} $为不同坐标系间的转换矩阵；$ {V_{{\text{coil }}i}} $为第$i$相线圈所占的空间，为线圈的有效作用区域；$ {\text{d}}\nu = {}^{\text{m}}r {\text{d}}{}^{\text{m}}r {\text{d}}{}^{\text{m}}\theta {\text{d}}{}^{\text{m}}{\textit{z}} $；$ {{\boldsymbol{J}}_i} $为线圈体积元内的电流密度，方向为通入的电流方向，大小为$|{{\boldsymbol{J}}_i}| = \dfrac{{N{I_i}}}{{{r_{\text{c}}}{h_{\text{c}}}}}$， $ N $为线圈匝数，$ {I_i} $为驱动电流，${r_{\text{c}}}$为线圈侧边宽度，${h_{\text{c}}}$为线圈厚度；${{\boldsymbol{r}}_i}$为线圈体积元到点$ {}^{\text{m}}O $的力臂. 此外，线圈的长边宽度和短边宽度分别为${l_{\text{c}}}$和${w_{\text{c}}}$.
可采用磁节点和高斯积分方法对式（2）进行计算[1]. 将${{\boldsymbol{F}}_i}、{{\boldsymbol{T}}_i}$与其驱动电流$ {I_i} $的关系式简记为
{Fi=[Fi,xFi,yFi,z]T=[Γi,xΓi,yΓi,z]TIi,Ti=[Ti,αTi,βTi,γ]T=[Γi,αΓi,βΓi,γ]TIi， (3) 式中：$ {F_{i,x}},{F_{i,y}},{F_{i,{\textit{z}}}} \in \mathbb{R} $，分别为第$i$相线圈提供的磁力在平动自由度$ x $、$ y $、$ {\textit{z}} $的分量，$ {\varGamma _{i,x}},{\varGamma _{i,y}},{\varGamma _{i,{\textit{z}}}} \in \mathbb{R} $，为相应的磁力系数；$ {T_{i,\alpha }},{T_{i,\beta }},{T_{i,\gamma }} \in \mathbb{R} $，分别为第$i$相线圈提供的磁力矩在转动自由度$ \alpha $、$ \beta $、$ \gamma $的分量，$ {\varGamma _{i,\alpha }}, {\varGamma _{i,\beta }},{\varGamma _{i,\gamma }} \in \mathbb{R} $，为对应的磁力矩系数，与动子的六自由度位移有关.
综合8相线圈提供的磁力和磁力矩，可得磁悬浮转台动子所受到的合驱动力F和力矩T为
w=[FT]=[Γ1,xΓ2,x⋯Γ8,xΓ1,yΓ2,y⋯Γ8,yΓ1,zΓ2,z⋯Γ8,zΓ1,αΓ2,α⋯Γ8,αΓ1,βΓ2,β⋯Γ8,βΓ1,γΓ2,γ⋯Γ8,γ][I1I2⋮I8]=ΓI, (4) 式中：$ {\boldsymbol{w}} \in {\mathbb{R}^6} $，为由动子所受的合磁力和磁力矩构成的向量；$ {\boldsymbol{\varGamma }} \in {\mathbb{R}^{6 \times 8}} $，为系统的电流-动力传递矩阵；${\boldsymbol{I}} \in {\mathbb{R}^8}$，为通入各相线圈的电流构成的向量.
1.2 系统动力学
在磁悬浮转台系统设计中，通过采用合理的电流分配方法，可以实现不同运动自由度上的力和力矩的动态解耦[13]，如式（5）.
I=Γ†w=ΓT(ΓΓT)−1w, (5) 式中：$ {{\boldsymbol{\varGamma }}^\dagger } $为$ {\boldsymbol{\varGamma }} $的伪逆.
记磁悬浮转台相对于初始状态的六自由度位移$ {\boldsymbol{x}} = {\left[ {xmymzmαmβmγm} \right]^{\text{T}}} $，对动子进行动力学分析可得
¨x=J−1(w−[00mg000]T)+d=J−1u+d, (6) 式中：$ {\boldsymbol{J}} = {\text{diag }}(m,{\text{ }}m,{\text{ }}m,{\text{ }}{J_\alpha },{\text{ }}{J_\beta },{\text{ }}{J_\gamma }) $，为动子的惯性矩阵，$m$为动子质量， $ {J_\alpha },{\text{ }}{J_\beta }, {\text{ }}{J_\gamma } \in \mathbb{R} $，分别为转动自由度$ \alpha $、$ \beta $、$ \gamma $的转动惯量；$ {\boldsymbol{u}} = {\boldsymbol{w}} - {[0\;\;0\;\;mg\;\;0\;\;0\;\;0]^{\mathrm{T}}} $；$ {\boldsymbol{d}} \in {\mathbb{R}^6} $，为系统集总扰动，包括模型误差、测量误差、六自由度间耦合项和其他外界干扰.
因此，可将$ {\boldsymbol{u}} $作为控制量设计闭环控制算法，根据式（5）得到电流向量$ {\boldsymbol{I}} $，从而调节线圈电流大小以实现磁悬浮转台运动控制.
假设1 扰动${\boldsymbol{d}}$缓慢变化，即$ ||{\boldsymbol{\dot d}}|{|_2} \leqslant \varepsilon $，$ \varepsilon \approx 0 $，其中||·||2表示二范数，$ \varepsilon $为${\boldsymbol{d}}$的二范数的上界.
由于系统实际工作环境较为稳定且能量有限，且工作台稳态运动时各自由度间的耦合项也随时间缓慢变化，故假设1合理.
设$ {{\boldsymbol{x}}_{\text{d}}} \in {\mathbb{R}^6} $，为磁悬浮转台的期望轨迹，$ {\boldsymbol{e}} = {\boldsymbol{x}} - {{\boldsymbol{x}}_{\text{d}}} = [ e1e2⋯en ]{{\text{ }}^{\text{T}}} $，为跟踪误差，则跟踪系统的动力学方程为
¨e=J−1u+d−¨xd. (7) 采用前向欧拉法将式（7）离散化，可得跟踪误差的离散域状态方程为
Δ2e(k)=J−1u(k)+d(k)−Δ2xd(k), (8) 式中：$ \Delta $为前向差分算子；$ \Delta {{\boldsymbol{x}}_{\text{d}}}(k) = \dfrac{{{{\boldsymbol{x}}_{\text{d}}}(k + 1) - {{\boldsymbol{x}}_{\text{d}}}(k)}}{h} $， $ {\Delta }^{{n}_{ + } + 1}{{\boldsymbol{x}}}_{\text{d}}(k)=\Delta ({\Delta }^{{n}_{ + }}{{\boldsymbol{x}}}_{\text{d}}(k))，{n}_{ + }\in {\mathbb{N}}^{ + } $，$ {{\boldsymbol{x}}_{\text{d}}}(k) = {{\boldsymbol{x}}_{\text{d}}}(kh) $，$h \in {\mathbb{R}^ + }$，为采样间隔，且满足$ 0 < h \ll 1{\text{ s}} $，$ k \in \mathbb{N} $，为采样步数.
2. 分数阶滑模控制器设计
2.1 非线性干扰观测器
采用非线性干扰观测器估计系统集总扰动${\boldsymbol{d}}$，以实现干扰补偿. 设计非线性干扰观测器为[15]
{ˆd=r+p(e),˙r=− L0(r+p(e)+J−1u−¨xd), (9) 式中：$ {\boldsymbol{\hat d}} \in {\mathbb{R}^6} $，为扰动估计值；$ {\boldsymbol{r}}，\dot{{\boldsymbol{r}}}\in {\mathbb{R}}^{6} $，为干扰观测器状态；$ {{\boldsymbol{L}}_0} \in {\mathbb{R}^{6 \times 6}} $，为观测器增益，是一个正定的常数矩阵；${\boldsymbol{p}}( {\text{•}} )$为非线性函数，且$\dot {\boldsymbol{p}}({\boldsymbol{e}}) = {{\boldsymbol{L}}_0}{\boldsymbol{\ddot e}}$.
由式（9）可知，扰动估计值的变化率为
˙ˆd=˙r+˙p(e)=−L0(r+p(e)+J−1u−¨xd)+L0¨e=L0[(¨e−J−1u+¨xd)−(r+p(e))]=L0(d−ˆd)=L0˜d， (10) 式中：$ {\boldsymbol{\tilde d}} = {\boldsymbol{d}} - {\boldsymbol{\hat d}} $，为扰动估计误差，其变化率如式（11）.
˙˜d=˙d−˙ˆd=˙d−L0˜d. (11) 定理1 在假设1的条件下，利用式（9）估计${\boldsymbol{d}}$，可以使估计误差$ {\boldsymbol{\tilde d}} $有界，且其上界可调节到任意小.
证明 构造Lyapunov函数为
U=12˜dT˜d. (12) 设$ {{\boldsymbol{L}}_0} $的最小特征值为$ {\lambda _{\min }}({{\boldsymbol{L}}_0}) $，$ \sigma $为任意小正常数且满足$ \varsigma = {\lambda _{\min }}({{\boldsymbol{L}}_0}) - \sigma > 0 $，使得杨氏不等式$ {{\boldsymbol{\tilde d}}^{\text{T}}}{\boldsymbol{\dot d}} \leqslant \sigma {{\boldsymbol{\tilde d}}^{\text{T}}}{\boldsymbol{\tilde d}} + \dfrac{1}{{4\sigma }}{{\boldsymbol{\dot d}}^{\text{T}}}{\boldsymbol{\dot d}} $成立. 则
˙U=˜dT˙˜d=˜dT(˙d−L0˜d)⩽−˜dTL0˜d+σ˜dT˜d+14σ˙dT˙d⩽−[λmin(L0)−σ]˜dT˜d+ε24σ=−2ςU+ε24σ， (13) 解得
U(t)⩽ε28σς+[U(0)−ε28σς]e−2ςt. (14) 由式（12）可得
||˜d||2⩽√ε24σς+[2U(0)−ε24σς]e−2ςt. (15) 因此，$ ||{\boldsymbol{\tilde d}}|{|_2} $最终将收敛到以$ \sqrt {{\varepsilon ^2}/4\sigma \varsigma } $为半径的球域内，即所设计的非线性干扰观测器对系统集总扰动${\boldsymbol{d}}$的估计误差$ {\boldsymbol{\tilde d}} $是有界的. 通过适当地选择$ {{\boldsymbol{L}}_0} $和$ \sigma $，可以使其上界任意小.
显然，采用式（10）中的干扰观测器，对于离散域中$ {\boldsymbol{\tilde d}}(k) $的每个分量$ {\tilde d_n}(k) $，存在常数$ {d^*} \in {\mathbb{R}^ + } $使得$ {\tilde d_n}(k) $的绝对值$ |{\tilde{d}}_{n}(k)|\leqslant {d}^{*}，k\in \mathbb{N} $，n = 1，2，…，6，即
−d∗⩽˜dn(k)⩽d∗. (16) 2.2 离散时间分数阶滑模控制器
定义1 非光滑分数幂函数定义为
[[e(k)]]q≜ |e(k)|qsgn(e(k)), (17) 式中：$ q $为正奇数之比， $ 0 < q < 1 $.
对于$ {\boldsymbol{e}}(k) = [ {{e_1}(k)}\;\,{{e_2}(k)}\;\, \cdots \;\,{{e_n}(k)} ]^{\text{T}} \in {\mathbb{R}^n} $，$ { [\kern-0.15em[ {\boldsymbol{e}}(k) ]\kern-0.15em] ^q} = [ {|e_1|^q{{\mathrm{sgn}}} ({e_1}(k))}\quad{|e_2|^q{{\mathrm{sgn}}} ({e_2}(k))\quad \cdots }\quad{|e_n|^q{{\mathrm{sgn}}} ({e_n}(k))} ]^{\text{T}} $.
注1 若$ {e_n}(k) > 0 $，则$ { [\kern-0.15em[ {e_n}(k) ]\kern-0.15em] ^q}{\text{ = }}e_n^q(k) $.
为削弱切换项带来的抖振问题，采用式（17）代替符号函数[7]，设计整数阶离散时间切换函数为
s(k)=Δe(k)+l1e(k)+l2[[e(k)]]q, (18) 式中：$ {\boldsymbol{s}}(k) \in {\mathbb{R}^6} $，为滑模状态，$ {l_1} $、$ {l_2} $为控制器参数，$ 0 < {l_1}h < 1 $，$ {l_2} > 0 $.
定义2 G-L形式的分数阶算子定义为
Δae(k)=1hak∑j=0(−1)j(aj)e(k−j), (19) 式中：$ a \in \mathbb{R} $，$ 0 < a < 1 $，$h \in {\mathbb{R}^ + }$，二项式系数$ \left( {aj} \right) $如式（20）.
(aj)={1,j=0,a(a−1)⋯(a−j+1)j!,j=1,2,3,⋯ (20) 对于$ {\boldsymbol{e}}(k) \in {\mathbb{R}^n} $，$ {\Delta ^a}{\boldsymbol{e}}(k) = [ {{\Delta ^a}{e_1}(k)}\quad{{\Delta ^a}{e_2}(k) \quad \cdots }\quad {{\Delta ^a}{e_n}(k)} ]^{\text{T}} $，$ {\Delta ^a}{\boldsymbol{e}}(k) = \Delta ({\Delta ^{a - 1}}{\boldsymbol{e}}(k)) $.
注2 在实际系统中，设$ L $为累加次数，简化计算式为
Δae(k)=1haL∑j=0(−1)j(aj)e(k−j)，L∈Z+. (21) 为减小跟踪误差，可以利用分数阶微积分的全局记忆性，将式（18）改进为
s(k)=Δe(k)+l1e(k)+l2Δa−1[[e(k)]]q. (22) 由式（8）可推出
s(k+1)=s(k)+hΔs(k)=s(k)+h(Δ2e(k)+l1Δe(k)+l2Δa[[e(k)]]q)=s(k)+h(J−1u(k)+d(k)−Δ2xd(k)+l1Δe(k)+l2Δa[[e(k)]]q). (23) 提出一种由指数到达律和等效控制律组成的滑模控制律. 其中，指数到达律表达式为
ure(k)=−J(k1s(k)+k2[[s(k)]]b), (24) 式中：$ {k_1} $、$ {k_2} $、$ b $为控制器参数，$ 0 < {k_1}h < 1 $，$ {k_2} > 0 $，$ b $为正奇数之比，$ 0 < b < 1 $.
令$ {\boldsymbol{s}}(k + 1) - {\boldsymbol{s}}(k) = {\boldsymbol{0}} $，并用扰动估计值$ {\boldsymbol{\hat d}}(k) $代替$ {\boldsymbol{d}}(k) $，可得等效控制律为
ueq(k)=J(Δ2xd(k)−l1Δe(k)−l2Δa[[e(k)]]q−ˆd(k)). (25) 因此，总的滑模控制律为
u(k)=ure(k)+ueq(k). (26) 将式（26）代入式（23）中，有
s(k+1)=s(k)−k1hs(k)−k2h[[s(k)]]b+h˜d(k). (27) 磁悬浮转台闭环系统控制框图如图3所示.
3. 稳定性分析
3.1 系统稳定性
磁悬浮转台闭环系统稳定性可由定理2构造.
定理2 对于式（8），若选择滑模状态$ {\boldsymbol{s}}(k) $的动态方程为式（22），选择干扰观测器为式（9），那么可以采用式（26）描述的分数阶滑模控制律来保证以下2个条件：
1） 当系统状态远离滑模面时，每一个滑模变量$ {s_n}(k) $，将在$ {k^ * } $步内进入区域$ \varLambda $，如式（28）、（29）.
{Λ=(sn(k):|sn(k)|⩽δ=ψ(b)max{(d∗k2)1b,(k2h1−k1h)11−b})，ψ(b)=1+bb1−b−b11−b, (28) {k∗=⌊s2n(0)−δ2μ2⌋+1，μ=k1hδ+(ψb(b)−1)d∗h, (29) 式中：$ \mu $为中间变量，$ \delta $为$ {s_n}(k) $的上界，$ \psi ({\text{•}}) $为自定义函数.
2） 滑模变量$ {s_n}(k) $进入区域${\rm{\varLambda}}$之后将始终保持在该区域中，即若$ |{s_n}(k)| \leqslant \delta $，则$ |{s_n}(k + 1)| \leqslant \delta $.
证明 证明分为到达段和滑动段两部分.
步骤1 （到达段） 首先证明当${s_n}(k) \notin \varLambda $，即$ |{s_n}(k)| > \delta $时，滑模状态将在$ {k^ * } $步内进入区域${\rm{\varLambda}}$.
在离散域构造Lyapunov函数为
V(k)=s2n(k). (30) V(k+1)−V(k)=s2n(k+1)−s2n(k)=−(k1hsn(k)+k2h[[sn(k)]]b−h˜dn(k))⋅(2s2n(k)−k1hs2n(k)−k2h[[sn(k)]]b+h˜dn(k)). (31) 引理1[7] 若$ \psi (\zeta ) = 1 + {\zeta ^{\frac{\zeta }{{1 - \zeta }}}} - {\zeta ^{\frac{1}{{1 - \zeta }}}} $，$ 0 < \zeta < 1 $，则有$ 1 < \psi (\zeta ) < 2 $，且对于任意$ 0 \leqslant \beta \leqslant 1 $有
βψ(ζ)−βζψζ(ζ)+ψ(ζ)−1⩾0. (32) 当$ |{s_n}(k)| > \delta $时，考虑以下2种情况：
情况1 $ {s_n}(k) < - \delta $.
由$ {s_n}(k) < - \delta = - \psi (b)\max \left\{ {{{\left( {\dfrac{{{d^ * }}}{{{k_2}}}} \right)}^{\frac{1}{b}}},{\text{ }}{{\left( {\dfrac{{{k_2}h}}{{1 - {k_1}h}}} \right)}^{\frac{1}{{1 - b}}}}} \right\} $得
{[[sn(k)]]b<−δb⩽−ψb(b)d∗k2<0[[sn(k)]]1−b>δ1−b⩾ψ1−b(b)(k2h1−k1h)>0. (33) 则有
{k2[[sn(k)]]b<−d∗ψb(b),(1−k1h)sn(k)<ψ1−b(b)k2h[[sn(k)]]b. (34) 由于$ - {\tilde d_n}(k) < {d^ * } $，由式（34）可得
k1hsn(k)+k2h[[sn(k)]]b−h˜dn(k)<−k1hδ−d∗hψb(b)+d∗h=−(k1hδ+(ψb(b)−1)d∗h)=−μ. (35) 由引理1可知，$ 1 < \psi (b) < 2 $，所以$ \mu $为正数.
由于$ {\psi ^{1 - b}}(b) > 1 $，$ { [\kern-0.15em[ {s_n}(k) ]\kern-0.15em] ^b} < 0 $，由式（34）可得
sn(k)<k1hsn(k)+ψ1−b(b)k2h[[sn(k)]]b<k1hsn(k)+k2h[[sn(k)]]b. (36) 由于$ {\tilde d_n}(k) < {d^ * } $，根据式（36）、（34）可推导出
2sn(k)−k1hsn(k)−k2h[[sn(k)]]b+h˜dn(k)<2(k1hsn(k)+k2h[[sn(k)]]b)−k1hsn(k)−k2h[[sn(k)]]b+d∗h=k1hsn(k)+k2h[[sn(k)]]b+d∗h<−k1hδ−d∗hψb(b)+d∗h=−(k1hδ+(ψb(b)−1)d∗h)=−μ. (37) 因此，有$ V(k + 1) - V(k) < - {\mu ^2} < 0 $.
情况2 当$ {s_n}(k) > \delta $时，同理可证明$ V(k + 1) - V(k) < - {\mu ^2} < 0 $.
综上所述，当$ {s_n}(k) < - \delta < 0 $或$ {s_n}(k) > \delta > 0 $时，即$ {s_n}(k) \notin \varLambda $时，均有$ V(k + 1) - V(k) < - {\mu ^2} < 0 $成立.
记滑模变量的初始值为$ {s_n}(0) $，由式（30）得
s2n(k)<s2n(k−1)−μ2<⋯<s2n(0)−kμ2. (38) 令$ s_n^2(0) - k{\mu ^2} = {\delta ^2} $，解得$ k = (s_n^2(0) - {\delta ^2})/{\mu ^2} $. 即系统状态远离滑模面时可在$ {k^ * } $步内进入区域${\rm{\varLambda}}$.
步骤2 （滑动段） 证明系统状态一旦进入区域${\rm{\varLambda}}$，之后将始终保持在区域${\rm{\varLambda}}$，即若$ |{s_n}(k)| \leqslant \delta $，则始终有$ |{s_n}(k + 1)| \leqslant \delta $成立.
推论1 记$ \omega = \max \left\{ {{{\left( {\dfrac{{{d^ * }}}{{{k_2}}}} \right)}^{\frac{1}{b}}},{\text{ }}{{\left( {\dfrac{{{k_2}h}}{{1 - {k_1}h}}} \right)}^{\frac{1}{{1 - b}}}}} \right\} $，则有
d∗h⩽k2hωb⩽(1−k1h)ω. (39) 证明 当$ \omega = {\left( {\dfrac{{{d^ * }}}{{{k_2}}}} \right)^{\frac{1}{b}}} \geqslant {\left( {\dfrac{{{k_2}h}}{{1 - {k_1}h}}} \right)^{\frac{1}{{1 - b}}}} $时，有$ (1 - {k_1}h){\omega ^{1 - b}} \geqslant {k_2}h $，故$ (1 - {k_1}h)\omega \geqslant {k_2}h{\omega ^b} = {d^ * }h $；
当$ \omega = {\left( {\dfrac{{{k_2}h}}{{1 - {k_1}h}}} \right)^{\frac{1}{{1 - b}}}} \geqslant {\left( {\dfrac{{{d^ * }}}{{{k_2}}}} \right)^{\frac{1}{b}}} $时，有$ {k_2}{\omega ^b} \geqslant {d^ * } $，$ (1 - {k_1}h){\omega ^{1 - b}} = {k_2}h $，故$ (1 - {k_1}h)\omega = {k_2}h{\omega ^b} \geqslant {d^ * }h $. 原式成立，推论1证毕.
设$ {s_n}(k) = \gamma \delta = \gamma \omega \psi (b) $，其中$ - 1 \leqslant \gamma \leqslant 1 $，故$ \delta > 0 $，由式（27）可得
sn(k+1)=sn(k)−k1hsn(k)−k2h[[sn(k)]]b+h˜dn(k)=(1−k1h)ωγψ(b)−k2h|ωγψ(b)|bsgn（γ）+h˜dn(k)=(1−k1h)ω|γ|ψ(b)sgn(γ)−k2hωb|γ|b×ψb(b)sgn(γ)+h˜dn(k). (40) 当$ |{s_n}(k)| < \delta $时，考虑以下2种情况：
情况1 当$ - 1 \leqslant \gamma < 0 $时，由式（40）得
sn(k+1)=−(1−k1h)ω|γ|ψ(b)+k2hωb|γ|bψb(b)+h˜dn(k). (41) 由于$ {\tilde d_n}(k) < {d^ * } $，所以
sn(k+1)<−(1−k1h)ω|γ|ψ(b)+k2hωb|γ|bψb(b)+d∗h. (42) 由推论1可知，${d^ * }h \leqslant {k_2}h{\omega ^b} \leqslant (1 - {k_1}h)\omega $，则有
sn(k+1)<−(1−k1h)ω|γ|ψ(b)+(1−k1h)ω|γ|bψb(b)+(1−k1h)ω=(1−k1h)ω(−|γ|ψ(b)+|γ|bψb(b)+1). (43) 因为$ 0 < {\text{|}}\gamma | \leqslant 1 $，由引理1可得
sn(k+1)<(1−k1h)ωψ(b)<ωψ(b)=δ. (44) 另外，由于$ {\tilde d_n}(k) > - {d^ * } $，由式（41）可得
sn(k+1)>−(1−k1h)ω|γ|ψ(b)+k2hωb|γ|bψb(b)−d∗h. (45) 由推论1可知$ - {d^ * }h \geqslant - {k_2}h{\omega ^b}$，则
sn(k+1)>−(1−k1h)ω|γ|ψ(b)+k2hωb|γ|bψb(b)−k2hωb=−(1−k1h)ω|γ|ψ(b)+k2hωb(|γ|bψb(b)−1). (46) 若$ \gamma \psi (b) \leqslant - 1 $，则$ {\text{|}}\gamma {{\text{|}}^b}{\psi ^b}(b) - {\text{1 > 0}} $，故
sn(k+1)>−(1−k1h)ω|γ|ψ(b)>−ωψ(b)=−δ. (47) 若$ - 1 < \gamma \psi (b) \leqslant 0 $，则$ {\text{1}} - {\text{|}}\gamma {{\text{|}}^b}{\psi ^b}(b){\text{ > 0}} $，由推论1可知$ - {k_2}h{\omega ^b} \geqslant - (1 - {k_1}h)\omega $，由式（46）可推出
sn(k+1)>−(1−k1h)ω|γ|ψ(b)−k2hωb(1−|γ|bψb(b))⩾−(1−k1h)ω|γ|ψ(b)−(1−k1h)ω(1−|γ|bψb(b))=−(1−k1h)ω(|γ|ψ(b)−|γ|bψb(b)+1). (48) 因为$ 0 < {\text{|}}\gamma | \leqslant 1 $，由引理1可得
sn(k+1)>−(1−k1h)ωψ(b)⩾−ωψ(b)=−δ. (49) 综上，当$ - 1 \leqslant \gamma < 0 $时，$ |{s_n}(k + 1)| \leqslant \delta $成立.
情况2 当$ 0 \leqslant \gamma < 1 $时，同理可证明有$ |s(k + 1)| \leqslant \delta $.
综上所述，当$ |{s_n}(k)| \leqslant \delta $时，$ |{s_n}(k + 1)| \leqslant \delta $恒成立，即一旦系统状态进入区域${\rm{\varLambda}}$，则之后任意时刻都将始终保持在区域${\rm{\varLambda}}$内.
根据步骤1和步骤2的证明结果，定理2得证.
3.2 跟踪误差界
引理2[10] 考虑有界序列${e_n}(k) $，若$\max |{e_n}(k - 1)| \leqslant \rho $，$\rho \in {\mathbb{R}^ + }$，则其分数阶微积分也有界，即
1ha−1L∑j=1(−1)j(a−1j)en(k−j)⩽(K−1)ρha−1，L∈Z+， (50) 式中：$K = \dfrac{{\varGamma (L + 2 - a)}}{{\varGamma (2 - a)\varGamma (L + 1)}} \geqslant 1$，$\varGamma (x) = \displaystyle\int_0^\infty {{\tau ^{x - 1}}e_n^{ - \tau }{\text{d}}\tau } $.
根据定理2，当滑模变量$ {s_n}(k) $进入界内后，有
sn(k)=Δen(k)+l1en(k)+l2Δa−1[[en(k)]]q⩽δ， (51) Δen(k)+l1en(k)+l2ha−1L∑j=0(−1)j(a−1j)[[en(k−j)]]q⩽δ, (52) 所以
Δen(k)+l1en(k)+l2hα−1[[en(k)]]q+l2ha−1L∑j=1(−1)j(a−1j)[[en(k−j)]]q⩽δ. (53) 由引理2可得
Δen(k)+l1en(k)+l2h1−a[[en(k)]]q⩽δ+l2(K−1)ρha−1， (54) 式中：$\rho \geqslant \max { [\kern-0.15em[ {e_n}(k - 1) ]\kern-0.15em] ^q} $.
由定理2可得，跟踪误差$ {e_n}(k) $最终有上界，即
|en(k)|⩽ψ(q)max{(ϑha−1l2)1q，(l2h2−a1−l1h)11−q}， (55) 式中：$ \vartheta = \delta + {l_2}\dfrac{{(K - 1)\rho }}{{{h^{a - 1}}}} $.
因此，当滑模变量进入区域${\rm{\varLambda}}$后，相应的跟踪误差将在跟踪误差界内运动. 当$ {\left( {\dfrac{{\vartheta {h^{a - 1}}}}{{{l_2}}}} \right)^{\frac{1}{q}}} \leqslant {\left( {\dfrac{{{l_2}{h^{2 - a}}}}{{1 - {l_1}h}}} \right)^{\frac{1}{{1 - q}}}} $，即$ {l_2} \geqslant \dfrac{{{\vartheta ^{1 - q}}{{(1 - {l_1}h)}^q}}}{{{h^{1 + q - a}}}} $时，跟踪误差界由后者决定. 由于$ 0 < h \ll 1{\text{ s}} $，所以
(l2h2−a1−l1h)11−q<(l21−l1h)11−q. (56) 此时，相较于整数阶滑模控制，分数阶滑模控制方法的跟踪误差界更小.
4. 磁悬浮转台控制系统验证
磁悬浮转台实验平台如图4所示，系统参数如表1所示. 动子由圆周型Halbach永磁阵列、背板以及测量目标物构成，定子由8个线圈、支承架和基座构成，6个位移传感器可实时测量动子的六自由度位移. 控制器采用NI的PXIe-8880和多功能IO模块PXIe-7856R. 首先，由PXIe-7856R对6路传感器信号进行A/D转换，经位移解算后数据被送入PXIe-8880控制器运算后输出控制信号，由PXIe-7856R模块产生模拟电压，通过PA12A进行功率放大并转换为驱动电流后送入8个线圈.
表 1 磁悬浮转台系统参数Table 1. Systematic parameters of maglev rotary table参数 数值 永磁体长（lm）、宽（wm）、高（hm）/mm 30、8、8 φm/（°） 7.5 r/mm 80 Br/T 1.2 线圈尺寸（lc、wc、hc、rc）/mm 60、10、10、10 线圈匝数 N/匝 300 m/kg 2.37 ${J_\alpha }$$({J_\beta })$、${J_\gamma }$/（kg·m2） 9.37×10−3 、1.87×10−2 采样间隔/ms 1 重力加速度/（m·s−2） 9.8 将所提方法与PID和整数阶离散时间滑模控制（DTSMC）的运动控制效果进行对比，取a = 0.5，N = 10. 其中，PID参数kp = 50.7，ki = 1，kd = 107.4，且各自由度的增益KPID如表2所示；DTSMC和所提方法参数中取b = 3/5，q = 3/5，各自由度的其余参数如表2所示.
表 2 控制器参数Table 2. Controller parameters自由度 PID DTSMC 和所提方法 KPID k1 k2 l1 l2 x,y,z 0.124 1 19.3 8.84 21.4 α,β 4.90 × 10−4 3.95 × 10−3 0.0763 0.0348 0.0846 γ 9.78 × 10−4 7.89 × 10−3 0.152 0.0694 0.169 各实验中动子的初始位移为（0, 0, 2 mm, 0, 0, 0）T. 设总采样点数为$ M $，跟踪误差的均方根和最大值分别为
{eRMS=√1MM−1∑k=0e2n(k), eMAX=maxk=0，⋯，M−1 en(k). 4.1 控制算法耗时
在1 s运行时间内，3种控制方法耗时的平均值分别为0.087 、0.220、0.244 s，满足采样周期1.000 ms的要求. 相对于DTSMC算法，分数阶微积分的引入使得所提方法的耗时增加了约10%，可通过改变式（21）中的L值来调节算法耗时.
4.2 阶跃信号跟踪
以1 s为时间间隔，在6个自由度上依次跟踪幅度为5.0 mm、5.0 mm、0.3 mm、0.87 mrad（0.05°）、0.87 mrad（0.05°）、5.24 mrad（0.3°）的阶跃信号，验证磁悬浮转台的阶跃响应性能， 跟踪效果如图5所示.
由图5可知，磁悬浮转台在所提出的分数阶滑模控制策略下，在保证调节时间短的条件下，超调量更小，综合来看，具有最佳的阶跃响应效果.
4.3 三角波跟踪
在x、y、γ自由度对频率为1 Hz，峰值分别为5 mm、5 mm、5.24 mrad（0.3°）的三角波进行跟踪，验证工作台的往复平移和旋转性能. 跟踪效果和跟踪误差ex、 ey和 ez的波形如图6所示，$ {e_{{\text{RMS}}}} $和$ {e_{{\text{MAX}}}} $如表3所示.
根据表3数据：磁悬浮转台在所提出的分数阶滑模控制策略下，对三角波的跟踪误差的均方根和最大值均最小； 与整数阶滑模控制方法相比，采用所提方法，x、y、γ自由度的跟踪误差均方根分别减小了12.8%、16.8%和23.7%，最大值分别减小了9.26%、13.00%和33.20%.
表 3 三角波跟踪误差Table 3. Tracking errors of triangular wave控制方法 x/mm y/mm γ/mrad eRMS eMAX eRMS eMAX eRMS eMAX PID 0.0148 0.0966 0.0155 0.0977 0.0706 0.212 DTSMC 0.0117 0.0875 0.0149 0.0953 0.0410 0.170 所提方法 0.0102 0.0794 0.0124 0.0829 0.0313 0.113 4.4 圆轨迹跟踪
令磁悬浮转台在x-y平面跟踪半径为5 mm的圆轨迹，跟踪效果如图7所示，对应跟踪误差的均方根$ {e_{{\text{RMS}}}} $和最大值$ {e_{{\text{MAX}}}} $如表4所示.
由表4数据可知：与整数阶滑模控制方法相比，采用所提方法，x、y自由度的跟踪误差均方值分别减小了6.39%和12.40%，最大值分别减小了9.90%和12.10%，磁悬浮转台在所提出的分数阶滑模控制策略下对圆轨迹的跟踪效果最优.
表 4 圆轨迹跟踪误差Table 4. Tracking errors of circular trajectorymm 控制方法 x y eRMS eMAX eRMS eMAX PID 0.0109 0.0686 0.0226 0.163 DTSMC 0.0101 0.0657 0.0251 0.158 所提方法 0.0091 0.0615 0.0220 0.138 5. 结 论
1） 为抑制系统非线性、强耦合性和不确定性对磁悬浮转台运动精度的影响，设计了基于扰动补偿的离散时间分数阶滑模控制策略以提高跟踪精度，并进行了实验验证.
2） 针对包含系统误差、六自由度间耦合项和外界干扰的集总扰动，设计了非线性干扰观测器进行扰动补偿，证明了估计误差有界且可调节到任意小.
3） 利用分数幂函数对抖振的抑制能力和分数阶微积分的全局记忆性设计运动控制器，证明了闭环系统稳定性和有限时间收敛特性.
4） 开展了磁悬浮转台系统六自由度轨迹跟踪实验，分别对阶跃信号、三角波和圆轨迹进行跟踪，系统在所提出的分数阶滑模控制策略下的跟踪误差最小，验证了所提控制策略的优越性.
-
图 6 磁悬浮轴承转子系统静态悬浮仿真响应
Figure 6. Static levitation simulation response of magnetic bearing-rotor system
图 7 磁悬浮轴承转子系统静态悬浮实验响应
Figure 7. Static levitation experiment response of magnetic bearing-rotor system
图 13 被控对象扫频试验原理框图
Figure 13. Sweep frequency experiment principle for controlled object
表 1 径向磁悬浮轴承结构参数
Table 1. Structural parameters of radial magnetic bearing
参数 取值 单个磁极匝数 N/匝 75 单个磁极面积 A/m2 4.05 × 10－4 偏置电流 I0/A 2 单边气隙 C0/mm 0.25 
下载: 导出CSV
表 2 理论和试验弯曲模态频率
Table 2. Theoretical and experimental bending mode frequencies
Hz 弯曲模态阶数 理论值 试验值 1 91.26 91.83 2 255.14 257.12 3 526.41 519.53 4 816.02 817.46 5 1046.77 1032.57
下载: 导出CSV
-
[1] SCHWEITZER G, MASLEN E, BLEULER H, et al. Magnetic bearing:theory, design and application to rotating machinery[M]. Berlin:Springer-Verilog, 2009:11-12. [2] YANNICK P,JEFFREY S,王小明. 应用API617第8版配备磁悬浮轴承的透平膨胀机-压缩机[J]. 风机技术,2016,58(5): 6-18. [3] HUTTERER M, SCHRÖDL M. Modeling and mu-synthesis control of a flexible rotor stabilized by active magnetic bearings including current free control[J]. Journal of Sound and Vibration, 2023, 546:117439.1-117439.18. [4] INOUE T, LIU J, ISHIDA Y, et al. Vibration control and unbalance estimation of a nonlinear rotor system using disturbance observer[J]. Journal of Vibration and Acoustics, 2009, 131(3): 031010.1-031010.8. [5] LI K X, PENG C, DENG Z Q, et al. Field dynamic balancing for active magnetic bearings supporting rigid rotor shaft based on extended state observer[J]. Mechanical Systems and Signal Processing, 2021, 158: 107801.1-107801.14. [6] TANG E Q, HAN B C, ZHANG Y. Optimum compensator design for the flexible rotor in magnetically suspended motor to pass the first bending critical speed[J]. IEEE Transactions on Industrial Electronics, 2016, 63(1): 343-354. doi: 10.1109/TIE.2015.2472534 [7] JIN C W, XIONG F, BAO Z P, et al. Research on modal vibration suppression of maglev steel strip based on phase compensator[J]. Journal of Sound and Vibration, 2018, 434: 78-91. doi: 10.1016/j.jsv.2018.07.043 [8] PENG C, ZHENG S Q, HUANG Z Y, et al. Complete synchronous vibration suppression for a variable-speed magnetically suspended flywheel using phase lead compensation[J]. IEEE Transactions on Industrial Electronics, 2018, 65(7): 5837-5846. doi: 10.1109/TIE.2017.2782204 [9] ZHENG S Q, HAN B C, WANG Y G, et al. Optimization of damping compensation for a flexible rotor system with active magnetic bearing considering gyroscopic effect[J]. IEEE/ASME Transactions on Mechatronics, 2015, 20(3): 1130-1137. doi: 10.1109/TMECH.2014.2344664 [10] ROY H K, DAS A S, DUTT J K. An efficient rotor suspension with active magnetic bearings having viscoelastic control law[J]. Mechanism and Machine Theory, 2016, 98: 48-63. doi: 10.1016/j.mechmachtheory.2015.11.012 [11] TANG E Q, FANG J C, HAN B C. Active vibration control of the flexible rotor in high energy density magnetically suspended motor with mode separation method[J]. Journal of Engineering for Gas Turbines and Power, 2015, 137(8): 082503.1-082503.10. [12] HAN X M, ZHOU J, ZHOU Y. Analysis and suppression of self-excited vibration of flexible rotor AMBs system[J]. Journal of Vibration Engineering & Technologies, 2021, 9(8): 1911-1922. [13] 魏彤,房建成. 基于双频Bode图设计磁悬浮弹性转子陷波器[J]. 光学精密工程,2008,16(5): 789-796.WEI Tong, FANG Jiancheng. Design of magnetically suspended elastic rotor Notch filter based on two-frequency Bode diagram[J]. Optics and Precision Engineering, 2008, 16(5): 789-796. [14] 崔培玲,赵光再,房建成,等. 基于相移陷波器的磁轴承不平衡振动全频自适应控制[J]. 振动与冲击,2015,34(20): 16-20,36.CUI Peiling, ZHAO Guangzai, FANG Jiancheng, et al. Adaptive control of unbalance vibration for magnetic bearings based on phase-shift Notch filter within the whole frequency range[J]. Journal of Vibration and Shock, 2015, 34(20): 16-20,36. [15] MUSHI S E, LIN Z, ALLAIRE P E. Stability analysis for a flexible rotor on active magnetic bearings subject to aerodynamic loads[C]//Proceedings of the 12th international symposium on magnetic bearings. Yokohama:[s.n.], 2010:22-25. [16] RAN S L, HU Y F, WU H C. Design, modeling, and robust control of the flexible rotor to pass the first bending critical speed with active magnetic bearing[J]. Advances in Mechanical Engineering, 2018, 10(2): 1-13. [17] 谭大力,陈进,廖明夫,等. 圆柱面配合激起的转子失稳振动研究[J]. 机械科学与技术,2014,33(12): 1786-1790.TAN Dali, CHEN Jin, LIAO Mingfu, et al. Instability caused by cylindrical surface fit in rotor system[J]. Mechanical Science and Technology for Aerospace Engineering, 2014, 33(12): 1786-1790. [18] 曾瑶,陈亮,蒋云帆,等. 圆柱面配合对柔性转子稳定性的影响[J]. 燃气涡轮试验与研究,2019,32(2): 42-48. doi: 10.3969/j.issn.1672-2620.2019.02.008ZENG Yao, CHEN Liang, JIANG Yunfan, et al. Analysis on the impact of cylindrical surface fit on flexible rotor[J]. Gas Turbine Experiment and Research, 2019, 32(2): 42-48. doi: 10.3969/j.issn.1672-2620.2019.02.008 [19] LIU Y, LIU H, YI J, et al. Investigation on the stability and bifurcation of a rod-fastening rotor bearing system[J]. Journal of Vibration and Control, 2015, 21(14): 2866-2880. doi: 10.1177/1077546313518817 [20] WU X L, JIAO Y H, CHEN Z B, et al. Establishment of a contact stiffness matrix and its effect on the dynamic behavior of rod-fastening rotor bearing system[J]. Archive of Applied Mechanics, 2021, 91(7): 3247-3271. doi: 10.1007/s00419-021-01963-9 [21] 魏彤,房建成. 磁悬浮控制力矩陀螺高速转子高频自激振动的抑制[J]. 宇航学报,2006,27(2): 291-296. doi: 10.3321/j.issn:1000-1328.2006.02.029WEI Tong, FANG Jiancheng. Self-excited vibration depression of high-speed rotor in magnetically suspended control moment gyroscope[J]. Journal of Astronautics, 2006, 27(2): 291-296. doi: 10.3321/j.issn:1000-1328.2006.02.029 [22] ZHOU Y, ZHOU J, MAHFOUD J, et al. Modelling and Validation of Rotor-Active Magnetic Bearing System Considering Interface Contact[C]//Proceedings of the 11th IFToMM International Conference on Rotordynamics. Mechanisms and Machine Science. Berlin:Springer, 2023:374-390. 期刊类型引用(0)
其他类型引用(1)
-
下载:
下载:
点击查看大图
计量
- 文章访问数: 311
- HTML全文浏览量: 201
- PDF下载量: 90
- 被引次数: 1