磁悬浮挠性转子系统模型复合辨识方法

张华强，王英广，赵学涛

（1. 山东理工大学 机械工程学院，淄博 255049；2. 北京控制工程研究所，北京 100094）

磁悬浮挠性转子系统模型复合辨识方法

张华强1，王英广2，赵学涛1

（1. 山东理工大学 机械工程学院，淄博 255049；2. 北京控制工程研究所，北京 100094）

针对磁悬浮挠性转子的纯有限元分析模型精度低的问题，提出一种磁悬浮挠性转子系统模型复合辨识方法。该方法首先采用有限元方法把磁悬浮挠性转子划分为多个Timoshenko梁单元，建立磁悬浮挠性转子系统模型，并对其挠性临界频率、阻尼、刚度和振型等关键参数进行分析，进而得到磁悬浮挠性转子的等效降阶数学模型；然后采用鲁棒自适应方法分析磁悬浮挠性转子系统的动态特性；最后，采用变LEVY方法从动态特性分析数据中对磁悬浮挠性转子进行系统辨识，校正有限元分析得到的降阶数学模型。实验结果表明，本方法可以得到磁悬浮挠性转子较为准确的系统模型。

挠性转子；有限元；鲁棒自适应；变LEVY方法；复合辨识

在磁悬浮离心机、涡轮机等超高速应用领域，转子一般设计成挠性的[1]。磁悬浮挠性转子系统是一个复杂的机电一体化系统，包括传感器、控制器、功率放大器、执行器(磁轴承)和挠性转子等组成部件[2-3]。获知磁悬浮转子系统模型是进行后续过临界转速控制的前提条件[4-5]。

磁悬浮转子刚性模型只是转子在低转速时的一种近似情况，是建立在转子自身无形变这一前提条件之下的，并只有六个自由度[6-7]。在进入挠性转速后，转子上任意两点的相对空间位置都是随时间变化的，其自身的形变不可忽略。挠性转子模型还受加工精度和装配条件影响较大，即便微小加工装配的误差，也会形成较大的模型差异[8]。

本文首先采用有限元方法把磁悬浮挠性转子划分为多个Timoshenko梁单元，建立磁悬浮挠性转子系统模型，并对其挠性临界频率、阻尼、刚度和振型等关键参数进行分析，得到转子系统的等效降阶数学模型；然后，使用鲁棒自适应方法对磁悬浮挠性转子系统进行动态特性分析；最后，采用变LEVY方法从动态特性分析数据中对磁悬浮挠性转子进行系统辨识，校正有限元分析得到的降阶数学模型，从而获得较为精确的磁悬浮挠性转子系统模型。

1 磁悬浮挠性转子系统建模

1.1 磁悬浮挠性转子有限元模型的建立及特性分析

采用有限元方法把磁悬浮挠性转子分割成一个个小单元，通过经典转子模型建模，再加上有限元之间的连接关系就组成单个有限元模型。通过一个个有限元模型的组成可以构建描述整个转子模型的矩阵，通过分析该矩阵即可获得磁悬浮挠性转子的特性参数。磁悬浮电机转子结构图如图1所示。

图1 磁悬浮电机转子结构图Fig.1 Rotor structure of magnetic levitated motor

单个有限元模型通常使用工程力学中的梁模型，目前的梁理论主要有Euler-Bernoulli梁和Timoshenko梁[9]。由于Timoshenko梁模型还涉及回转惯性和梁截面剪切变形，其表达式包含梁横截面回转所产生的动能和剪切变形所引起的弹性势能，因此采用Timoshenko梁模型可以得到很好的分析结果。

从磁悬浮挠性转子中划分 Timoshenko梁单元及其受力分析示意图如图2所示。

从图2中得到转轴挠曲q(z,t)动力学模型如下：

式中，E为弯曲弹性模量，I为有限元截面对中心轴的惯性矩，ρ为有限元材料密度，A为截面面积，Gshear为剪切弹性模量。式(1)中，第一项为弯曲变形势能项，第二项为径向运动动能项，第三项为主转动惯量动能项，第四项为主剪切变形势能项，第五项为合并的剪切变形和转动惯量耦合项。

图2 单个有限元划分及其受力分析Fig.2 Single finite element division and its stress analysis

图3 AR4模型描述单个有限元Fig.3 Single finite element described by AR4 model

式(2)后两项分别为磁轴承项和不平衡扰动项。根据Timoshenko梁理论的AR4模型，可把式(2)变换成下列表示形式：

传感器输出为

式中，S为观测矩阵。设式(3)特解的形式为

代入式(3)得其特征方程为

因陀螺耦合阵G不能对角化，从式(6)解得的特征值为复数形式，虚部对应转子模态频率。

对应的特征向量Φn1、Φn2(n=1,2,…,N)亦为复数形式，表现为空间三维曲线。

使用RotFE软件对4 kW磁悬浮电机转子进行有限元分析。将转子划分为55个有限元，前两阶挠性模态分析结果如图4、图5所示，一阶弯曲模态频率为691 Hz，二阶弯曲模态频率为1 442 Hz。如式(6)所示，随转速变化，转子的模态频率因陀螺效应出现分叉现象，升速时与前向模态频率(高频)相交，降速时与后向模态频率(低频)相交。使用RotFE软件绘制Compbell图如图6所示，升速临界转速43 000 r/min，降速临界转速在41 000 r/min。

图4 挠性转子一阶弹性模态空间运动轨迹Fig.4 Space motion trajectory of the first-order elastic modal

图5 挠性转子二阶弹性模态空间运动轨迹Fig.5 Space motion trajectory of the second-order elastic modal

图6 4 kW电机转子Campbell图Fig.6 Campbell diagram of 4 kW motor rotor

1.2 磁悬浮挠性转子系统模型的等效降阶处理

首先，将特征向量组成矩阵Φ：

高阶弹性模态频率一般远高于系统额定转速，而且因众多低通延时环节的存在，系统会对高阶模态失去控制能力。另外，真实系统中往往只有前两到三阶弹性模态对系统稳定性有影响，取Φ前r列组成rΦ对式(8)所描述系统进行降阶处理：

得到系统动力学方程：

观测方程为：

上述转换矩阵除Gr外均为对角形式，对式(10)进行拉氏变换，得到系统频域方程：

The Development of County Rural Tourism under the Background of the All-For-One Tourism——A Case Study of Fengning Manchu Autonomous County,Hebei Province____________________________XU Xinguo,WANG Wenxuan 7

从式(14)可见，转子特征频率处运动方程为：

图7所示为全阶系统和降阶系统不平衡相应比较，激励点旋转在最左端，响应点选择在转子中心，系统模型阶数从55降为8，降阶后只包含两个刚性模态和两个挠性模态。在额定转速的两倍转速120000r/min 以内，不平衡响应函数几乎完全重合，而且功放带宽有限，系统对高阶模态几乎不响应。因此，使用降阶模型已经足够精确，但计算量大大降低。

图7 全阶转子和降阶转子模型的不平衡响应Fig.7 Unbalance response of full-order rotor and reduced-order rotor model

系统模型写成状态方程的形式为：

2 磁悬浮挠性转子系统辨识方法

上文对磁悬浮挠性转子的模型进行了理论分析。挠性转子模型受加工精度和装配条件影响较大，最终导致基于线性—弹性模型的有限元分析精度较低，造成分析结果与实际结果有较大误差。因此，需要使用系统辨识方法对有限元模型进行校正。

系统辨识采用频域方法，并使用频率分析手段区分各个模态，对包含多模态的复杂系统有较好的效果。频率响应方法（FRF）使用扫频激励方式，可选择性地在各个频率给转子提供激励，给予各个模态以更多细节上的反应[10]。FRF模态辨识方法，由两部分组成，一是测得系统的频率响应，二是从频率响应中提出模态参数。下面将对其分别进行介绍。

2.1 磁悬浮挠性转子系统动态特性分析

如图8所示，磁悬浮转子频率响应主要包括激励和辨识两个环节。激励采用正弦扫频信号形式，频率变化可选择递乘形式。

在激励接入点之后的被测转子系统两端引出两个受激励的信号，送入自适应辨识环节。自适应辨识环节采用基于鲁棒自适应的方法对磁悬浮挠性转子系统的动态特性进行分析，其结构图如图9所示。

自适应辨识环节通过两个引出信号与激励信号间的幅值之比和相位之差，得到被测系统的动态特性。采用LMS自适应算法提取被辨识信号中的同频成分，拟合公式如下：

图8 磁悬浮转子频率响应分析结构图Fig.8 Analysis structure of magnetic flexible rotor frequency response

图9 鲁棒自适应辨识结构图Fig.9 Robust adaptive identification structure

LMS是一种递归收敛方法，式(17)的自适应步长采用固定步长形式。μ取值大，收敛速度快，但收敛误差大；μ取值小，收敛误差小，但收敛速度慢。而且不同频率端对收敛误差的要求不同。挠性模态附近频率响应函数变化快，要求收敛速度尽量快，但频率响应幅值大，对收敛精度敏感度低，要求μ取大。而在其他频率段，频率响应幅值小，辨识结果对收敛精度很敏感，要求μ取小。为平衡收敛速度和辨识误差间的关系，在此使用鲁棒变步长 LMS算法，以辨识误差均方差最小为原则，实时调整步长。调整公式为

式中，εmin为拟合误差的均方差。

在原始信号中叠加有大量噪声，影响辨识精度，需要对其进行预滤波处理。滤波器在二阶Butterworth陷波器进行改进，激励信号频率f和通带宽度控制参数sharp计算得到二阶带通滤波器系数(A,B)，计算公式为

其中，A=(A[0],A[1],A[2])为二阶带通滤波器分母多项式的系数，B=(B[0],B[1],B[2])为二阶带通滤波器分子多项式的系数。

得到被辨识模块输入输出信号在激励正余弦信号构造的直角坐标系中的投影坐标后，计算相频幅频变化公式如下:

得到控制电流到转子位移的幅频相频曲线如图10所示。为更好地对比鲁棒自适应和所设计单通滤波器的效果，图10将普通自适应、鲁棒自适应、带通滤波+鲁棒自适应进行对比。普通自适应方法噪声大，所产生的累积误差导致高频幅值偏大。鲁棒自适应在跟踪速度上虽然比普通自适应方法略差，但噪声明显减小，高频累积误差大幅降低。从图 10可见，加入带通滤波器后，噪声进一步降低。

图10 磁悬浮转子频率响应结果Fig.10 Frequency response of magnetic flexible rotor

2.2 基于变LEVY方法的系统辨识方法

得到转子系统频率特性后，获知转子模型还需从频率响应函数中提取转子的模态参数，得到转子系统模型，本文使用变LEVY方法实现这一过程[11]。转子系统传递函数如下：

上式简写作：

因[Z(s)]非奇异可逆，由式(24)可得：

式中：

式(28)略去了下标i j，式中，

为使误差方程线性化，将式(30)乘以Dk，称为加权误差：

Levy法矩阵方程变为

求解如上方程得到{a}、{b}之后即可解得转子模态参数。

3 系统模型复合辨识方法的效果及其分析

使用变LEVY方法对上述测试结果进行辨识。为提高辨识精度，采用分段拟合方法。首先对一阶模态频率附近区域进行辨识，辨识结果如图11所示。测试结果与辨识结果相差较大，主要是传感器与磁轴承轴向位置不同所致。对初始相位角修正-213°后，得到如图12所示的辨识结果，修正后的辨识结果与测试结果吻合度很好。一阶模态辨识结果如下：

同理，对二阶模态频率附近幅相曲线进行辨识得到如下结果，拟合曲线见图13。

图11 未修正初始相角时转子一阶挠性模态辨识效果图Fig.11 The first-order flexible modal identification effect without correcting the initial phase angle

图12 修正初始相角后转子一阶挠性模态辨识效果图Fig.12 The first-order flexible modal identification effect after correcting the initial phase angle

图13 磁悬浮转子二阶挠性模态辨识效果图Fig.13 The second-order flexible modal identification effect of magnetic flexible rotor

对低频段附近幅相曲线进行辨识得到如下结果，拟合曲线见图14。

将低频段，一阶模态、二阶模态进行综合得到如下传递函数：

图14 磁悬浮转子低频段辨识效果图Fig.14 Identification effect in low-frequency range of the magnetic flexible rotor

测试幅相特性曲线与辨识幅相特性曲线比较见图 15。本文目的是穿越磁悬浮挠性转子的临界转速，关心挠性模态频率附近的模型精度，而对其它频率段的辨识精度要求不高。从图15辨识结果看，挠性模态频率附近辨识模型与测试结果拟合程度很高。

图15 磁悬浮转子辨识结果图Fig.15 Identification effect of the magnetic flexible rotor

4 结 论

获知磁悬浮转子系统模型是进行后续过临界转速控制的前提条件。本文针对磁悬浮挠性转子的纯有限元分析模型精度低的问题，提出了一种磁悬浮挠性转子系统模型复合辨识方法。该方法首先采用有限元方法对磁悬浮挠性转子进行了多个 Timoshenko梁单元的划分，并建立磁悬浮挠性转子系统模型，对其挠性临界频率、阻尼、刚度和振型等关键参数进行分析，通过深入分析得到了磁悬浮挠性转子的等效降阶数学模型；然后采用鲁棒自适应方法分析了磁悬浮挠性转子系统的动态特性；最后，采用变LEVY方法从动态特性分析数据中对磁悬浮挠性转子进行了系统辨识，校正有限元分析得到的降阶数学模型。实验结果表明，该方法可以得到较为准确的磁悬浮挠性转子系统模型，具有很强的理论研究价值和工程意义。

（References）：

[1] Sun Zhen, He Ying, Zhao Jing-jing, et al. Identification of active magnetic bearing system with a flexible rotor[J]. Mechanical Systems and Signal Processing, 2014, 49(1-2): 302-316.

[2] 房建成, 张会娟, 刘虎. 磁悬浮刚性转子系统振动机理分析与动力学建模[J]. 控制理论与应用, 2014, 31(12): 1707-1713. Fang Jian-cheng, Zhang Hui-juany, Liu Hu. Vibration mechanism analysis and dynamic model development of magnetically suspended rigid rotor system[J]. Control Theory & Applications, 2014, 31(12): 1707-1713.

[3] Zhao Jie, Zhang Hai-tao, Fan Ming-can, et al. Control of a constrained flexible rotor on active magnetic bearings [J]. IFAC-Papers On Line, 2015, 48(28): 156-161.

[4] 黄梓嫄, 韩邦成, 周银锋. 非线性接触下磁悬浮电机柔性转子系统模态分析[J]. 中国电机工程学报, 2014, 34(15): 2438-2444. Huang Zi-yuan, Han Bang-cheng, Zhou Yin-feng. Modal analysis of the flexible rotor system of magnetic levitation motors under nonlinear contact[J]. Proceedings of the CSEE, 2014, 34(15): 2438-2444.

[5] Wang Y G, Fang J C, Zheng S Q. A field balancing technique based on virtual trial-weights method for a magnetically levitated flexible rotor[J]. Journal of Engineering for Gas Turbines and Power, 2014, 136(9): 1-7.

[6] 刘超, 刘刚, 盖玉欢. 基于磁力等效原理的刚性磁悬浮转子系统高精度在线动平衡[J]. 振动与冲击, 2016, 35(4): 67-71. Liu Chao, Liu Gang, Ge Yu-huan. Field balancing for a magnetically suspended rigid rotor based on magnetic forces equivalence principle[J]. Journal of vibration and shock, 2016, 35(4): 67-71.

[7] Tiwari R, Chougale A. Identification of bearing dynamic parameters and unbalance states in a flexible rotor system fully levitated on active magnetic bearings[J]. Mechatronics, 2014, 24(3): 274-286.

[8] Fang jian-cheng, Zheng Shi-qiang. AMB vibration control for structural resonance of double-gimbal control moment gyro with high-speed magnetically suspended rotor [J]. IEEE-ASME Transactions on Mechatronics, 2013, 18(1): 32-43.

[9] Shafiei N, Kazemi M, Ghadiri M. Comparison of modeling of the rotating tapered axially functionally graded Timoshenko and Euler-Bernoulli microbeams[J]. Physica E: Low-dimensional Systems and Nanostructures, 2016, (83): 74-87.

[10] 万文军, 李军. 一种基于LCR发散振荡响应的控制系统频率特性辨识方法[J]. 电机与控制学报, 2014, 18(11): 113-119. Wan Wen-jun, Li Jun. A method of frequency domain identification for control system based on process response in LCR diverging oscillation[J]. Electric Machines and Control, 2014, 18(11): 113-119.

[11] 朱萌, 曹国武, 张志伟, 等. 基于LEVY法的气动舵机系统辨识[J]. 弹箭与制导学报, 2011, 31(6): 69-72. Zhu Meng, Cao Guo-wu, Zhang Zhi-wei, et al. The system identification of pneumatic actuator based on levy method[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2011, 31(6): 69-72.

[12] 范达, 范春石, 贺杨, 等. 磁化悬浮-感应驱动型动量球系统[J]. 中国惯性技术学报, 2016, 24(4): 524-530. Fan Da, Fan Chun-shi, He Yang, et al. Momentum sphere system of magnetization suspension and induction driving [J]. Journal of Chinese Inertial Technology, 2016, 24(4): 524-530.

Compound identification method of magnetic flexible rotor model

ZHANG Hua-qiang1, WANG Ying-guang2, ZHAO Xue-tao1
(1. School of Mechanical Engineering, Shandong University of Technology, Zibo 255049, China;
2. Beijing Institute of Control Engineering, Beijing 100094, China)

Aiming at the problem of low accuracy of the pure finite element analysis model for the magnetic flexible rotor, a new compound identification method of magnetic flexible rotor model is proposed. First, a finite element method is used to divide the magnetic flexible rotor into a plurality of Timoshenko beam elements. Then, a magnetic flexible rotor system model is built, the key parameters of the flexible rotor such as critical frequency, damping, stiffness and vibration mode are analyzed, and then the equivalent reducedorder mathematical model is obtained. Based on this, the dynamic characteristics of the flexible rotor system are analyzed using a robust adaptive method. At last, a correction LEVY method is used to identify the magnetic flexible rotor system by analyzing the dynamic characteristics data of the rotor, and the equivalent reduced-order mathematical model obtained by the finite element method is calibrated. Experimental results show that this method can obtain a more accurate system model of the magnetic flexible rotor.

flexible rotor; finite element; robust adaptive; correction LEVY method; compound identification

V414

A

1005-6734(2017)02-0249-07

10.13695/j.cnki.12-1222/o3.2017.02.021

2017-01-011；

2017-03-24

山东省自然科学基金（ZR2015FL012）；国家自然科学基金（51605031）；山东省自然科学基金教育厅联合专项（ZR2014JL027）

张华强（1982—），男，博士，讲师，从事检测与导航技术研究。E-mail: huaqiang.zhang@163.com