液浮陀螺儀的流固熱仿真及特性分析

古 旻, 張 晶, 樊景松, 蘇 巖

（1. 南京理工大學機械工程學院, 南京 210094;2. 北京航天控制儀器研究所, 北京 100039）

0 引言

液浮陀螺儀具有體積小、 精度高、 動態(tài)性能好、 過載能力強等特點, 廣泛應用于船舶導航、導彈的精確打擊以及航天航空領域［1-2］。 在實際工況下, 由于內部熱源發(fā)熱和外部加熱片的保溫加熱等以熱傳導和對流換熱方式向其他部件散發(fā)熱量, 使液浮陀螺儀內部存在著溫度梯度, 陀螺儀內部溫度不同導致浮油密度差異, 進而引起浮液流動, 影響了液浮陀螺儀的精度［3］。 隨著液浮陀螺儀精度要求的不斷提高, 儀表內部溫度分布對精度的影響愈發(fā)明顯, 研究人員越來越重視液浮陀螺儀的溫度特性。 在掌握儀表溫度特性的基礎上,對結構進行合理的改進, 并采用精準的溫控系統(tǒng),便能提高陀螺儀溫度分布的均勻性, 減小由溫度因素帶來的干擾力矩, 提高儀表的精度［4-5］。 就目前而言, 相關技術已經成為制約我國慣導精度提高的關鍵問題之一, 而有限元分析是研究此類問題的關鍵工具［6-7］。

針對以上問題, 國內外學者進行了大量研究。楊盛林等［8］利用有限元軟件ANSYS 得到了陀螺儀殼體的溫度分布。 白永杰等［9］運用有限元方法對液浮陀螺儀浮子進行了熱應力分析, 得到了工作環(huán)境下浮子結構的最大熱應力值, 為陀螺儀設計時材料的選擇提供了重要依據(jù), 但其未考慮內部流體對浮子的應力情況。 卜石等［10］使用ANSYS Fluent 有限元分析軟件得到陀螺儀整體的溫度梯度,但將浮液視為固體, 并未考慮其流場對溫度分布的影響。 牛文韜等［11］基于ANSYS/CFX 對液浮陀螺儀進行流固熱耦合仿真, 分別計算了浮筒外側的儀表部件-浮油耦合模型以及浮子-內部氣體耦合模型, 得到了陀螺儀的溫度場以及內部流體的流場。

本文使用多物理場有限元仿真軟件對陀螺儀整表進行計算, 分別考慮了其內部流體及固體的特點, 首先介紹了液浮陀螺儀的工作原理, 其次分析了流體及傳熱問題有限元仿真計算的理論基礎, 然后建立了液浮陀螺儀流固熱多物理場耦合的仿真模型, 根據(jù)實際工作情況分析或理論計算得到模型的邊界條件, 最終計算并分析了內部流體的流場以及儀表整體的溫度場特性, 為實現(xiàn)液浮陀螺儀的分區(qū)溫控以及儀表改進提供依據(jù)。

1 液浮陀螺儀工作原理

液浮陀螺儀利用高速旋轉剛體具有的定軸性和進動性測量角速度, 本文所研究的液浮陀螺儀馬達為動壓馬達, 采用氣浮軸承支持, 馬達轉子高速旋轉運動產生較大角動量, 形成陀螺效應。該陀螺儀由三大組件組成: 殼體組件、 波紋管組件以及浮子組件。 浮子組件包括馬達和浮子框架,其內部充滿氣體。 在工作狀態(tài)下, 馬達轉子與定子通過高速的相對運動在兩者之間形成氣膜, 實現(xiàn)非接觸懸浮。 殼體上的法蘭環(huán)通過固定螺釘與鈦合金平臺固緊, 殼體與浮子之間的間隙充滿浮油, 殼體外周均勻纏繞電阻加熱絲, 浮子右轉子上同樣纏有發(fā)熱線圈, 通過調節(jié)這兩熱源的電流大小對陀螺儀進行溫度控制, 進而使浮油的浮力等于浮子的重力, 最大程度地減小摩擦, 提高陀螺儀精度［12］。 圖1 為液浮陀螺儀內部結構示意圖,由于陀螺儀內部結構復雜, 固體、 浮油和馬達內部氣體相互影響, 考慮陀螺儀內部散熱的同時又要保證浮子的重力等于浮油的浮力, 因此對液浮陀螺儀進行整體有限元仿真分析是陀螺儀設計過程中至關重要的一步。

圖1 液浮陀螺儀內部結構示意圖Fig.1 Schematic diagram of liquid floated gyroscope internal structure

2 仿真計算基礎

本文的仿真計算主要涉及陀螺儀內部流體的流場以及儀器整體的溫度場, 結合多物理場有限元軟件仿真過程中的實際情況與陀螺儀的工作原理, 現(xiàn)做以下假設:

1）對結構進行簡化, 忽略零件的倒角、 不重要的孔等特征, 將非關鍵組件如接線柱、 螺釘?shù)扰c基礎零件合為一體, 看成一個整體（假設1）;

2）由于陀螺儀的整體溫度范圍較低, 不考慮熱輻射的影響（假設2）;

3）只考慮流體的重力, 不考慮結構的重力、 熱變形和熱應力（假設3）。

根據(jù)以上假設, 建立如下的控制方程。

2.1 流體流動控制方程

在流動過程中, 流體需要始終滿足連續(xù)性方程、 動量守恒定律和能量守恒定律, 它們是陀螺儀內部流體進行有限元仿真的理論基礎［13-14］。

（1）連續(xù)性方程

式（1）中,ρ為流體密度,t為時間,Δ為拉普拉斯算子,U為流體的流速。

（2）動量守恒定律

式（2） 中,ui、uj為速度張量,p為靜壓力,xi、xj為坐標張量,τij為應力張量,ρgi為重力體積力,Fi為體積力。

（3）能量守恒定律

式（3）中,ht為流體的總能,ΔT為溫度差,λf為等效導熱系數(shù),Sh為外部熱源,τ為黏性流體的剪切應力。

2.2 陀螺儀內部傳熱方程

傳熱是指當有溫度梯度出現(xiàn)時所發(fā)生的能量轉移。 在自然界中, 熱傳導、 熱對流和熱輻射被稱為最基本的傳熱途徑, 任何傳熱問題都可由其中某種或多種組合而成。 基于假設2 以及陀螺儀實際工作情況, 其內部的傳熱途徑主要包括: 固體與固體之間的熱傳導、 固體與氦氣和浮油之間的對流換熱。以下是兩種傳熱途徑的數(shù)學描述, 它們是陀螺儀內部傳熱問題仿真計算的理論基礎［15-16］。

（1）陀螺儀內部固體間熱傳導

式（4）中,Q為熱流量,k為導熱系數(shù),ΔT為溫度差,A為傳熱面積。

（2）對流換熱

式（5）中,q為對流換熱量,Ch為對流換熱系數(shù),Tss為固體表面溫度,Tsf為接觸面流體溫度,c為流體熱容,u為近壁面流體流速,T+為接觸面無量綱溫度［17］。

3 有限元模型的建立

上述方程是仿真計算的理論基礎, 對于簡單情況可以通過上述公式直接求得, 但對于本文所研究的陀螺儀而言, 其結構不規(guī)律, 流體運動與傳熱過程十分復雜, 無法求得解析解。 因此, 需要借助有限元法對上一節(jié)中的方程進行處理, 方可求得結果。 液浮陀螺儀有限元模型的建立過程如下。

3.1 結構建模與網格劃分

基于假設1, 在多物理場有限元仿真軟件中建立簡化后的三維模型, 如圖2 所示。 其中, 重力方向為-z方向。

陀螺儀各零組件的材料類型及物理屬性直接決定了其溫度場, 由前文的流體流動控制方程以及傳熱方程可知, 在儀表的仿真計算過程中需要知道材料的密度、 動力黏度等參數(shù)。 經統(tǒng)計整理后, 計算中需要使用到的材料及其相關參數(shù)如表1所示。

表1 材料參數(shù)表Table 1 Parameters of material

由于三維模型結構復雜, 因此主要以自由四面體為基本網格單元對整個模型進行網格劃分。同時, 對關鍵部位進行局部網格加密, 使用多物理場有限元仿真軟件劃分網格, 得到浮油的網格劃分模型, 如圖3 所示。 得到網格的單元數(shù)為27320394 個, 平均網格質量為0.7。

圖3 浮油網格劃分Fig.3 Schematic diagram of oil slick meshing

3.2 流場及熱邊界條件設置

陀螺儀內部流場包括兩部分: 氣體流場和浮油流場。 由于高速旋轉的馬達定子會帶動周圍氣體運動, 基于流體密度、 動力黏度以及流速等參數(shù)通過下式可計算得到氣體的流場雷諾數(shù)（Re ＞2000）, 因此使用湍流模型進行計算。 浮油主要受溫度梯度影響發(fā)生流動, 流動速率較小, 通過下式可計算得到浮油的流場雷諾數(shù)（Re ＜2000）, 采用層流模型進行計算。

式（6） 中,ρ、μ為流體的密度和動力黏性系數(shù),U、L為流場的特征速度和特征長度。

液浮陀螺儀的熱源主要包括三部分: 馬達熱源、 主動溫控熱源和黏性損耗。 馬達熱功耗約為3.5W; 主動溫控熱源包括馬達線圈發(fā)熱、 浮子右轉子線圈發(fā)熱以及殼體外加熱片的保溫加熱, 各部分的熱功率如表2 所示; 黏性損耗產生的熱功率只需在有限元仿真軟件中打開, 系統(tǒng)會自動計算。

表2 熱源功率表Table 2 Power of heat sources

考慮到陀螺儀工作狀態(tài)下的外罩與大氣接觸,故可將陀螺儀外罩處的對流換熱視為大空間自然換熱。 根據(jù)傳熱學相關公式［17］, 可計算得到格拉曉夫數(shù)Gr 和努塞爾數(shù)Nu

式（8）、 式（9）中,g為重力加速度,a為空氣體積膨脹系數(shù), Δt為溫度差,v為動力黏度,l為特征尺度,D為特征直徑,L為特征長度。

進一步得到陀螺儀表面對流換熱系數(shù)

式（9）中,ε為空氣導熱系數(shù),L為特征長度。

3.3 計算條件

設置模型的初始溫度為實際工作環(huán)境溫度328K, 在啟動一段時間后, 陀螺儀內部流場及溫度場會趨于穩(wěn)定狀態(tài), 因此設置求解類型為凍結轉子（穩(wěn)態(tài)）, 同時為保證計算結果的精確程度, 相對容差設為0.001 作為收斂條件。

4 流場與溫度場計算結果

計算結束后累計迭代總次數(shù)為160 次, 通過溫度實驗測量得到某點的溫度為329K, 通過對多物理場有限元仿真模型的流場及溫度場進行穩(wěn)態(tài)計算, 得到該點的溫度為330.176K, 誤差為0.36%。兩者溫度相近, 誤差在可接受的范圍之內, 可認為陀螺儀溫度仿真分析的結果合理可信。 進一步針對模型網格進行優(yōu)化后, 網格細化20%, 陀螺儀整體溫度范圍如表3 所示, 表明仿真模型的網格單元加密20%后, 仿真溫度梯度變化量為0.05%。同時, 選取關鍵部位, 如圖4 所示, 得到其溫度曲線。 由圖4 可知, 仿真模型在關鍵部位的50μm 空間間隔內, 溫度梯度分辨率達到0.01K。

表3 網格細化前后溫度對比Table 3 Temperature comparison before and after grid refinement

圖4 關鍵部位及其溫度曲線Fig.4 Schematic diagram of key part and its temperature curve

4.1 氣體與浮油流場

由式（5）可知, 固體與流體之間的熱傳導需要知道流體的流速, 因此有必要對陀螺儀流體的流場進行仿真計算。 計算得到陀螺儀內部流體的流場如圖5 和圖6 所示。

圖5 氣體流場流速圖Fig.5 Flow velocity diagram of gas flow field

圖6 浮油流場流速圖Fig.6 Flow velocity diagram of oil slick flow field

由圖5 可知, 高速旋轉的馬達轉子帶動周圍的氣體一起運動, 形成渦流。 氣體的最大流速為67.40m/s, 在轉子距軸線半徑最大處。

由圖6 可知, 中段浮油即浮筒與殼體之間微小間隙內的浮油運動相對比較穩(wěn)定, 而左右段浮油在非均勻熱場中的運動十分復雜。 氣體的最大流速為3.04 ×10-6m/s, 在右端蓋過油孔處。

4.2 溫度場

對整表進行計算, 得到液浮陀螺儀垂直輸入軸方向的溫度分布切面, 如圖7 所示。 由圖7 可知, 儀表在正常工作條件下的溫度范圍為328K ～341.38K, 溫度最高點位于馬達定子與轉子高速運動時形成的氣膜處, 溫度最低點位于法蘭環(huán)與外界相接觸的位置。 馬達內部氣體的溫度云圖如圖8所示, 浮油的溫度云圖及部分關鍵位置的溫度如圖9 所示。 由圖8、 圖9 可知, 計算得到的氣體溫度范圍為338.43K ～341.32K, 溫差為2.89K; 浮油的溫度范圍為337.23K ～338.84K, 溫差為1.61K, 系統(tǒng)自動計算得到流體流動過程中產生的黏性損耗占總熱量的13.90%。 若不考慮馬達轉子與定子之間氣體體積變化產生的熱量, 陀螺儀整體溫度范圍為328K ～335.49K, 可知馬達內部氣體收縮膨脹產生的熱量約占總熱量的56%。 圖9 中的兩個曲線圖為靠近殼體部位的浮油溫度分布圖,上面的曲線圖位置為遠離法蘭環(huán)掛耳的位置, 法蘭環(huán)正對點溫度為338.43K; 下面的曲線圖位置為靠近法蘭環(huán)掛耳的位置, 法蘭環(huán)正對點溫度為338.29K。 從中可知在關鍵位置處, 溫度精度可達0.01K 以下, 且法蘭環(huán)掛耳位置對浮油溫度梯度影響較大, 靠近掛耳位置與遠離掛耳位置的溫度相差0.14K。

圖7 整表的溫度梯度圖Fig.7 Temperature gradient diagram of the whole liquid floated gyroscope

圖9 浮油溫度梯度云圖及關鍵位置溫度Fig.9 Nephogram of temperature gradient for oil slick and temperatures at key positions

5 結論

本文通過使用多物理場有限元仿真軟件, 建立了液浮陀螺儀整表的流固熱耦合多物理場計算模型, 充分考慮了陀螺儀固體與流體之間的相互作用, 結合實際工作環(huán)境, 得到模型準確的計算邊界條件, 最終計算得到陀螺儀內部溫度分布情況以及內部流體的運動傳熱狀態(tài), 為后續(xù)陀螺儀溫控系統(tǒng)設計提供依據(jù)。 具體可得到以下結論:

1）法蘭環(huán)掛耳對液浮陀螺儀浮筒外側浮油溫度梯度影響顯著, 可有針對性地設計控溫系統(tǒng);

2）流體在固體接觸表面由于摩擦產生的熱量約占總熱量的13.90%, 馬達內部氣體收縮膨脹產生的熱量約占總熱量的56%, 其余熱量為主動熱源產生;

3）浮油溫度梯度較明顯, 最高溫與最低溫相差1.61K, 會造成浮油密度分布不均, 引起浮油流動, 進而影響浮子外側的應力分布和陀螺儀精度。