控制力矩陀螺(CMG)是空间站、高分辨卫星等航天器实现姿态调节和指向控制的重要执行机构[1]，具有工作效率高和输出力矩大的特点。CMG主要由高速转子、低速框架和轴承组成，低速框架通过改变高速转子的角动量方向，产生陀螺力矩，实现扭矩放大，进而完成航天器姿态调控[2]。高速转子的柔性特点和轴承磨损失效是CMG内部微扰动产生的源头，制约着航天器姿态高精度控制的进一步发展[3]。受到柔性转子和外框架耦合运动的影响，高速轴承往往处于时变速度和时变载荷的航天高工况，容易因滚动体异常滑动发生润滑失效和磨损，是引起CMG微振动和失效的重要因素。因此，建立CMG框架-轴承-柔性转子的耦合动力学模型，实现滚动体运动状态分析和测量验证，揭示滚动体空间运动规律，是找到轴承磨损抑制方法，提升控制力矩陀螺动力学特性的关键。
CMG动力学研究主要分为2方面。一部分学者将CMG当作力矩输出装置，利用相对简化的动力学模型进行奇异性分析[4]、抑振控制[5]、分配策略升级[6]、控制方法优化[7]、隔振器设计[8]等。上述研究中，CMG内部结构间的相互作用被弱化，CMG被简化为1~3个自由度并具备力矩输出功能的质点。另一部分学者主要通过逐步完善CMG动力学模型，分析CMG受到内部扰动和外部干扰时的动态响应[9-13]。Zhuo等[9]通过理论建模确定了多个谐波扰动的主要来源。Luo等[10]通过能量法提出较为系统的CMG动力学分析模型，考虑了轴承弹性及飞轮和框架的耦合运动，并进行了实验验证。Wang等[11]细化了轴承的作用，给出了具有非线性滚动轴承支承的CMG系统动力学建模思路。此外，还有学者考虑了转子的安装误差、不平衡问题、隔振器影响，建立了相应的动力学分析模型[12]。即使CMG中的许多因素得到了考虑，但现有模型主要将转子简化为质点，忽略了转子柔性的影响。工程实践表明：作为系统中主要柔性组件的高速轴承和转子的动态响应会改变CMG的内部扰动和输出特性，并影响航天器姿态调节的精度[13]。因此如何将转子柔性引入CMG动力学分析，构建框架自由度和转子自由度间的耦合关系，建立相应动力学模型是需要解决的问题。
相比CMG动力学的发展，轴承接触特性的分析方法已经有较为成熟的基础，经历了集中参数模型、准静态分析方法、准动态分析方法和动力学模型4个阶段[14]。这4种方法在适用范围、计算成本、运动描述等方面各具优点，被广泛用于故障建模、滚动体运动分析、保持架稳定性分析等[15]。针对CMG中的轴承，Han等[16]以高速轴承非线性力和力矩为中间量建立了分析模型，重点研究了不平衡力、接触表面波纹度和安装倾斜角等因素对CMG输出力矩、非线性支撑力和系统稳定度的影响。研究表明：轴承会受到系统动力学行为的影响，提供不同的支撑力，影响系统的稳定性和输出特性。然而，现有学者主要研究了高速轴承对CMG系统动力学行为的影响，针对CMG系统对高速轴承动态特性影响机理的研究还不深入。因此为了实现CMG中高速滚动体接触与运动特性的精准分析，需要将轴承分析方法与CMG柔性系统动力学实时耦合，建立CMG动力学-高速滚动体接触特性与运动学分析模型。
本文采用有限单元法和Lagrange法建立了考虑转子柔性、框架-轴承-转子耦合效应、滚动体接触行为的CMG动力学模型。基于CMG原理实验台和双目视觉测量方法对滚动体运动行为进行测量，验证了模型的准确性。进一步，分析了转子柔性对CMG中高速滚动体接触与空间运动特性的影响。
1 CMG柔性转子动力学建模
CMG的结构如图 1所示。
CMG动力学建模分为3部分：柔性转子动力学建模、轴承动力学建模和系统动力学建模。
其中：M是质量矩阵，C是阻尼矩阵，G是陀螺矩阵，K是刚度矩阵，x是广义位移向量，F是合外力向量。
为了便于将高速转盘的影响纳入整个系统，将其视为固连在转轴第j个节点位置，即转盘与第j个节点具有相同的运动行为。因此在系统方程第j个节点处考虑了转盘质量、陀螺效应及偏心力的影响。对于高速轴承，内滚道的位移等于当前有限单元节点的位移，外圈位移等于低速框架位移。
2 CMG轴承动力学建模
2.1 滚动体和保持架的相互作用
滚动体和保持架的相互作用见图 3，主要为接触压力和摩擦力。FpbjfXYp是OpXP和OpYP方向摩擦力的合力，Fbpjnp是接触压力。
滚动体j和保持架兜孔j之间的相对位置关系为
其中：$\overrightarrow{\boldsymbol{r}}_{\mathrm{b} j}^{\mathrm{I}}$ 是滚动体$j$ 在轴承惯性坐标系下的位置，$\overrightarrow{\boldsymbol{r}}_{\mathrm{c}}^{I}$ 是保持架在惯性坐标系下的位置，$\overrightarrow{\boldsymbol{r}}_{\mathrm{p} j}^{\mathrm{c}}$ 是保持架兜孔$j$ 在保持架本体坐标系下的位置，$\boldsymbol{T}_{\mathrm{Ic}}$ 和$\boldsymbol{T}_{\mathrm{pI}}$ 分别为保持架坐标系$O_{\mathrm{c}}-X_{\mathrm{c}} Y_{\mathrm{c}} Z_{\mathrm{c}}$ 到轴承惯性坐标系$O_{\mathrm{I}}-X_{\mathrm{I}} Y_{\mathrm{I}} Z_{\mathrm{I}}$ 和$O_{\mathrm{I}}-X_{\mathrm{I}} Y_{\mathrm{I}} Z_{\mathrm{I}}$ 到兜孔坐标系$O_{\mathrm{p}}-X_{\mathrm{p}} Y_{\mathrm{p}} Z_{\mathrm{p}}$ 的转换矩阵。式(2)中$\overrightarrow{\boldsymbol{r}}_{\mathrm{bp} j}^{\mathrm{p}}$ 在$X_{\mathrm{p}}$ 和$Y_{\mathrm{p}}$ 方向的坐标值为$\left(r_{\mathrm{bp} j X}^{\mathrm{p}}, r_{\mathrm{bp}, jY}^{\mathrm{p}}\right)$ ，令：
滚动体j与保持架兜孔j间的接触压力为
其中：Kl是润滑剂刚度，Kl=11/Cp；Kn是滚动体-兜孔接触变形系数；Cp是保持架兜孔游隙。滚动体j与兜孔j的碰撞力夹角为
保持架对滚动体j的接触压力在Xp和Yp方向的分量为
在分析滚动体受力时，需要将其转换到滚动体本体坐标系Ob-XbYbZb下。在保持架平动位移较小时，有(QpbjXb QpbjYb)=(QpbjXp QpbjYp)。
为确定滚动体j与兜孔j间摩擦力的方向，计算二者的滑动速度为
其中：$\overrightarrow{\boldsymbol{r}}_{\mathrm{b} j}^{\mathrm{I}}$ 和$\overrightarrow{\boldsymbol{r}}_{\mathrm{p} j}^{\mathrm{I}}$ 分别是滚动体j和兜孔j的平动速度，Obj是滚动体j质心，Cbj为滚动体j和兜孔j的接触点，TIb是滚动体本体坐标系到惯性坐标系的转换矩阵，TbI是惯性坐标系到滚动体本体坐标系的转换矩阵。滚动体j的自旋转速为
则滚动体j与兜孔j的相对滑动速度与Zpj的夹角为
因此在兜孔坐标系下，保持架对滚动体j的摩擦力及在各个方向的分量为
其中μc是滚动体和兜孔之间的摩擦系数。本文参考文[18]，将其设置为0.08。将摩擦力转化到滚动体坐标系中：
2.2 滚动体和滚道的相互作用
任一时刻下内滚道、保持架和滚动体j的广义自由度分别为
根据文[19]可以求解出滚动体j与内、外滚道的相对位置关系和接触力Qij、Qoj。此外，滚动体与内、外滚道之间存在牵引力和牵引力矩，与接触区域应力分布、润滑油特性、滚动体与滚道的滑动速度有关。
滚动体j与滚道接触界面的滑动速度由相对平动和滚动体自旋组成，其运动分析如图 4所示。下标i和o分别表示内和外滚道接触，滚动体相对于外滚道接触点(xo, yo)和内滚道接触点(xi, yi)的滑动半径分别为：
其中：Ro=foDb/(2fo-1)，Ri=fiDb/(2fi-1)，fo、fi分别为内、外滚道接触曲率半径系数，ai、ao分别为内、外滚道接触椭圆半长轴。Db为滚动体直径。滚动体与外滚动接触点的相对平动速度为
其中：Dm为轴承节圆直径，ωbjo为滚动体公转角速度，αo为滚动体-外滚道接触角。滚动体-外滚道接触自旋速度为
同理，滚动体与内滚道接触点的相对运动平动和自旋转速为
其中：ωi为轴承内圈转速，αi为滚动体-内滚道接触角。
由式(3)—(5)得到滚动体-滚道接触界面的滑动速度分布为
其中：$\boldsymbol{i}^{\prime}、\boldsymbol{j}^{\prime}$ 分别为$o^{\prime} x^{\prime}、o^{\prime} y^{\prime}$ 方向的单位向量。计算滚动体-外滚道接触状态和滚动体-内滚道接触状态时，$\Delta V^{x^{\prime}}、\Delta V^{y^{\prime}}、\omega_{\mathrm{s}}^{z^{\prime}}$ 分别取$\Delta V_{\mathrm{o}}^{x^{\prime}}、\Delta V_{\mathrm{o}}^{y^{\prime}}、\omega_{\mathrm{so}}^{z^{\prime}}$ 和$\Delta V_{\mathrm{i}}^{x^{\prime}}、\Delta V_{\mathrm{i}}^{y^{\prime}}、\omega_{\mathrm{si}}^{z^{\prime}}$ 。
接触椭圆坐标系下滚道对滚动体的牵引力和牵引力矩为
其中：$h$ 为中心油膜厚度，$\eta$ 为润滑油黏度分布，两者的求解方法参考文[20]，$a、b$ 分别为滚动体-滚道接触椭圆长半轴，短半轴，$\boldsymbol{k}^{\prime}$ 为接触平面$o^{\prime}-x^{\prime} y^{\prime}$ 法向$o^{\prime} z^{\prime}$ 的单位向量。由式(6)得到内滚道对滚珠的牵引力和牵引力矩分别为FXji、FYji和MZji，外滚道对滚珠的牵引力和牵引力矩分别为FXjo、FYjo和MZjo。
2.3 CMG轴承动力学方程
轴承动力学方程分为内滚道、滚动体、保持架和外滚道4部分。CMG高速轴承的内滚道与所在转轴节点固连，外滚道与低速框架固连，所以其动力学行为由系统动力学方程求解。
滚动体j除了受到来自滚道和保持架的作用力外，还存在由自身惯性产生的离心力及润滑剂的黏滞阻力：
其中：Cd是阻力系数，ρb和ρoil分别为滚动体和润滑剂密度。考虑滚动体j和兜孔j间摩擦力对滚动体转动的影响，得到滚动体j所受摩擦力矩为
其中：滚动体j的受力分布见图 5，最终动力学方程为
保持架与滚动体j之间的作用力在惯性坐标系下为
考虑保持架与引导面的相互作用力Fcoy、Fcoz，相互作用力矩Mcox[18]和滚动体对保持架兜孔的相互作用力，保持架的动力学方程为
其中: Nb是滚动体总数, mc是保持架质量，Jc是保持架转动惯量，yc和zc分别是保持架质心Oc沿OIYI和OIZI方向的平动位移，θc是保持架沿OIXI方向的转动角。
3 CMG系统动力学模型
3.1 轴承耦合方法
根据轴承的接触特性，高速轴承与转轴(内滚道)和低速框架(外滚道)的支撑刚度分别为[19]
其中：Nj是滚动体j的接触特性向量，与接触角和角位置有关；k0是预紧力下的线性化近似刚度。因此轴承支撑刚度矩阵为
转子-高速轴承系统与CMG中的转子坐标系不同，因此需要进行坐标转换，转换矩阵为
此外，CMG框架和基座间也存在轴承支撑，其刚度矩阵KG结构及取值参考文[21]。
3.2 系统动力学方程
CMG系统运动分析见图 6，其中Oin是CMG系统惯性坐标系原点，Og是框架本体坐标系原点。第i个节点的自由度定义为
框架为系统的节点(n+1)，框架位移为$\overrightarrow{\boldsymbol{r}}_{g}$ ，在惯性坐标系中，框架的运动速度包括平移和旋转：
在框架本体坐标系中研究柔性转子-高速轴承系统的运动学特性。转子节点i的位移为
其中：$\overrightarrow{\boldsymbol{r}}_{i 1}$ 是转子节点i在框架本体坐标下的初始位置，由转子系统质心位置决定；$\vec{r}_{i 2}$ 是转子节点i的位移矢量。
节点i的平动速度为
节点i的转动速度为：
框架的动能为
其中mg是框架质量。
转子节点i的动能为
其中mi是转子节点i的质量。
由此得到CMG系统的动能为
将CMG系统动能带入Lagrange方程中：
其中：f是系统广义坐标总个数。CMG系统动力学方程的质量矩阵和陀螺矩阵可通过Lagrange方程中与动能相关的2项来获取。每个柔性转子节点与低速框架之间的质量矩阵为
陀螺矩阵为
其中M11和G11与高速转子中每个节点的运动有关，是大小为5×5维的矩阵；M22和G22与框架运动有关，是大小为5×5维的矩阵；M33和G33是一个与框架旋转有关的数值；其他为耦合矩阵。本节构建系统动力学方程时，与柔性转子节点相关的M11和G11将被式(1)所导出的矩阵M和G取代。
高速轴承的刚度矩阵如3.1节所示，需要利用式(7) 进行坐标变换。
第1和第2章建立了考虑转子柔性的转子-轴承系统动力学方程。3.2节中利用Lagrange方程构建了低速框架与各有限元节点间的联系。
最终得出考虑转子柔性的CMG动力学方程为
其中：Mcmg是质量矩阵；Kcmg是刚度矩阵；Gcmg是陀螺矩阵；Ccmg=10-5Mcmg+10-6Kcmg；Fcmg是系统时变力向量；Xcmg是系统位移向量。系统刚度矩阵结构见图 7，有限元转轴刚度见式(1)，转子轴承刚度、框架轴承刚度见3.1节。其中，节点1、n、(n+1)为刚度耦合节点，以节点1为例，其刚度为有限元转轴刚度和转子前轴承刚度的叠加。
利用CMG动力学方程和轴承分析方法可实现CMG高速滚动体接触特性与运动学分析，算法流程如图 8所示。
4 模型验证
5 结果对比及分析
5.1 转子柔性对滚动体接触特性的影响
为分析转子柔性对CMG高速滚动体空间运动的影响，首先对滚动体接触力和接触角进行研究，工况设置如表 2所示。高速转子以恒定转速转动，低速框架转速以时间t为变量进行正弦运动，运动周期为T0=16 s，转速定义为
表 2 工况参数 |
| 参数类型 | 数值 | 参数类型 | 数值 |
| 转子极转动惯量 | 0.123 kg·m2 | 高速轴承型号 | B7005 |
| 转轴直径 | 32.0 mm | 总节点数 | 43 |
| 转轴长度 | 206.0 mm | 转盘节点序号 | 21 |
| 转子质量 | 7.5 kg | 转子转速 | 9 500 r/min |
| 轴承位置节点序号 | 1、41 | 框架极转动惯量 | 0.16 kg·m2 |
| 框架运动类型 | 正弦转动 | 框架径转动惯量 | 0.19 kg·m2 |
| 正弦周期 | 16 s | 框架质量 | 10.0 kg |
| 轴承预载 | 300 N |
此外，由图 13c可知，滚动体同时存在高频波动，波动的频率约为66.67 Hz，是轴承转速9 500 r/min对应频率158.73 Hz的0.42倍，这一频率对应轴承保持架转速的特征频率。这是因为转子在生产加工或装配过程中存在一定的静不平衡量和动不平衡量，对应离心力和力矩的变化周期由转子转速决定，轴承中滚动体的接触特性变化表现为保持架的特征频率。以滚动体和内滚道接触力为例对滚动体接触特性进行分析，值得注意的是，滚动体其他的动态行为特征如接触角、自旋转速、旋滚比等也表现出相同的变化性质。
进一步分析转子柔性对滚动体-滚道接触角的影响，结果如图 15所示。
滚动体与内、外滚道接触角在轴承初始接触角附近波动。为平衡离心力的作用，滚动体与内滚道接触角往往大于外滚道接触角。通过对比“柔性转子”模型和“刚性转子”模型可以看出，转子柔性会让CMG中高速轴承滚动体-滚道接触角变化更加极端，即内滚道接触角的最大值更大，外滚道接触角的最小值更小。这也和CMG柔性转子轴承需要承担更多的径向力有关。径向力的作用使轴承中“承载区”和“非承载区”的分布更加极端，所以接触角变化波动更大。
5.2 转子柔性对滚动体运动特性的影响
进一步分析转子柔性对CMG中高速滚动体运动特性的影响。计算了2种模型下滚动体俯仰角的变化趋势，如图 16所示。受到转子柔性的影响，滚动体俯仰角会更小，变化更加剧烈。滚动体俯仰角的变化趋势与滚动体-外滚道接触角十分相似。对于高速转子轴承，滚动体与外滚道的接触力更大，接触区域的润滑油牵引力和阻力矩也会更大，因此会呈现更小的相对滑动和自旋转速，在拟静力学轴承分析中往往将其近似为“外滚道控制”状态。本文中转子高速运动，采用轴承动力学的分析方法得到的结果也能用这一现象进行解释。在“外滚道控制”状态下，滚动体俯仰角与外滚道接触角直接相关，因此与外滚道接触角表现出相似的变化趋势，在考虑转子柔性时得到的滚动体俯仰角则会更小。
按照图 5坐标系的定义，对滚动体自旋转速进行分解，并选取其中较大的分量即“Xb方向自旋”和“Zb方向自旋”进行分析，结果见图 17。滚动体俯仰角由3个方向的自旋转速决定，但Yb方向自旋主要由滚动体陀螺力矩产生，转速远远小于Xb方向的自旋转速，因此在一定程度上可以认为滚动体俯仰角主要受到Xb方向自旋转速和Zb方向自旋转速的影响。与“刚性转子”模型相比，考虑转子柔性影响后的滚动体俯仰角更小，表明Xb方向的自旋转速更大，Zb方向自旋转速更小。由图 4可知，为平衡滚动体与外滚道间的相对运动，当滚动体与外滚道接触角更小时，需要滚动体表现出更大的Xb方向自旋转速和更小的Zb方向自旋转速。因此与刚性转子相比，在转子柔性作用下，滚动体Xb方向自旋转速更大，Zb方向自旋转速更小。
6 结论
为揭示CMG高速滚动体空间运动特性，本文建立了考虑转子柔性、框架-轴承-转子耦合效应、轴承运动行为的CMG动力学模型。通过搭建CMG高速轴承滚动体运动测量实验台，基于双目视觉方法对滚动体运动进行了测量，验证了理论模型的正确性。进一步讨论了转子柔性对滚动体运动的影响，主要结论如下：
1) CMG中高速轴承滚动体的动态响应同时存在低频波动和高频波动。低频波动主要受到框架转动所产生陀螺力矩的影响，变化周期为框架转动周期的一半。高频波动主要受到转子转速和轴承参数的影响，表现为保持架特征频率。
2) “柔性转子”模型考虑了系统径向转动和径向平动的耦合影响，将陀螺力矩更多转化为径向力，改变了轴承中“承载区”滚动体和“非承载区”的动态行为，使滚动体-滚道接触力和接触角的计算结果表现出更极端的变化。
3) CMG中转子的高速转动使滚动体的接触行为特征接近“外滚道控制”状态，因此滚动体俯仰角和外滚道接触角的变化趋势相似。考虑转子柔性时滚动体的俯仰角的计算结果会更小。
4) 与“刚性转子”模型相比，考虑转子柔性时计算得到的滚动体在Xb方向的自旋转速更大，Zb方向的自旋转速更小。随着滚动体-内滚道接触角增大，滚动体与内滚道接触点的自旋分量会变大，滚动分量会减小，因此滚动体与内滚道旋滚比的计算结果在考虑转子柔性时会更大。
