半球型動壓氣浮軸承陀螺儀的靜態誤差模型

李巖， 段富海

(大連理工大學 機械工程學院， 大連 116024)

平臺慣導系統的精度和成本瓶頸主要取決于陀螺儀[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.1 氣膜壓力場控制方程

以可壓縮氣體的Reynolds潤滑方程計算氣膜壓力場。計算中所涉及的氣膜最小厚度和最大厚度分別約為0.1和4.9 μm，取環境壓力為標準大氣壓，此時氣體分子的平均自由程λ=6.63×10-8m，克努森數Kn∈(0.014, 0.66)，需要考慮氣體稀薄效應。采用Fukui和Kaneko[19]提出的廣義氣體Reynolds潤滑方程作為該動壓氣浮軸承壓力場的控制方程，其量綱一化形式為

(1)

由于轉子位移u尺度為微米級，遠小于軸承半徑R，氣膜厚度可近似為

H=1+u·nξ/c+hg/cξ=1,2

式中：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，以將角度α控制在很小的范圍內，殼體與慣導平臺固聯。

對陀螺組件(浮筒、轉子和氣浮軸承)繞o軸運用動量矩定理，得到

(3)

Hrωi-Ki1+Mo=0

進而得到理想無干擾情況下的陀螺儀輸出i1=(Hr/K)ωi,以及干擾力矩引起的陀螺儀角速度誤差ωd=Mo/Hr。這里僅討論由動壓氣浮軸承氣膜變形引起的誤差，不計機械零件的加工誤差。

當環境存在比力f時，氣膜發生變形以產生對轉子的載荷F，即F=mf，m為轉子質量；同時，由于轉子產生位移u，浮液的浮心與陀螺組件的重心不再重合。而動壓氣浮軸承的剛度矩陣具有交叉項，因此轉子位移u和比力f不一定共線，如圖1所示，此時將產生力矩

M=-u×(mf)=-u×F

(4)

力矩M在o軸上的分量就是干擾力矩Mo，所以，動壓氣浮軸承的氣膜變形引起的陀螺儀角速度誤差為

ωd=(-u×F·eo)/Hr

(5)

式中：eo為沿o方向的單位向量。

2.2 回歸分析過程的簡化

圖4為各主要變量的計算關系，對給定的轉子位移u，通過潤滑計算求出載荷F，并由式(5)求出陀螺儀角速度誤差ωd，然后通過回歸分析得到陀螺儀角速度誤差ωd與比力f的關系，即靜態誤差模型。

由于給定參數是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的數據進行擬合，得到由于氣膜變形造成的陀螺儀角速度誤差

ωd=sgn(fz)Mρ(fρ,|fz|/Hr)sin(γf+ψ(fρ,|fz|))

(9)

3 計算結果與分析

根據第1節和第2節分析，使用MATLAB編寫計算程序，所使用的參數基于某型三浮陀螺儀的相關設計參數，如表1所示。

表1 計算參數

3.1 氣膜壓力分布

圖5為轉子位移u=(-0.7，0，-0.5)μm時，ξ=1的軸承在有溝槽區域和無溝槽區域的周向壓力分布對比。其中，φ1=0.8 rad處為無溝槽區域，φ2=0.5 rad處為有溝槽區域。可以看出，在有溝槽的區域存在較大的壓力波動。在高壓區，壓力波動更為劇烈，因為高壓區的氣膜厚度較小，即溝槽處和背脊處的氣膜厚度對比更加明顯。這種壓力波動是由溝槽引起的，并將潤滑氣體泵入較窄的通道，發生阻塞現象，從而提高了高壓區的壓力。壓力波動還會傳遞到無溝槽的區域，但會明顯減小。壓力波動現象與Feng[16]和Sahu[21]等得到的結果相吻合。

3.2 陀螺儀誤差模型回歸分析

圖6為干擾力矩的徑向分量Mρ和2個比力分量的關系。可以看出，當徑向比力fρ或軸向比力fz接近0時，Mρ也幾乎為0；只有當2個方向的比力均比較大時，才會產生較大的徑向干擾力矩，才有可能產生較大的陀螺誤差。如果只有軸向比力，則轉子只產生軸向位移，即兩者平行，因此其向量積為0，即干擾力矩為0。如果只有徑向比力，轉子只產生徑向位移，雖然兩者因偏位角θa的存在而不平行，但其向量積是軸向的，因此徑向干擾力矩Mρ仍為0，仍不會引起陀螺儀的測量誤差。另外，當fρ較小時，Mρ隨著fρ的增大而增大；當fρ較大時，Mρ反而隨著fρ的增大而減小。因為較大的徑向比力會引起較大的轉子偏心，隨著轉子偏心距的增大，氣膜高壓區會更靠近氣膜厚度最小點，即偏位角θa減小，導致位移與比力的向量積減小，即干擾力矩減小。上述偏位角隨偏心距增大而減小的結論已在本文計算的位移-載荷數據中得到驗證，限于篇幅不再展示結果。

由于計算采用參數方程，輸入參數是轉子位移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)給出插值得到的結果，通過觀察，構造回歸模型

b11fρ|fz|+b21fρ2|fz|+b40fρ4+b31fρ3|fz|

(10)

采用非線性最小二乘法擬合確定式(10)中的系數，得到結果如表2所示，擬合過程的均方根誤差為9.91×10-8N·m，擬合函數曲面如圖6(b)所示。

圖7為比力與干擾力矩的周向夾角ψ與2個比力分量fρ和fz的關系，其中圖7(a)為根據數據插值得到的曲面，圖7(b)為擬合函數曲面。可見，在不同的比力下，ψ幾乎不隨軸向比力變化，且整體上變化較小，在1.35～1.55 rad之間，即干擾力矩與比力相比超前1.35～1.55 rad。

用本節的方法對數據進行插值和擬合，構造回歸模型

(11)

擬合得到式(11)各系數的值如表3所示，擬合的均方根誤差為4.595×10-3rad。

將式(10)和式(11)代入式(9)，最終可得三自由度比力作用下由動壓氣浮軸承氣膜變形引起的陀螺儀靜態誤差模型

ωd=(b00+b10fρ+b01|fz|+b20fρ2+b11fρ|fz|+

b30fρ3+b21fρ2|fz|+b40fρ4+b31fρ3|fz|)·

sgn (fz)sin(γf+a00+a10fρ+a01|fz|)/Hr

(12)

表2 式(10)各系數擬合結果

系 數a00/rada10/(rad·s2·m-1)a01/(rad·s2·m-1)數 值1.482-3.746×10-23.608×10-3

對于任意給定的比力f=(fi,fo,fz)，將其化為柱坐標，并同表2和表3的各系數一起代入式(12)，即可得到陀螺儀的角速度誤差。

4 結 論

1) 通過引入2個中間參數——干擾力矩與比力的周向夾角和徑向干擾力矩，成功地將陀螺儀靜態誤差模型的三元回歸分析問題轉化為2個二元回歸分析問題。

2) 徑向干擾力矩隨著軸向比力的增大而增大，隨著徑向比力的增大呈現先增大后減小的趨勢。當比力的徑向或軸向分量為0時，轉子只有軸向或徑向的位移，不會產生陀螺儀靜態誤差。干擾力矩在周向上與比力相比較超前1.35～1.55 rad。

3) 建立的陀螺儀靜態誤差模型克服了潤滑計算中必須以轉子位移為自變量的不便，可直接預測300 m/s2以內的任意三自由度環境比力下由轉子位移引起的陀螺儀靜態誤差。