CN108534996A

CN108534996A - 一种基于微振动瀑布图的飞轮刚度参数提取方法

Info

Publication number: CN108534996A
Application number: CN201810142688.5A
Authority: CN
Inventors: 罗睿智; 张激扬; 李林峰; 翟百臣; 伏蓉; 冯洪伟; 杨震; 张阿妮; 高岩; 李建普
Original assignee: Beijing Institute of Control Engineering
Current assignee: Beijing Institute of Control Engineering
Priority date: 2018-02-11
Filing date: 2018-02-11
Publication date: 2018-09-14
Anticipated expiration: 2038-02-11
Also published as: CN108534996B

Abstract

一种基于微振动瀑布图的飞轮刚度参数提取方法，包括如下步骤：(1)建立飞轮转子的径向动力学模型和径向微振动传递模型；(2)对飞轮升速时的微振动进行测试并生成瀑布图；(3)得到飞轮微振动瀑布图上的局部峰值坐标；(4)通过对坐标数据进行分类和鲁棒回归，拟合得到飞轮在全转速范围内的模态及涡动曲线的函数表达式，将四个表达式移项相乘得到构造特征多项式；(5)得到优化的目标函数；(6)通过最优化的方法求解目标函数，得到飞轮径向动力学模型和振动传递模型中的各刚度参数。本发明可应用于设计飞轮时优化配置转子结构刚度，涡动及模态，也可作为飞轮批量生成过程中判定产品结构一致性与合格性的依据。

Description

一种基于微振动瀑布图的飞轮刚度参数提取方法

技术领域

本发明提出了一种飞轮刚度参数提取方法，属于制模型参数辨识技术领域。

背景技术

飞轮是航天器上应用最为广泛的惯性姿态执行机构，主要由轴承组件、电机组件、轮体组件和壳体组件构成。其中，轮体组件、轴承组件和电机组件中的旋转零部件共同构成飞轮转子。鉴于转子中存在众多的活动部件，加之各个零部件不可避免存在一定的加工制造误差和转子装配误差。一方面，使得转子存在一定的静动不平衡量，在高速旋转过程中将会产生一定的离心力和力偶；另一方面，轴承沟道、滚动体等的加工误差使得处于预紧状态的轴承在旋转过程中将会产生预紧力的波动。这些力和力矩作为振源共同作用于飞轮转子，当这些与转速相关的振源和转子的结构模态及涡动耦合时，可能会产生较大的共振。这些振动力和力矩再通过基座传递至星体，它们相对于发射段火箭导致的振动幅值小得多，因此对飞轮而言，其在运转过程中产生的振动通常称之为“微振动”；而对于航天器结构和其它载荷而言，这种振动通常称之为“扰振”。飞轮微振动是一系列频率离散且分布范围较宽的微幅复合振动。它具有微幅、固有、宽频、难控等一系列属性。

国内外研究表明，以飞轮为代表的高速旋转的惯性执行机构是卫星姿态稳定和精度保持的最大干扰源。这些振动力和力矩将对飞轮自身及星体产生一定的影响。随着现代科技的发展，尤其是航天科技的飞速发展，对卫星等航天器平台的姿态稳定性提出了越来越高的要求。飞轮转子在高速旋转的过程中产生的宽频微振动日益凸显，成为我国高分辨率卫星性能进一步提升的一个障碍。

经对现有技术文献的检索发现，目前对于飞轮的微振动特性研究的一个重要方向是研究飞轮的动力学模型，分析其动力学特征，模拟飞轮的微振动输出，但是这些模型参数没有试验基础，模型的准确性需要提高；另一方面，模型仅停留于理论分析，未对飞轮产品的科研生产提供积极的指导作用。

发明内容

本发明解决的技术问题是：本发明针对现有技术的不足，提供了一种基于微振动瀑布图的飞轮刚度参数提取方法，通过本发明的方法，结合微振动传递模型，从飞轮的微振动力和力矩瀑布图中提取了飞轮的刚度参数，刚度参数作为飞轮生产过程中的一个指标，判定飞轮转子的结构刚度的一致性与合格性。

本发明的技术解决方案是：一种基于微振动瀑布图的飞轮刚度参数提取方法，包括如下步骤：

(1)在惯性坐标系下，建立飞轮转子沿径向平动和绕径向转动的动力学方程为：

其中，矩阵G为:

其中，m_w表示飞轮转子的质量、I_wd表示飞轮转子径向转动惯量、I_wp表示飞轮转子轴向转动惯量，Ω表示转子转速，s为复变数；

x_w表示飞轮转子质心沿x方向的线位移、y_w表示飞轮转子质心沿y方向的线位移；

θ_w表示飞轮转子绕ox轴转动的欧拉角,表示飞轮转子绕oy轴转动的欧拉角；

f_wx(s)、f_wy(s)分别表示飞轮转子质心处的等效激励力在x轴、y轴的分量；

g_wx(s)、g_wy(s)分别表示飞轮转子质心处的等效激励力矩在x轴、y轴的分量；

K_wfx,K_wfy,K_wgθ,分别为飞轮定子对旋转体的支撑刚度；

K_wfθ,K_wgx,K_wgy分别为飞轮定子对旋转体的耦合刚度；

C_wx,C_wy分别为径向线振动的阻尼系数；

C_wθ,分别为径向角振动的阻尼系数；

建立飞轮转子产生的振动力和力矩的输出方程：

其中，矩阵P和矩阵C分别为:

其中，F_wx(s),F_wy(s)分别为飞轮微振动输出力在x轴、y轴的分量；G_wx(s),G_wy(s)分别为飞轮微振动输出力矩在x轴、y轴的分量；h表示参考坐标系和惯性坐标系的之间的距离；

(2)将步骤(1)建立的动力学方程代入输出方程得到飞轮内部激振力和力矩到基座上的振动力和力矩的传递函数阵T_H(s)为：

T_H(s)＝PCG^-1

振动传递函数阵T_H(s)分母为转子的特征多项式：E^t(s,Ω)＝|G|；忽略E^t(s,Ω)中的阻尼项，将s＝j2πf代入E^t(s,Ω)得到理论特征多项式：

其中，f表示频率，j为虚数单位；

(3)测得飞轮转子升速时的输出的微振动力和力矩，并生成微振动力和力矩瀑布图；

(4)在微振动力和力矩瀑布图中搜索飞轮微振动的局部峰值点，得到多组局部峰值坐标，局部峰值坐标由频率f和转速Ω组成；

(5)对局部峰值数据对进行分类和鲁棒回归，拟合得到飞轮在工作转速内的涡动及模态曲线频率其中，γ为第α条曲线的阶次；

对f_α(Ω)进行鲁棒回归得到第k次估计系数得到涡动及模态曲线的显函数，得到构造特征多项式：

将构造特征多项式E^c(f,Ω)展开并降幂排列得到：

其中，观测矩阵X_k＝[1,Ω,…,Ωⁿ]，Y_k为局部峰值坐标由频率f组成的向量，Λ_k为权矩阵；k为正整数；

(6)根据理论特征多项式E^t(f,Ω)和构造特征多项式E^c(f,Ω)的各阶次系数在任意转速下都有飞轮转子的质量m_w和转动惯量I_wd、I_wp，转子动力学径向刚度参数对等定律：K_wfθ＝K_wgy以及频率f，求得飞轮微振动传递函数阵T_H(s)中的刚度参数K_wfx,K_wfy,K_wgθ,K_wfθ,K_wgx,K_wgy的初值；

其中，v∈{0,2,4,6}；

(7)将经鲁棒回归后得到的涡动及模态曲线频率f_α(Ω)代入E^t(f,Ω)，并让特征多项式E^t[f_α(Ω)]加权平方和在全转速范围内积分最小，建立优化目标函数

求解得到刚度参数K_wfx,K_wfy,K_wgθ,K_wfθ,K_wgx,K_wgy；

其中，权函数ρ(Ω)为转速的指数函数：

Ω_u表示工作转速的上限；Ω_d表示工作转速的下限。

本发明与现有技术相比的优点在于：

(1)本发明从单端固支的飞轮转子的结构出发，建立了准确的飞轮转子径向微振动传递模型，相较于其他动力学模型，可更为准确地反映了转子的径向动力学特征，包括模态及涡动曲线分布等；

(2)本发明从飞轮转子缓慢升速过程的微振动瀑布图中提取结构参数，瀑布图中包含的信息反映了飞轮转子的结构特征，以及转子运转过程中表现出来的动态特征，因此微振动瀑布图中提取参数将更加全面更加准确；

(3)本发明采用鲁棒回归拟合了飞轮转子升速时表现出来的涡动曲线，再将拟合曲线代入到理论的特征多项式中，让特征多项式加权平方和在全转速范围内的积分最小，并利用优化方法得到了飞轮的动力学模型参数，也即结构参数，该参数的提取源于微振动试验数据，所得参数更为准确。

(4)本发明所提取的参数可作为一个判定依据，用于飞轮轴承的加载量判定，或者作为飞轮转子结构一致性的一个判定准则，可以直接指导生产。

附图说明

图1为本发明中飞轮的结构及坐标系定义示意图；

图2为本发明中飞轮基座上输出的径向微振动力矩瀑布图；

图3为本发明中模型参数辨识时所用的鲁棒回归流程图；

图4为本发明中飞轮基座上输出的径向微振动瀑布图等高线与用辨识的参数计算出的涡动曲线对比图。

具体实施方式

本发明方法首先建立了单端固定支撑飞轮转子的动力学模型，得到飞轮转子的微振动传递函数。其次，利用鲁棒回归拟合了微振动瀑布图中频率-转速涡动曲线，为4个显性的频率-转速分段函数；最后，将分段函数代入传递函数的理论特征多项式中，通过令其加权平方和在全转速范围内的积分最小，辨识出了最佳逼近模型参数。

本发明提供了一种基于微振动瀑布图的飞轮刚度参数提取方法，包括如下步骤：

(1)建立飞轮坐标系

如图1所示，惯性坐标系F_I：oxyz。其原点o位于支撑轴的根部中心处，以主轴根部平面的法线向上为z轴，y轴位于平面内且向右，x轴由右手法则确定。

本体坐标系F_b：o_bx_by_bz_b。固连于旋转体质心o_b，随旋转体同步旋转，其三个轴为旋转体的三个惯量主轴。

在零初始时刻，两坐标系F_b与F_I之间三轴对应平行，原点相距h，它们在飞轮结构体中的位置如图1所示。F_b可通过F_I依次进行如下旋转得到：F_I先绕ox轴转动θ_w；再绕旋转后的坐标系的oy′转动最后，绕上次旋转后的坐标系的oz″转动Ωt。其中，Ω表示转子转速，t表示转子运行时间。θ_w和表示了转子的径向角振动，是极小角度，因此它们的三角函数皆可近似线性化。

(2)建立飞轮转子的径向动力学模型

在惯性坐标系下，建立飞轮转子沿径向平动和绕径向转动的动力学方程为：

其中，矩阵G为:

其中，m_w表示转子的质量、I_wd表示转子径向转动惯量、I_wp表示转子轴向转动惯量，亦即转子的极惯量，s为复变数；

K_wfx,K_wfy,K_wgθ,为飞轮定子对旋转体在x和y方向的线刚度和角刚度；

K_wfθ,K_wgx,K_wgy为飞轮定子对旋转体在x和y方向的耦合刚度；

C_wx,C_wy分别为x和y方向的线振动的阻尼系数；

C_wθ,分别为x和y方向的角振动的阻尼系数。

飞轮转子产生的外传振动力和力矩：

其中，矩阵P和C分别为:

其中，F_wx(s),F_wy(s)：为飞轮基座输出的x和y方向的微振动力；

G_wx(s),G_wy(s)：为飞轮基座输出的x和y方向的微振动力矩。

(3)建立飞轮转子的径向微振动传递模型

将动力学方程(1)代入输出方程(2)得到了飞轮内部振源(力和力矩)到基座上输出的微振动(力和力矩)之间的传递函数矩阵为：

T_H(s)＝PCG^-1 (3)

微振动传递函数矩阵T_H(s)是刚度、阻尼、质量及转动惯量等结构参数和转子转速的函数，它是高速转子的结构本质的数学描述。它从频域内反映了飞轮转子径向振源与基座上径向微振动输出间的传递关系。

振动传递函数阵T_H(s)分母为转子的特征多项式：E^t(s,Ω)＝|G|。忽略E^t(s,Ω)中的阻尼项，并将s＝j2πf代入E^t(s,Ω)，可得到：

将(4)展开并按降幂排列成如下多项式形式：

(4)测试飞轮在缓慢升速过程中的微振动

测得飞轮转子缓慢升速时的微振动力和力矩，并生成瀑布图，如图2所示。

(5)搜索飞轮微振动瀑布图中的局部峰值点

从飞轮的微振动瀑布图中搜索飞轮微振动的局部峰值点，得到一系列的飞轮微振动瀑布图上的局部峰值坐标(频率f，转速Ω)，这些点主要属于飞轮转子的倍频激励和飞轮的模态及涡动耦合导致的局部共振点。

(6)鲁棒回归飞轮的涡动曲线

假设待估参数β_k为n维列向量，则

Y_k＝X_kβ_k+E_k (6)

其中，X_k是m_k×n观测矩阵X_k＝[1,Ω,…,Ωⁿ]；Y_k为由局部峰值坐标的频率f组成的向量，E_k为m_k维误差向量，k为正整数；则残差加权最小二乘为：

m_k为正整数；

其中，Λ_k为权矩阵，第k次估计系数：

若令(8)的权函数为单位阵Λ_k＝I，即为普通最小二乘回归，其回归系数：观测向量估计值则帽子矩阵：

其主对角元素即为观测值的杠杆值，由此可计算出第i个中心化杠杆值：

为矩阵H_k主对角线上的第i个分量；

残差中位数M_k：

其中，median为中位数函数。

误差项的标准离差估计s_k：

标准化残差：

其中，为上次回归的残差的第i个分量；c_t为调节常数。

鲁棒回归的权函数是上次回归残差的函数，对残差较大的数据赋予较小的权重，反之亦然。常用的权函数有Bisquare核函数、Cauchy核函数、Welsch核函数和Huber核函数等。其中Bisquare核函数为：

权矩阵：将式(14)代入(8)中即可得到该次估计系数观测向量估计值残差率为为：

若和上次回归的估计值之间的差异不超过参数估计的最大允许误差ξ，即则该次估计值即为最终的回归值；若差异较大，则应重新计算残差和权值，再次进行鲁棒回归。本发明暂且将阈值ξ定为0.05。

若残差率不超过残差率上限τ_u，即满足i∈{1,2,…m_k}，则将这些数据构成新数据对(X_k+1,Y_k+1)，重新回归，计算流程如图3所示，图中ξ为参数估计的最大允许误差。本发明暂且将阈值τ_u依次设定为0.1和0.03。图4为本发明中飞轮基座上输出的径向微振动瀑布图等高线与用辨识的参数计算出的涡动曲线对比图。

某型号50Nms飞轮的涡动曲线及模态曲线(频率f-转速函数Ω)的鲁棒回归结果f_α(Ω),α∈{1,2,3,4}如表1所示。显然，这四个多项式是转速Ω的显函数。

表1特征频率f-转速Ω拟合函数

名称	函数
		f₁(Ω)	87.62-1.37×10^-2Ω+7.75×10^-7Ω²
f₂(Ω)	87.31+1.31×10^-2Ω+1.49×10^-6Ω²-2.8×10^-10Ω³
		f₃(Ω)	183.42-7.1×10^-4Ω
f₄(Ω)	189.47+1.9×10^-3Ω-5.1×10^-9Ω²+2.29×10^-10Ω³

用该拟合曲线构造特征多项式如下：

将表1中的各个函数代入构造多项式E^c(ω,Ω)，展开并降幂排列成：

(7)模型参数辨识

理论模型和实际模型吻合，则有理论特征多项式E^t(f,Ω)和构造特征多项式E^c(f,Ω)几乎处处相等，即：

则，这两个特征多项式的各阶次系数在任意转速下都应近似相等，即：

表示理论特征多项式E^t(f,Ω)的各阶次系数，表示构造特征多项式E^c(f,Ω)的各阶次系数；

结合转子动力学径向耦合刚度参数对等定律：K_wfθ＝K_wgy；可得到上述某50Nms飞轮的频率如下：

根据飞轮转子的质量和转动惯量等结构参数，并结合频率计算公式：

可求得飞轮微振动传递模型中的刚度参数的初值如表2第3列所示。

表2某50Nms飞轮的结构参数

(8)模型参数优化

为了进一步减小拟合曲线和实际的涡动曲线之间的微弱误差。将鲁棒回归模型f_α(Ω),α∈{1,2,3,4}代入理论模型E^t(f,Ω)，并让特征多项式E^t[f_α(Ω)]加权平方和在全转速范围内积分最小，从而得到优化的目标函数；

其中，权函数ρ(Ω)不妨构造为转速的指数函数：

其中，Ω_u表示工作转速的上限；Ω_d表示工作转速的下限。

利用Matlab工具可求解上述最优问题，从而得到更为准确的模型刚度参数K_wfx,K_wfy,K_wgθ,K_wfθ,K_wgx,K_wgy，上述50Nms飞轮的刚度参数最优辨识结果如表2第4列所示。将这些刚度参数与飞轮刚度参数指标比对，判定上述某50Nms飞轮为合格产品，并且该飞轮转子结构和其他结构基本一致。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims

1.一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，包括如下步骤：

(1)在惯性坐标系下，建立飞轮转子沿径向平动和绕径向转动的动力学方程；建立飞轮转子产生的振动力和力矩的输出方程；

(2)将步骤(1)建立的动力学方程代入输出方程得到飞轮内部激振力和力矩到基座上的振动力和力矩的传递函数阵T_H(s)；

振动传递函数阵T_H(s)分母为转子的特征多项式E^t(s,Ω)，忽略E^t(s,Ω)中的阻尼项，将s＝j2πf代入E^t(s,Ω)得到理论特征多项式E^t(f,Ω)；

其中，f表示频率，j为虚数单位；

(5)对局部峰值数据对进行分类和鲁棒回归，拟合得到飞轮在工作转速内的涡动及模态曲线频率

对f_α(Ω)进行鲁棒回归得到第k次估计系数得到涡动及模态曲线的显函数，并得到构造特征多项式E^c(f,Ω)：

其中，γ为第α条曲线的阶次，α∈{1,2,3,4}；观测矩阵X_k＝[1,Ω,…,Ωⁿ]，Y_k为局部峰值坐标由频率f组成的向量，Λ_k为权矩阵；k为正整数；

其中，v∈{0,2,4,6}；

求解得到刚度参数K_wfx,K_wfy,K_wgθ,K_wfθ,K_wgx,K_wgy；

其中，权函数ρ(Ω)为转速的指数函数，Ω_u表示工作转速的上限；Ω_d表示工作转速的下限。

2.根据权利要求1所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，所述步骤(1)中，飞轮转子沿径向平动和绕径向转动的动力学方程为：

其中，矩阵G为:

K_wfx,K_wfy,K_wgθ,分别为飞轮定子对旋转体的支撑刚度；

K_wfθ,K_wgx,K_wgy分别为飞轮定子对旋转体的耦合刚度；

C_wx,C_wy分别为径向线振动的阻尼系数；

C_wθ,分别为径向角振动的阻尼系数。

3.根据权利要求1或2所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，所述步骤(1)中，飞轮转子产生的振动力和力矩的输出方程为：

其中，矩阵P和矩阵C分别为:

其中，F_wx(s),F_wy(s)分别为飞轮微振动输出力在x轴、y轴的分量；G_wx(s),G_wy(s)分别为飞轮微振动输出力矩在x轴、y轴的分量；h表示参考坐标系和惯性坐标系的之间的距离。

4.根据权利要求3所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，所述步骤(2)中，传递函数阵T_H(s)为：

T_H(s)＝PCG^-1；

转子的特征多项式E^t(s,Ω)＝|G|。

5.根据权利要求4所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，所述步骤(2)中，理论特征多项式E^t(f,Ω)为：

6.根据权利要求5所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，所述步骤(5)中，将得到的构造特征多项式展开并降幂排列得到构造特征多项式E^c(f,Ω)降幂排列后的表达式：

7.根据权利要求5或6所述的一种基于微振动瀑布图的飞轮刚度参数提取方法，其特征在于，

所述权函数