平台惯导系统的精度和成本瓶颈主要取决于陀螺仪[1]。用于平台惯导的三浮陀螺仪是高精度陀螺仪，因采用动压气浮轴承来支承陀螺电机，实现转子在无接触状态下的高速旋转，在旋转速度、连续工作时间和能耗等方面均有很大优势，从而大幅延长使用寿命和提高精度[2-3]。在飞机、舰船、潜艇、战略导弹和载人航天等领域得到广泛应用[4]。
然而，动压气浮轴承的承载能力小、刚度低，这些缺点严重限制了三浮陀螺仪精度的进一步提高[5]。轴承气膜在比力作用下产生变形，导致整个陀螺组件的偏心，进而造成浮心与陀螺组件的重心不重合，产生干扰力矩和陀螺仪误差。目前有关陀螺仪误差模型的研究主要依靠实验[6-7]，很难单独考虑每个误差源对总误差的贡献。被广泛认可的单自由度陀螺仪静态误差模型[1, 8]是基于将陀螺组件视为线性弹性体提出的，而动压气浮轴承的刚度随位移发生显著变化，在平衡位置附近的刚度很小[9-10]。因此这种误差模型不能很好地描述带有动压气浮轴承的陀螺仪，有必要用润滑数值计算方法对其造成的误差进行单独研究。
动压气浮轴承的力学性能与其造成的误差有直接联系，已有许多学者对气浮轴承的力学性能展开了深入研究。这些研究大多数是针对径向轴承[11-13]和止推轴承[14-15]，而三浮陀螺仪中常使用结构更为复杂、可同时承受径向和轴向载荷的半球型轴承，与之相关的研究较少。其中，Feng等[16]研究了带有螺旋槽的半球型动压气浮轴承的力学性能，在一阶滑移条件下计算了其径向刚度和阻尼系数。Deheri和Patel[17]计算了球面与半球支座之间磁流体挤压膜的承载能力，但是并未考虑径向载荷。另外，可同时承受径向和轴向载荷的还有锥台形轴承和全球形轴承等，有一些与之相关的研究同时考虑了这2种载荷，例如，Cui等[18]自主设计了带有节流孔和狭缝的全球型气浮轴承，并给出其轴向与径向承载力和刚度矩阵的计算方法；陈广强等[4]用CFD数值仿真方法计算了锥台型动压气浮轴承的承载能力和刚度，并进行了实验验证。这些研究大多以求解轴承力学参数为目的，即使以陀螺仪中的气浮轴承为研究对象，其计算结果也均未涉及陀螺仪的误差模型。另外，润滑计算往往以转子位移或偏心率为输入条件，而在陀螺仪使用过程中，科研人员并不关心转子的位移，只关心环境对陀螺仪输出的影响。因此，有必要建立以环境比力为自变量，角速度误差为因变量的误差模型。
本文首先用有限差分法求解Reynolds润滑方程，计算带有螺旋槽的半球型动压气浮轴承轴向和径向的载荷-位移关系；然后使用位移和载荷数据计算因气膜变形造成的干扰力矩和陀螺仪静态误差；最后拟合了空间任意方向比力作用下的陀螺仪静态误差模型。
1 半球型动压气浮轴承的静态特性本文研究的半球型动压气浮轴承的结构如图 1所示。2个对置的半球组成的轴承固定于陀螺仪浮筒上，并支承着转子。轴承半径为R，宽度为b，轴承上刻有螺旋形沟槽，沟槽方向角为βg，2个对置半球的间距为d，轴承间隙c远小于其他结构尺寸。转子以角速度ω绕z轴旋转，在比力f的作用下，转子中心Or相对于轴承中心Ob的位移为u=(ux, uy, uz)，偏位角为θa。直角坐标系Obioz与轴承固联，i和o分别为陀螺仪的输入轴和输出轴，为便于研究，引入另一直角坐标系Obxyz，满足u始终在xz平面内。
|
| 图 1 半球型动压气浮轴承的结构 Fig. 1 Structure of gas-dynamic hemispherical bearing |
以可压缩气体的Reynolds润滑方程计算气膜压力场。计算中所涉及的气膜最小厚度和最大厚度分别约为0.1和4.9 μm，取环境压力为标准大气压，此时气体分子的平均自由程λ=6.63×10-8 m，克努森数Kn∈(0.014, 0.66)，需要考虑气体稀薄效应。采用Fukui和Kaneko[19]提出的广义气体Reynolds润滑方程作为该动压气浮轴承压力场的控制方程，其量纲一化形式为
|
(1) |
式中：θ和φ均为气膜在图 2坐标系中的坐标变量；P=p/Pa，p为气膜压力，Pa为环境压力；H=h/c为量纲一化参数，h为气膜厚度；值得注意的是，半球型轴承的润滑方程中，由于截面半径r=Rcos φ是变化的，轴承数Λ=6μω(Rcos φ)2/(Pac2)不再是常数，μ为气体黏度；Q=Qp/Qcon，Qp为Poiseuille流量系数，Qcon为连续Poiseuille流量系数，两者描述了气体稀薄效应的影响，Qcon=D/6，D为逆Knudsen数，
|
|
| 图 2 计算域的局部坐标系与网格划分 Fig. 2 Local coordinate system and meshing of calculation domain |
由于转子位移u尺度为微米级，远小于轴承半径R，气膜厚度可近似为
|
式中：hg为沟槽深度；nξ为轴承表面的单位法向量，以指向转子的方向为正，ξ=1表示z轴正半轴穿过的轴承，ξ=2表示z轴负半轴穿过的轴承，即nξ=(cos θcos φ, sin θcos φ, -(-1)ξsin φ)
1.2 控制方程的求解计算域所采用的局部坐标系和网格如图 2所示，容易得到Oθφ坐标系与Oxyz坐标系的变换关系为
|
(2) |
将计算域划分为180×31的网格进行后续计算。由于轴承间隙很小，在最小膜厚处，轴承数Λ约为3 000，此时如果用传统的方法，把含有Λ的剪切项作为差分的辅助项计算，容易导致迭代过程不稳定。因此，采用迎风格式求解Reynolds方程[20]，以保证得到收敛解。
假设气膜周向压力分布连续，且环境压力为标准大气压，两端面压力均为环境压力，可得到如下边界条件：
|
通过在计算域上的积分得到轴承载荷：
|
式中：Ωξ为轴承的2个半球型表面对应的计算域；dA为微元面积。
根据上述分析，可由任意的转子位移计算对应的轴承载荷，得到两者的变化关系，即动压气浮轴承的静态特性。这些载荷-位移数据将用于建立陀螺仪的静态误差模型。
2 静态误差模型的建立 2.1 陀螺仪的静态误差三浮陀螺仪的结构简图如图 3所示。转子在工作时高速旋转，并由气浮轴承和转子之间的气膜支承，气浮轴承与浮筒固联，浮筒和壳体之间的间隙内充满了浮液，使浮筒悬浮在壳体内，并用磁悬浮轴承限制除绕输出轴o以外的5个自由度的运动，浮筒绕o轴的转动角度α可使传感器产生感应电流i1，并经过放大电路反馈给力矩器，产生反馈控制力矩Mcmd，以将角度α控制在很小的范围内，壳体与惯导平台固联。
|
| 1—传感器；2—壳体；3—浮筒；4—转子；5—力矩器；6—气浮轴承。 图 3 三浮陀螺仪的结构简图 Fig. 3 Sketch of three-floated gyroscope structure |
对陀螺组件(浮筒、转子和气浮轴承)绕o轴运用动量矩定理，得到
|
(3) |
式中：Io
在稳态时有
进而得到理想无干扰情况下的陀螺仪输出i1=(Hr/K)ωi, 以及干扰力矩引起的陀螺仪角速度误差ωd=Mo/Hr。这里仅讨论由动压气浮轴承气膜变形引起的误差，不计机械零件的加工误差。
当环境存在比力f时，气膜发生变形以产生对转子的载荷F，即F=mf，m为转子质量；同时，由于转子产生位移u，浮液的浮心与陀螺组件的重心不再重合。而动压气浮轴承的刚度矩阵具有交叉项，因此转子位移u和比力f不一定共线，如图 1所示，此时将产生力矩
|
(4) |
力矩M在o轴上的分量就是干扰力矩Mo，所以，动压气浮轴承的气膜变形引起的陀螺仪角速度误差为
|
(5) |
式中：eo为沿o方向的单位向量。
2.2 回归分析过程的简化图 4为各主要变量的计算关系，对给定的转子位移u，通过润滑计算求出载荷F，并由式(5)求出陀螺仪角速度误差ωd，然后通过回归分析得到陀螺仪角速度误差ωd与比力f的关系，即静态误差模型。
|
| 图 4 主要变量的计算关系 Fig. 4 Calculation relationship among main variables |
由于给定参数是u，而拟合自变量是f，需要通过调整u的取值范围来控制f的取值范围，使其充满需要的定义域并尽可能地分布均匀。
比力f具有3个方向的分量，即需要拟合一个三元函数。由于数据量随着自变量数量呈指数型增加，拟合难度也随之大大提高，这里考虑利用对称性减少自变量。对轴承和转子而言，x轴和y轴是对称的，但对陀螺仪而言，i轴和o轴没有对称性。假设已知比力f引起的位移为u，干扰力矩为M，设另一情形的比力为f′，且f′等于f绕z轴转过了β角，即
|
(6) |
由轴承的回转对称性可知，f′引起的位移u′将满足：
|
(7) |
将式(4)、式(6)与式(7)联立并化简可得
|
即f′引起的干扰力矩M′也等于M绕z轴转过了β角。
为利用该对称性简化拟合过程，可在柱坐标系Obργz中研究该问题，并引入中间参数：干扰力矩与比力的周向夹角ψ=γM -γf和径向干扰力矩Mρ。其中，γM与γf分别为M和f在Ob ργz中的方位角。上述结论可表述为ψ=γM-γf与γf相互独立，所以不必考虑γf的变化，不妨规定转子位移的方位角γu=0。迭代计算时，取不同的uρ和uz，求出对应的f和M，通过拟合把ψ与Mρ表示成fρ和fz的函数ψ=ψ(fρ, fz)和Mρ=Mρ(fρ, fz)，则陀螺仪角速度误差为
|
(8) |
这样可将三元函数ωd(fi, fo, fz)的回归分析问题转化为2个二元函数ψ(fρ, fz)和Mρ(fρ, fz)的拟合。
另外，由对称性容易得到ωd(fi, fo, -fz)=-ωd(fi, fo, fz)，故可仅取fz＞0的数据进行拟合，得到由于气膜变形造成的陀螺仪角速度误差
|
(9) |
根据第1节和第2节分析，使用MATLAB编写计算程序，所使用的参数基于某型三浮陀螺仪的相关设计参数，如表 1所示。
| 参数 | 数值 |
| 轴承半径R/mm | 6 |
| 轴承宽度b/mm | 5 |
| 轴承间隙c/mm | 2 |
| 两轴承间距d/mm | 8 |
| 沟槽深度hg/μm | 1 |
| 沟槽数量Ng | 6 |
| 沟槽方向角βg/(°) | 45 |
| 转子质量m/g | 60 |
| 转子角动量Hr/(kg·m2·s-1) | 0.016 7 |
| 气体黏度μ/(Pa·s) | 1.79×10-5 |
| 转速n/(r·min-1) | 30 000 |
| 环境压力Pa/Pa | 1.013×105 |
3.1 气膜压力分布
图 5为转子位移u=(-0.7，0，-0.5)μm时，ξ=1的轴承在有沟槽区域和无沟槽区域的周向压力分布对比。其中，φ1=0.8 rad处为无沟槽区域，φ2=0.5 rad处为有沟槽区域。可以看出，在有沟槽的区域存在较大的压力波动。在高压区，压力波动更为剧烈，因为高压区的气膜厚度较小，即沟槽处和背脊处的气膜厚度对比更加明显。这种压力波动是由沟槽引起的，并将润滑气体泵入较窄的通道，发生阻塞现象，从而提高了高压区的压力。压力波动还会传递到无沟槽的区域，但会明显减小。压力波动现象与Feng[16]和Sahu[21]等得到的结果相吻合。
|
| 图 5 气膜压力分布 Fig. 5 Pressure distribution of gas film |
图 6为干扰力矩的径向分量Mρ和2个比力分量的关系。可以看出，当径向比力fρ或轴向比力fz接近0时，Mρ也几乎为0；只有当2个方向的比力均比较大时，才会产生较大的径向干扰力矩，才有可能产生较大的陀螺误差。如果只有轴向比力，则转子只产生轴向位移，即两者平行，因此其向量积为0，即干扰力矩为0。如果只有径向比力，转子只产生径向位移，虽然两者因偏位角θa的存在而不平行，但其向量积是轴向的，因此径向干扰力矩Mρ仍为0，仍不会引起陀螺仪的测量误差。另外，当fρ较小时，Mρ随着fρ的增大而增大；当fρ较大时，Mρ反而随着fρ的增大而减小。因为较大的径向比力会引起较大的转子偏心，随着转子偏心距的增大，气膜高压区会更靠近气膜厚度最小点，即偏位角θa减小，导致位移与比力的向量积减小，即干扰力矩减小。上述偏位角随偏心距增大而减小的结论已在本文计算的位移-载荷数据中得到验证，限于篇幅不再展示结果。
|
| 图 6 径向干扰力矩与2个比力分量的关系 Fig. 6 Radial interference torque versus specific force in two directions |
由于计算采用参数方程，输入参数是转子位移u而不是比力f，因此不能控制每个数据点的精确位置，只能通过调整输入位移间接地调整数据点的位置。为了提高回归分析的精度，先粗略调整转子位移的分布，使得到的数据点在Mρ(fρ, fz)的曲率较大区域分布较密集。为了避免数据点的疏密程度影响拟合误差的权重分配，先对数据点进行线性插值，得到均匀分布的插值点，再用这些插值点进行拟合。考虑到数据点围成的凸集以外的点不能作为插值点，先取一个较大的范围fρ＜400 m/s2且fz＜400 m/s2的数据点作为插值的源数据，而拟合时使用fρ＜300 m/s2且fz＜300 m/s2的插值结果，则得到的误差模型适用于比力|f|＜300 m/s2的情况。图 6 (a)给出插值得到的结果，通过观察，构造回归模型
|
(10) |
采用非线性最小二乘法拟合确定式(10)中的系数，得到结果如表 2所示，拟合过程的均方根误差为9.91×10-8 N·m，拟合函数曲面如图 6(b)所示。
| 系数 | 数值 |
| b00/(N·m) | 2.48×10-6 |
| b10/(N·s2) | -4.88×10-7 |
| b01/(N·s2) | 1.60×10-6 |
| b20/(N·s4·m-1) | -4.41×10-7 |
| b11/(N·s4·m-1) | -8.13×10-8 |
| b30/(N·s6·m-2) | 4.56×10-7 |
| b21/(N·s6·m-2) | -4.54×10-7 |
| b40/(N·s8·m-3) | -1.11×10-7 |
| b31/(N·s8·m-3) | 1.59×10-7 |
图 7为比力与干扰力矩的周向夹角ψ与2个比力分量fρ和fz的关系，其中图 7(a)为根据数据插值得到的曲面，图 7(b)为拟合函数曲面。可见，在不同的比力下，ψ几乎不随轴向比力变化，且整体上变化较小，在1.35~1.55 rad之间，即干扰力矩与比力相比超前1.35~1.55 rad。
|
| 图 7 比力与干扰力矩的周向夹角与2个比力分量的关系 Fig. 7 Circumferential angle between specific force and interference torque versus specific force in two directions |
用本节的方法对数据进行插值和拟合，构造回归模型
|
(11) |
拟合得到式(11)各系数的值如表 3所示，拟合的均方根误差为4.595×10-3 rad。
| 系数 | a00/rad | a10/(rad·s2·m-1) | a01/(rad·s2·m-1) |
| 数值 | 1.482 | -3.746×10-2 | 3.608×10-3 |
将式(10)和式(11)代入式(9)，最终可得三自由度比力作用下由动压气浮轴承气膜变形引起的陀螺仪静态误差模型
|
(12) |
对于任意给定的比力f=(fi, fo, fz)，将其化为柱坐标，并同表 2和表 3的各系数一起代入式(12)，即可得到陀螺仪的角速度误差。
4 结论1) 通过引入2个中间参数——干扰力矩与比力的周向夹角和径向干扰力矩，成功地将陀螺仪静态误差模型的三元回归分析问题转化为2个二元回归分析问题。
2) 径向干扰力矩随着轴向比力的增大而增大，随着径向比力的增大呈现先增大后减小的趋势。当比力的径向或轴向分量为0时，转子只有轴向或径向的位移，不会产生陀螺仪静态误差。干扰力矩在周向上与比力相比较超前1.35~1.55 rad。
3) 建立的陀螺仪静态误差模型克服了润滑计算中必须以转子位移为自变量的不便，可直接预测300 m/s2以内的任意三自由度环境比力下由转子位移引起的陀螺仪静态误差。
| [1] |
严恭敏, 李四海, 秦永元.
惯性仪器测试与数据分析[M]. 北京: 国防工业出版社, 2015: 10-55.
YAN G M, LI S H, QIN Y Y. Test and data analysis of inertial meter[M]. Beijing: National Defense Industry Press, 2015: 10-55. (in Chinese) |
| [2] | DELLACORTE C, RADIL K C, BRUCKNER R J, et al. Design, fabrication, and performance of open source generation Ⅰ and Ⅱ compliant hydrodynamic gas foil bearings[J]. Tribology Transactions, 2008, 51 (3): 254–264. DOI:10.1080/10402000701772579 |
| [3] |
秦冬黎. 一种球形气浮气动陀螺仪的设计方法及误差分析研究[D]. 哈尔滨: 哈尔滨工业大学, 2009: 11-27.
QIN D L. Research on design method and error analysis of a spherical gas-floated and driven gyroscope[D]. Harbin: Harbin Institute of Technology, 2009: 11-27(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10213-2010031070.htm |
| [4] |
陈广强, 杨云军, 雷娟棉, 等. 锥台型气体润滑动压轴承动力学数值模拟研究[J].
机械工程学报, 2016, 52 (4): 185–191.
CHEN G Q, YANG Y J, LEI J M, et al. Numerical simulation research on cone self-acting gas lubrication bearing dynamics[J]. Chinese Journal of Mechanical Engineering, 2016, 52 (4): 185–191. (in Chinese) |
| [5] | ARMENISE M N, CIMINELLI C, DELL'OLIO F, et al. Advances in gyroscope technologies[M]. New York: Springer Science & Business Media, 2010: 52-133. |
| [6] |
武丽花, 凌林本. 三浮陀螺仪漂移模型的建立及MATLAB实现[J].
中国惯性技术学报, 2004, 12 (6): 77–80.
WU L H, LING L B. Model of gyro drift and realizing in MATLAB[J]. Journal of Chinese Inertial Technology, 2004, 12 (6): 77–80. (in Chinese) |
| [7] | GOLIKOV A N, IGNATOVSKAYA A A. The statement of design and application questions for the gyroscope with a gas-dynamic suspension of ball rotor in the navigation support drilling system[J]. Journal of Physics, 2016, 671 : 012020. |
| [8] |
刘晶石. 气浮陀螺仪干扰力矩影响因素研究[D]. 哈尔滨: 哈尔滨工业大学, 2011: 18-31.
LIU J S. Research on influencing factors of interference torque of gas-floated gyroscope[D]. Harbin: Harbin Institute of Technology, 2011: 18-31(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10213-1012000318.htm |
| [9] | YANG Q, ZHANG H, LIU Y. Improved modified Reynolds equation for thin-film gas lubrication from an extended slip velocity boundary condition[J]. Microsystem Technologies, 2016, 22 (12): 2869–2875. DOI:10.1007/s00542-015-2667-4 |
| [10] | ZHANG W, MENG G, WEI X, et al. Slip flow and heat transfer in microbearings with fractal surface topographies[J]. International Journal of Heat and Mass Transfer, 2012, 55 (23-24): 7223–7233. DOI:10.1016/j.ijheatmasstransfer.2012.07.045 |
| [11] | CHENG F, JI W. A new model of water-gas turbulent lubrication for analysis of the static and dynamic characteristics in a journal bearing[J]. Proceedings of the Institution of Mechanical Engineers, Part J:Journal of Engineering Tribology, 2016, 230 (12): 1439–1451. DOI:10.1177/1350650116635927 |
| [12] | GERTZOS K P, NIKOLAKOPOULOS P G, PAPADOPOULOS C A. CFD analysis of journal bearing hydrodynamic lubrication by Bingham lubricant[J]. Tribology International, 2008, 41 (12): 1190–1204. DOI:10.1016/j.triboint.2008.03.002 |
| [13] | BHORE S P, DARPE A K. Investigations on characteristics of micro/meso scale gas foil journal bearings for 100-200 W class micro power systems using first order slip velocity boundary conditions and the effective viscosity model[J]. Microsystem Technologies, 2013, 19 (4): 509–523. DOI:10.1007/s00542-012-1639-1 |
| [14] | LIU R, WANG X L, ZHANG X Q. Effects of gas rarefaction on dynamic characteristics of micro spiral-grooved thrust bearing[J]. Journal of Tribology, 2012, 134 (2): 222011–222017. |
| [15] | GAD A M, KANEKO S. A new structural stiffness model for bump-type foil bearings:Application to generation Ⅱ gas lubricated foil thrust bearing[J]. Journal of Tribology, 2014, 136 : 0417014. |
| [16] | FENG K, LI W, XIE Y, et al. Theoretical analysis of the slip flow effect on gas-lubricated micro spherical spiral groove bearings for machinery gyroscope[J]. Microsystem Technologies, 2016, 22 (2): 387–399. DOI:10.1007/s00542-015-2487-6 |
| [17] | DEHERI G M, PATEL S J. Combined effect of slip velocity and surface roughness on a magnetic squeeze film for a sphere in a spherical seat[J]. Indian Journal of Materials Science, 2015, 1155 (10): 1–9. |
| [18] | CUI D P, YAO Y X, QIN D L. Study on the dynamic characteristics of a new type externally pressurized spherical gas bearing with slot-orifice double restrictors[J]. Tribology International, 2010, 43 (4): 822–830. DOI:10.1016/j.triboint.2009.11.009 |
| [19] | FUKUI S, KANEKO R. A database for interpolation of Poiseuille flow-rates for high Knudsen number lubrication problems[J]. Journal of Tribology-Transactions of the ASME, 1990, 112 (1): 78–83. DOI:10.1115/1.2920234 |
| [20] |
黄平.
润滑数值计算方法[M]. 北京: 高等教育出版社, 2012: 93-103.
HUANG P. Lubrication numerical calculation methods[M]. Beijing: Higher Education Press, 2012: 93-103. (in Chinese) |
| [21] | SAHU M, SARANGI M, MAJUMDAR B C. Thermo-hydrodynamic analysis of herringbone grooved journal bearings[J]. Tribology International, 2006, 39 (11): 1395–1404. DOI:10.1016/j.triboint.2005.11.022 |