阅读说明：本技术 控制力矩陀螺输出力矩的预示方法 (Indication method for controlling moment gyro output moment ) 是由 韩勤锴 褚福磊 于 2020-06-15 设计创作，主要内容包括：本发明公开了一种控制力矩陀螺输出力矩的预示方法,包括：采用理论分析和数值仿真相结合的方式,考虑控制力矩陀螺CMG(control moment gyros)高速转子微振动的输出力矩模型,明确高速转子微振动源与输出力矩的耦合传递特性,定量揭示关键设计参数对整机输出力矩的影响规律,最后通过试验验证,理论计算干扰频率的误差小于5％,并为整机输出力矩干扰成分的优化提供理论依据。(The invention discloses a method for indicating the output torque of a control moment gyro, which comprises the following steps: by adopting a mode of combining theoretical analysis and numerical simulation, an output torque model of the micro-vibration of a control torque gyro CMG (control moment gyros) high-speed rotor is considered, the coupling transfer characteristic of a micro-vibration source and the output torque of the high-speed rotor is determined, the influence rule of key design parameters on the output torque of the whole machine is quantitatively disclosed, and finally, the error of theoretically calculated interference frequency is less than 5% through test verification, and a theoretical basis is provided for the optimization of the interference component of the output torque of the whole machine.)

控制力矩陀螺输出力矩的预示方法

技术领域

本发明涉及航天器结构动力学技术领域，特别涉及一种控制力矩陀螺输出力矩的预示方法。

背景技术

控制力矩陀螺(Control Moment Gyro,CMG)，是敏捷卫星等航天器实现快速姿态机动的理想执行机构。如何对CMG输出力矩进行优化，尽量降低CMG输出的干扰力和干扰力矩是高性能航天器亟待解决的关键问题。目前，国内外学者已开展了针对空间CMG输出力矩干扰因素的理论和实验研究工作，多数研究集中在分析高速转子系统的微振动机理方面，但仍存在一些问题尚未解决：(1)CMG高速转子微振动模型不准确；(2)缺乏考虑转子微振动的输出力矩模型；(3)转子微振动与输出力矩的耦合传递特性不清楚。

发明内容

本发明旨在至少在一定程度上解决相关技术中的技术问题之一。

为此，本发明的目的在于提出一种控制力矩陀螺输出力矩的预示方法，该方法有助于理清CMG扰振动力学特性，并为整机输出力矩干扰成分的优化提供理论依据。

本发明的目的在于提出一种控制力矩陀螺输出力矩的预示方法。

为达到上述目的，本发明实施例提出了控制力矩陀螺输出力矩的预示方法，包括以下步骤：步骤S1，将控制力矩陀螺CMG高速转子系统划分为五个子系统，包括高速转子、轴承组件、挠性支承、隔振装置和低速框架，并针对所述高速转子、所述轴承组件、所述挠性支承和所述隔振装置建立整机系统微振动分析模型，通过数值方法求解所述整机系统微振动分析模型，得到飞轮转子系统微振动响应；步骤S2，建立低速框架坐标系，计算所述高速转子自转和公转引起的绕低速框架三个坐标轴的陀螺力矩向量，处理所述陀螺力矩向量得到控制力矩陀螺CMG输出力矩的数学表达式；步骤S3，在所述飞轮转子系统微振动响应的基础上，预设所述低速框架的转速，并代入所述控制力矩陀螺CMG输出力矩的数学表达式，获得控制力矩陀螺CMG输出力矩的理论预测值。

本发明实施例的控制力矩陀螺输出力矩的预示方法，采用理论分析和数值仿真相结合的方式，先提出考虑CMG高速转子微振动的输出力矩模型，其次明确高速转子微振动源与输出力矩的耦合传递特性，定量揭示关键设计参数对整机输出力矩的影响规律，最后通过试验测试结果修正相关理论模型，并为整机输出力矩干扰成分的优化提供理论依据。

另外，根据本发明上述实施例的控制力矩陀螺输出力矩的预示方法还可以具有以下附加的技术特征：

进一步地，在本发明的一个实施例中，所述低速框架包括框架主轴、框架轴承和框架支撑结构，其中，所述框架主轴通过四点接触球滚动轴承支撑，通过所述挠性支承中的安装支架与所述高速转子连接，所述框架支撑结构与卫星本体连接，在电机的驱动下，所述框架主轴低速旋转，输出控制力矩，以控制卫星姿态。

进一步地，在本发明的一个实施例中，所述步骤S1进一步包括：

分别推导高速转子扰振方程组、轴承组件非线性支承力向量、挠性支承运动微分方程组和隔振装置运动微分方程组，通过整理、组集，得到整机系统微振动分析模型：

其中，

表示系统自由度向量，M,K₁,C₁,G分别为系统的质量、刚度、阻尼和陀螺矩阵，F_e,F_g表示不平衡质量激励和自身重力激励，F_b表示左右轴承引起的非线性扰振力向量。

进一步地，在本发明的一个实施例中，实际安装过程中，所述高速转子除考虑静不平衡引起的激励力f_u1，f_u2和动不平衡引起的激励力矩T_u1，T_u2外，还考虑旋转轴线与轴承中间线之间存在夹角所引起的不平衡扭矩T_s1，T_s2，计算公式为：

T_s1＝-(I_dfw-I_pfw)βΩ²cosΩt

T_s2＝(I_dfw-I_pfw)βΩ²sinΩt

其中，I_dfw为所述高速转子的直径转动惯量，I_pfw为所述高速转子的极转动惯量。

进一步地，在本发明的一个实施例中，所述陀螺力矩向量为：

G_y＝T_gfJ_fω_f

其中，G_y为陀螺力矩向量，J_f为高速转子的转动惯量向量，T_gf表示角速度矩阵，ω_f为瞬态角速度向量。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：

图1为本发明一个实施例的控制力矩陀螺输出力矩的预示方法的流程图；

图2为本发明一个实施例的小型控制力矩陀螺CMG结构示意图；

图3为本发明一个实施例的整机微振动分析模型示意图；

图4为本发明一个实施例的转子分析自由度示意图；

图5为本发明一个实施例的具有静、动不平衡质量的飞轮转子系统示意图；

图6为本发明一个实施例的角接触球轴承分析坐标系示意图；

图7为本发明一个实施例的轴承内外圈以及滚动体波纹度示意图；

图8为本发明一个实施例的受载前后滚道沟曲率中心及滚珠中心示意图；

图9为本发明一个实施例的受载状态下滚动体受力和力矩示意图；

图10为本发明一个实施例的轴承组件分析坐标系示意图；

图11为本发明一个实施例的整机微振动分析模型求解流程图；

图12为本发明一个实施例的高速转子角速度和框架坐标系示意图；

图13为一种CMG实验系统示意图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

下面参照附图描述根据本发明实施例提出的控制力矩陀螺输出力矩的预示方法。

图1是本发明一个实施例的控制力矩陀螺输出力矩的预示方法的流程图。

如图1所示，该控制力矩陀螺输出力矩的预示方法包括以下步骤：

在步骤S1中，将控制力矩陀螺CMG高速转子系统划分为五个子系统，包括高速转子、轴承组件、挠性支承、隔振装置和低速框架，并针对高速转子、轴承组件、挠性支承和隔振装置建立整机系统微振动分析模型，通过数值方法求解所述整机系统微振动分析模型，得到飞轮转子系统微振动响应。

具体地，如图2所示，小型CMG系统的结构主要包括高速转子、轴承组件、舱板壳体、安装支架、低速框架及其支承结构，其中，挠性支承包括舱板壳体和安装支架。一对高精度角接触球轴承面对面安装，内圈固定在主轴上，通过内外加载套筒施加一定的轴向预载。转子通过轴承安装壳与轴承外圈固定，工作状态下电机驱动绕主轴高速旋转。主轴一端与安装支架相连，而另一端支承在舱板壳体上。低速框架主轴四点接触球轴承支承，通过安装支架与高速转子连接。框架支承结构与卫星本体连接，在电机的驱动下框架主轴低速旋转，进而输出控制力矩，用于控制卫星姿态。CMG工作状态下，高速转子转速一般超过5000rpm，而低速框架转速最高不超过10rpm。二者转速相差悬殊，从振动的角度分析，二者的耦合较弱，低速框架对高速转子主要起支承作用。本发明实施例将低速框架视为刚体，而转子系统视为一个复杂的刚柔耦合系统，分析其振动响应对输出力矩的影响。

如图3所示，在振动分析建模时，分别针对转子、轴承组件、舱板壳体和安装支架四个子系统，建立整机非线性弹簧-质量集中参数分析模型，并分别针对上述子系统推导各自的扰振分析方程，最后通过方程组集，得到转子系统振动分析模型。对所建立的系统可以描述为：

(1)高速转子为复杂的轮辐式结构，忽略高速转子自身的结构振动，将转子视为集中质量，考虑其空间五个自由度运动，包括三个方向的平动和两个方向的摆动(忽略绕轴线的扭转振动)，建立分析坐标系X_fw-Y_fw-Z_fw。

(2)工作状态下的不平衡扰振力通过轴承支承力传递至主轴，进而通过安装支架、低速框架传递至卫星本体。在传递不平衡扰振力的同时，由于轴承内部Hertz接触、弹流润滑、波纹度等微观因素，也会产生内部激振力。在发明实施例通过建立角接触轴承五个方向的非线性扰振力(即除了绕轴线的扭矩外，三个方向的扰振力和两个方面扰振力矩)来描述上述两方面的扰振激励。对左、右轴承，扰振力分别为f_bLi,f_bRi，其中i＝1,2,3,4,5。

(3)分别采用集中质量m_s,m_f模拟固定主轴、安装支架的惯性。与转子类似，在X_fw-Y_fw-Z_fw坐标系下，集中质量m_s具有三个方向的平动和两个方向的摆动，与转子之间通过轴承非线性扰振力f_bLi,f_bRi进行力耦合。为了简化分析，集中质量m_f只考虑其沿三个方向的平动自由度，与集中质量m_s之间通过弹簧耦合，耦合刚度为k_x1,k_y1,k_z1。集中质量n_f通过挠性支架与卫星本体连接，挠性支架可简化为三个方向线性弹簧，其刚度用k_fx,k_fy,k_fz表示。此外，固定主轴的上端通常也支承在整机壳体上，两个方向的支承刚度为k_x2,k_y2。

首先，建立高速转子扰振方程组，如下：

(如图4所示)，高速转子为复杂的轮辐式结构，忽略转子自身的结构振动，将转子视为集中质量，建立分析坐标系X_fw-Y_fw-Z_fw。考虑转子空间五个自由度运动，包括沿三个方向的平动(u_fw,v_fw,w_fw)和绕X_fw,Y_fw轴的摆动

根据力的平衡关系，可得到高速转子扰振方程：

式中，m_fw,I_dfw,I_pfw分别为的质量、直径转动惯量和极转动惯量，f_bLi,f_bRi(i＝1,2,…,5)分别表示左右角接触轴承的非线性扰振力和力矩，具体求解将在轴承非线性扰振力中给出。f_u1,f_u2,T_u1,T_u2,T_s1,T_s2为不平衡质量、安装偏斜引起的激振力和力矩。

如图5所示，实际的轮辐和轮毂之间是通过加工装配而组合在一起的，因此不可避免的存在质量不平衡。本发明实施例中考虑两类不平衡激励：静不平衡和动不平衡。

对于静不平衡，引起时变力激励，沿X_fw,Y_fw轴的分量为f_u1,f_u2，具体表示为：

其中，U_s＝m_sr_ms表示静不平衡质量，

表示静不平衡量的初始相位角。对于动不平衡，将引起时变力矩激励。同样将其沿X_fw,Y_fw轴分解为T_u1,T_u2，具体为：

其中，U_d＝m_dr_{md f}表示动不平衡质量，表示动不平衡量的初始相位角。

在实际安装情况下，不可避免出现偏斜，即旋转轴线与轴承中心线之间存在夹角β。通过理论分析，认为偏斜同样会引起不平衡扭矩T_s1,T_s2，具体为：

T_s1＝-(I_dfw-I_pfw)βΩ²cosΩt

T_s2＝(I_dfw-I_pfw)βΩ²sinΩt (4)

接下来，获得轴承组件非线性支承力向量，即轴承非线性扰振力，本发明实施例通过载荷分析确定滚动体与内外滚道之间的接触变形和接触角，从而获得轴承非线性扰振力。对于角接触球轴承，在工作状态下通常需要施加预载，以保持滚动体与滚道之间的有效接触，进而提高轴承回转精度和支承刚度。基于Hertz接触理论和弹流润滑理论，先分析预载作用下的轴承载荷，再确定轴承非线性扰振力，具体如下：

(1)分析预载作用下的轴承载荷

如图6所示，外圈以Ω为转速匀速旋转，图6为轴承整体坐标系(X-Y-Z)和第j个滚动体局部坐标系(x_j-y_j-z_j)，其中，X-Y-Z中心固定于轴承轴心处，x_j-y_j-z_j则绕Z轴以ω_cj旋转(ω_cj表示滚动体公转速度)。不考虑轴承外圈绕Z轴的扭转，则外圈其余五个方向的自由度向量为：u＝[u,v,w,θ^u,θ^v]^T。相应的，五个方向所受的预载向量为：f_p＝[f_u,f_v,f_w,m_u,m_v]^T，其中，f_u,f_v,f_w分别为三个方向的力载荷，而m_u,m_v则为施加于X,Y轴上的转矩。

某时刻，第j个滚动体的位置角其中，φ₀初始位置角，N_b为滚动体个数。不失一般性，令φ₀＝0。该时刻，载荷使得外圈出现移动，则外圈位于第j个滚动体处的沟曲率中心沿y,z轴平移以及绕x轴的摆动位移向量分别为：δ_j＝[ξ_j,η_j,θ_j]^T。在小变形条件下，δ_j可由外圈位移u通过坐标变换得到：

δ_j＝T_ju+Δw_j (5)

式中T_j为变换矩阵，可表示为

其中，R_o＝0.5d_e-(r_go-0.5d_b)cosβ₀，表示轴承旋转轴心到外圈沟道曲率中心的距离。d_e表示轴承节圆直径，r_go表示外圈沟曲率半径，d_b表示滚珠直径，β₀表示轴承名义接触角。

公式(5)中Δw_j表示因轴承内、外圈以及滚动体波纹度引起的附加位移向量。如图7所示，p_ij,p_oj表示与滚动体j接触的内、外圈沿滚道周向的波纹度，而q_ij,q_oj则表示的是沿轴向的波纹度。w_ij,w_oj分别表示滚动体j与内滚道接触、外滚道接触时的波纹度。这些波纹度均为时间的周期函数，可分别采用谐波函数描述为：

p_ij＝∑A_incos(n_inω_cjt+2π(j-1)/N_b+ψ_in)

p_oj＝∑A_outcos(n_out(Ω-ω_cj)t+2π(j-1)/N_b+ψ_out)

式中，A_in,B_in,A_out,B_out,C_n，n_in,n_out,n_ba和ψ_in,

ψ_out,分别表示轴承各个部件波纹度的幅值、阶次和相位，ω_sj为滚动体自旋速度。纯滚动条件下，考虑到内圈固定、外圈旋转的结构形式，则滚动体公转速度可表示为：

自转速度为：

结合图6和7所示的几何关系可知：

如图8所示，不受载时，滚动体与内外滚道只是接触，但未出现变形，滚动体质心与内外滚道曲率中心共线，图8虚线所示，接触角为名义接触角β₀。外圈一旦承受载荷，首先使得外圈产生移动，由于滚动体与内外滚道接触变形，进而导致滚动体质心随之移动，三者重新达到平衡，图8实线所示。由于内圈固定，则内圈沟曲率中心始终保持不变，而外圈沟曲率中心的位移由式(5)可得为(ξ_j,η_j)。滚珠中心的位置为待定量，以(v_jy,v_jz)表示。由于载荷作用，使得内外接触角不等，且外滚道接触角减小为β_oj，而内滚道接触角增大为β_ij。

如图8所示，变形前，内圈沟曲率中心与滚珠中心的距离为L_ij，变形后该距离变为l_ij。对于外圈沟曲率中心与滚珠中心的距离，类似的，由变形前的L_oj变化为l_oj。若内外滚道沟曲率半径为r_gi,r_go，考虑到轴承各个部件的波纹度，则有：

L_ij＝r_gi+p_ij-h_ci-(0.5d_b+w_ij) (11)

L_oj＝r_go+p_oj-h_co-(0.5d_b+w_oj) (12)

其中，h_ci,h_co分别为滚动体与内外滚道之间的润滑油膜厚度。根据Hamrock和Dowson的研究结果，润滑油膜厚度可由如下经验公式计算得到

h_c＝2.69R_xU^0.67G^0.53W^-0.067(1-0.61e^-0.73κ) (13)

式中，表示无量纲速度，其中，η₀表示大气压力下、温度为20℃时润滑油粘度，

表示滚动体滚动方向的等效半径，外滚道接触取“+”，内滚道接触取“-”，

为等效旋转速度，且d_ro,d_ri分别为轴承外滚道和内滚道直径。G＝E′c_ηp表示无量纲弹性模量，其中c_ηp为黏压系数，E′＝E/(1-v²)表示有效弹性模量(E为材料弹性模量，v为泊松比)。

表示无量纲载荷，其中，Q_j为滚动体与滚道之间的接触载荷。

根据图8所示的几何关系，可得：

式(14-17)为载荷分析中所需满足的几何方程。除了几何条件外，还需满足载荷平衡条件，分别针对滚动体和轴承外圈。第j个滚动体的受力情况如图9所示。图中F_cj和M_gj表示滚动体旋转所受到的离心力和陀螺力矩，可知

和M_gj＝I_bω_sjω_cjsinα_j，其中，m_b,I_b分别为滚动体质量和转动惯量，α_j为自转轴与z轴夹角。纯滚动条件下，可得

滚动体与内外滚道的接触力分别由Q_ij,Q_oj表示。根据Hertz接触理论，可得：其中，K_i,K_o和δ_ij,δ_oj分别为滚动体与内外滚道的接触刚度系数和接触变形。由前面的几何分析可知，δ_ij＝l_ij-L_ij和δ_oj＝l_oj-L_oj。当δ_ij,δ_oj>0时，χ_ij＝1,χ_oj＝1；当δ_ij,δ_oj≤0时，χ_ij＝0,χ_oj＝0。对于接触刚度系数，由接触区的几何尺寸、材料弹性等因素共同确定，有其中

为当量主曲率半径，κ为椭圆率，ζ和∈分别为第一类和第二类椭圆积分，与k有关。根据Harris的专著即可分别计算对应于滚动体与内滚道、滚动体与外滚道的当量曲率半径和椭圆率，进而可得内外滚道接触刚度系数K_i,K_o。

如图9所示，在受力分析中，忽略平面外的摩擦力，假设陀螺力矩恰好被y-z平面内的摩擦力所产生的力矩所平衡。λ_ij,λ_oj表示滚动体与内外滚道接触的摩擦系数，有λ_ij,＝λ_oj＝1。根据图9所给出的滚动体受力关系，可得如下滚动体的力平衡方程：

式(18)和(19)为一对非线性代数方程。对于每一个滚动体，均存在类似的力平衡方程。

在载荷分析中，除了滚动体受力平衡外，还需考虑外圈的载荷平衡。第j个滚动体对外圈所施加的载荷大小为Q_oj和方向与图9所示相反。将这些载荷向外圈位于第j个滚动体处的沟曲率中心转化，有

则所有滚动体对外圈的作用合力：

其中，f_b＝[f_bu,f_bv,f_bw,m_bu,m_bv]^T表示所有滚动体对外圈作用合力向量，即轴承非线性支承力向量。对于轴承外圈，存在如下的受力平衡方程：

f_p-f_b＝0 (22)

联立方程(18)、(19)和(22)，采用Newton-Raphson方法迭代求解该非线性代数方程组。显然，该方程组的维数为：2N_b+5。进而，根据式(14-17)可得到滚动体与内外滚道之间的接触变形δ_ij,δ_oj和接触角β_ij,β_oj。同时也可得到给定预载下外圈的变形位移向量u_o。

(2)确定轴承非线性扰振力

如图10所示，X_fw-Y_fw-Z_fw表示轴承组件分析模型的坐标系，而X_bL-Y_bL-Z_bL和X_bR-Y_bR-Z_bR分别表示左右轴承分析坐标系。左右支承轴承距离质心的距离为l₁,l₂。由高速转子扰振方程组的分析可知，在X_fw-Y_fw-Z_fw坐标系下的位移向量为：主轴位移向量为：

由于左右轴承内圈与主轴固结，可通过主轴位移通过坐标换算得到轴承内圈位移。

根据图10可知，通过坐标变换并考虑预载引起的预变形u_o，可分别得到左右轴承外圈的相对变形向量为：

其中，

和分别表示左右轴承外圈的位移向量，而T_L和T_R表示相应的变换矩阵，根据图10中三个坐标系的定义可得：

将旋转体处的振动位移转换至轴承外圈的位移之后，通过前述中的非线性迭代方程精确求解轴承的非线性扰振力。

下面以左轴承为例对轴承非线性支承力的计算过程进行详细说明。

已知左轴承外圈位移向量

代入式(5)中可得左轴承外圈第j个滚动体处的扰振位移向量δ_jL。类似的，也可获得左轴承各个滚动体所满足的力平衡方程，如式(18)和(19)所示。联立式(18)和(19)，采用Newton-Raphson方法迭代求解该非线性代数方程组(维数为：2N_b)，求解得到各个滚动体在随体坐标系中的位移，进而，根据式(14-16)以及式(20)可得到第j个滚动体对轴承外圈所施加的载荷向量为：[Q_yjL,Q_zjL,Q_θxjL]^T。因此，在整体坐标系下，所有滚动体对外圈的作用合力(即左轴承对固定主轴的扰振力)可表示为：

同理，对于右轴承同样可得扰振力向量为：

由于Hertz接触非线性、润滑油膜以及滚动体与滚道之间的单向接触(仅存在压缩而不存在拉伸)，轴承扰振力与位移之间是非线性关系。

最后分别建立挠性支承运动微分方程组和隔振装置的运动微分方程组。

将固定主轴简化为集中质量m_s，具有五个方向的自由度，包括三个方向的平动u_s,v_s,w_s和两个方向的摆动在X_fw-Y_fw-Z_fw坐标系内，固定主轴与转子之间通过轴承非线性扰振力实现五个方向自由度运动的耦合。固定主轴一端通过耦合刚度k_x1,k_y1,k_z1与安装支架的质量连接，另一端通过舱板壳体来支承，其刚度为k_x2,k_y2。根据力的平衡关系，可得固定主轴的运动微分方程为：

其中，k_sxz,k_syz分别表示主轴横向与轴向运动的耦合刚度，I_sd表示主轴横向转动惯量。

对于隔振装置的集中质量m_f，除了与固定主轴通过刚度耦合外，自身通过挠性支架与卫星本体连接。挠性支架三个方向的连接刚度为k_fx,k_fy,k_fz。同样根据力的平衡关系，也可得到其运动微分方程组为

以上分别考虑固定主轴、隔振装置、舱板支架的柔性和质量惯性，建立各自的运动微分方程，其中固定主轴与转子之间通过滚动轴承施加非线性力耦合作用。

进一步地，对上述得到各个子系统的扰振分析方程进行整理、组集，可以得到整机系统微振动分析模型：

其中，

表示系统自由度向量，M,K₁,C₁,G分别为系统的质量、刚度、阻尼和陀螺矩阵，F_e,F_g表示不平衡质量激励和自身重力激励，F_b表示左右轴承引起的非线性扰振力向量。通过组集各个子系统的系数矩阵和外激励向量可分别得到上述矩阵和向量。具体表达式如下

M＝diag([m_fw,m_fw,m_fw,I_dfw,I_dfw,m_s,m_s,m_s,I_sd,I_sd,m_f,m_f,m_f]) (31)

G＝diag([ΩG_fw,0]) (34)

由于考虑转子与主轴集中质量之间通过非线性轴承扰振力耦合，因此刚度矩阵K₁中的子矩阵K_fw＝0，其余子矩阵表示如下：

K_f＝diag([(k_x1+k_fx) (k_y1+k_fy) (k_z1+k_fz)]) (38)

陀螺矩阵G的子矩阵可表示为：

由于考虑了滚动轴承非线性扰振力，该微分方程组为非线性耦合方程组，需通过数值方法求解。求解流程如图11所示，首先求解给定预载下外圈的预变形δ_oj和接触角β_oj；通过迭代求解非线性代数方程，分别得到左右轴承的非线性扰振力；代入整机系统微振动微分方程组公式(30)中，采用MATLAB中的微分方程求解器(如ODE45)计算当前时刻下的系统扰振响应。在每一个计算时间步，均需进行非线性迭代分析，以求解精确的轴承非线性扰振力。若计算时刻未超过设定的时长，则重复上述求解过程，直至计算时刻等于设定的计算时长，输出计算结果，即得到飞轮转子系统微振动响应。

在步骤S2中，建立低速框架坐标系，计算高速转子自转和公转引起的绕低速框架三个坐标轴的陀螺力矩向量，处理陀螺力矩向量得到控制力矩陀螺CMG输出力矩的数学表达式。

具体地，通过对整机系统微振动分析模型的分析求解，可得到系统五个自由度方向的振动位移和速度响应，即平动位移响应(u_fw,v_fw,w_fw)、平动速度响应和角位移响应

角速度响应

如图12所示，基于坐标变换理论，可以得到考虑转子振动影响后，瞬态角速度向量为

式中

表示绕w_fw的自转角速度，显然

表示坐标变换矩阵：

将式(41-42)代入式(40)中，经过整理可得

由上式可以看出，在不考虑高速转子振动时，即

瞬态角速度分量ω_fx＝0,ω_fy＝0,ω_fz＝Ω。一旦考虑振动引起的偏转运动，不仅导致绕Z_fw轴的角速度出现变化，而且在X_fw,Y_fw轴引起附加角速度分量。而这种现象将对CMG输出力矩产生显著影响。

进一步地，图12还给出了低速框架坐标系X_fg-Y_fg-Z_fg。电机驱动框架绕X_fg旋转，角速度为ω_cg。使得高速转子坐标系相对低速框架坐标公转，进而输出控制力矩。实际中，低速框架主轴难免存在安装偏斜，使得高速转子坐标系的公转速度不仅绕X_fg轴，而且在其他两坐标轴也可能存在分量。若假定安装偏斜在X_fg-Z_fg平面内，且与X_fg轴的夹角为β_gz，则高速转子的公转角速度分量为：

类似的，若安装偏斜在X_fg-Y_fg平面内，且与X_fg轴的夹角为β_gy，则高速转子的公转角速度分量为：

根据陀螺力矩的定义，因高速转子自转和公转引起的绕低速框架三个坐标轴的陀螺力矩向量可通过如下关系得到：

G_y＝T_gfJ_fω_f (46)

式中J_f＝diag([I_dfw,I_dfw,I_pfw])表示转子转动惯量向量，T_gf表示角速度矩阵

考虑到实验中输出力矩测试是固定方向的测试，而低速框架是不断旋转的，因此式(46)计算得到的陀螺力矩G_y需通过坐标变换，才能与实际测试的陀螺力矩G_yt进行对比，即为

因此，将式(43-45)代入式(46,48)中，整理后可得考虑高速转子振动的CMG输出力矩的表达式。

对于存在安装偏斜β_gz的情况：

类似的，对于存在安装偏斜β_gy的情况：

在步骤S3中，在飞轮转子系统微振动响应的基础上，预设低速框架的转速，并代入控制力矩陀螺CMG输出力矩的数学表达式，获得控制力矩陀螺CMG输出力矩的理论预测值。

也就是说，通过以上的分析、推导，可以得到CMG输出力矩的数学表达式，分别如式(49-54)所示。在得到转子整机振动响应的基础上，通过给定低速框架的转速，即可通过式(49-54)获得CMG输出力矩的理论预测值。

进一步地，本发明提出一种对CMG输出力矩模型动态试验测试以及与理论模型的预测值对比分析和验证的详细方法，包括：

(1)构建试验装置

如图13所示，振动试验系统包括：CMG通过三轴气浮台支承，放置在Kistler9253B型测力平台上；信号测试系统主要包括电荷放大器、OR35数据采集仪和测试计算机等；测试中以框架坐标系为定义坐标系，即X_fg为框架旋转轴线方向(也是垂直于测力台方向)，Y_fg,Z_fg分别为垂直于框架轴线方向，也即测力台的水平宽度方向和水平长度方向。测试中CMG转子系统将绕框架轴线旋转，初始时假定转子自转矢量方向与Y_fg重合，测力平台将实时测量绕框架三个坐标轴的动态力矩。

(2)模型参数分析方法

本发明实施例涉及的一种CMG高速系统的参数由表1列出，包括、支承结构以及不平衡激励参数等。本发明实施例给定内外滚道沿周向和轴向的波纹度幅值为2微米，阶数取为12阶(即为滚动体个数)，相位设置为零。由式(8)和(9)可知，内外滚道波纹度引起的特征频率是不同的。对于内滚道波纹度，其特征频率为：

对于外滚道波纹度，其特征频率为

表1一种CMG高速转子系统参数

(3)模型对比与验证方法

试验中需测试框架转速为一定值时，三轴输出力矩并绘制其相应的FFT频谱图。为了与实验结果进行对照，理论模型也需计算对应框架转速时，三轴输出力矩的时域波形图和相应的FFT结果，并对比实验测试和理论预测的高速转子振动频率值。

进一步地，本发明实施例根据建立的CMG输出力矩模型，确定影响CMG输出力矩的关键参数的方法，具体包括：

首先分析转子静不平衡量的影响。计算不同U_s值时，CMG的输出力矩及其频谱。进而分析动不平衡量的影响，分析不同U_d值时，CMG的输出力矩及其频谱。

分析轴承波纹度的影响，将轴承波纹度幅值增加到一定值，计算输出力矩的时域变化曲线和相应的FFT频谱。

分析安装偏斜角的影响，表征安装偏斜的指标是偏斜角β。根据式(4)可知，该角度的存在将引起附加的不平衡扭矩。计算β为一定值时，输出力矩的时域变化曲线和相应的FFT频谱。

分析CMG重力对输出力矩的影响，以模拟实际空间环境与地面模拟环境的差别。

分析转子轴承的预载的影响，以分析预载对输出力矩的作用机制。

进一步地，根据本发明实施例可以提出一种输出力矩干扰成分优化的方法，包括干扰成分溯源和干扰成分优化方法两部分。

明确高速转子系统设计参数、微振动和输出力矩之间的传递特性，对输出力矩高频干扰成分进行溯源，首先需要表征输出力矩干扰成分的大小。对于小型CMG，不考虑高速的振动，且认为安装良好、框架转速稳定时，式(49-51)可退化为理想状态下的输出力矩，可表示为：

Y轴和Z轴输出力矩较为相似，只是相差90度的相位。此时将Y轴输出力矩的干扰成分的均方根值来表征干扰成分的大小，即为：

式中，RMS(·)表示取干扰成分的均方根值。

根据式(60)所定义的干扰成分指标ΔG，计算不平衡量、轴承波纹度、安装误差(偏斜角)、考虑重力与否、轴承轴向预载以及支承挠性(主轴刚度、支架刚度)对输出力矩干扰成分的影响规律。或者针对不同干扰成分，结合材料、加工、装配等环节，给出抑制输出力矩干扰成分的优化策略。

再根据前述干扰成分溯源和影响分析结果，判断影响较为显著的因素。

上述的示例，经判断影响较为显著的因素为动不平衡量、轴承波纹度、安装误差、支承挠性，而静不平衡量、重力、轴承预载则几乎没有影响。因此，从设计角度，优化策略如下：

从输出力矩优化的角度，静不平衡量、轴承预载、重力参数可不必作为重点控制参数。从微振动角度，静不平衡量还是有很大影响。因此，从工程实践角度，应慎重考虑；

因动不平衡量、轴承波纹度、安装误差显著影响输出力矩的干扰成分，且为正相关，因此在实际设计、监控中，应严格控制三者的大小，尽量降低其对干扰成分的影响；

支承挠性对于输出力矩干扰成分的影响较为复杂，是优化分析的重要因素。在实际设计中，通过相应的优化设计，使支承挠性处于一个合理的范围内，设计出的CMG的输出力矩干扰成分最小化。

综上，本发明实施例提出的控制力矩陀螺输出力矩的预示方法，采用理论分析和数值仿真相结合的方式，先提出考虑CMG高速转子微振动的输出力矩模型，其次明确高速转子微振动源与输出力矩的耦合传递特性，定量揭示关键设计参数对整机输出力矩的影响规律，最后通过试验测试结果修正相关理论模型，并为整机输出力矩干扰成分的优化提供理论依据。

此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。