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

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

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

0 引言

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

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

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

1 液浮陀螺儀工作原理

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

圖1 液浮陀螺儀內(nèi)部結(jié)構(gòu)示意圖Fig.1 Schematic diagram of liquid floated gyroscope internal structure

2 仿真計(jì)算基礎(chǔ)

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

1）對(duì)結(jié)構(gòu)進(jìn)行簡化, 忽略零件的倒角、 不重要的孔等特征, 將非關(guān)鍵組件如接線柱、 螺釘?shù)扰c基礎(chǔ)零件合為一體, 看成一個(gè)整體（假設(shè)1）;

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

3）只考慮流體的重力, 不考慮結(jié)構(gòu)的重力、 熱變形和熱應(yīng)力（假設(shè)3）。

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

2.1 流體流動(dòng)控制方程

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

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

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

（2）動(dòng)量守恒定律

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

（3）能量守恒定律

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

2.2 陀螺儀內(nèi)部傳熱方程

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

（1）陀螺儀內(nèi)部固體間熱傳導(dǎo)

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

（2）對(duì)流換熱

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

3 有限元模型的建立

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

3.1 結(jié)構(gòu)建模與網(wǎng)格劃分

基于假設(shè)1, 在多物理場(chǎng)有限元仿真軟件中建立簡化后的三維模型, 如圖2 所示。 其中, 重力方向?yàn)?z方向。

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

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

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

圖3 浮油網(wǎng)格劃分Fig.3 Schematic diagram of oil slick meshing

3.2 流場(chǎng)及熱邊界條件設(shè)置

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

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

液浮陀螺儀的熱源主要包括三部分: 馬達(dá)熱源、 主動(dòng)溫控?zé)嵩春宛ば該p耗。 馬達(dá)熱功耗約為3.5W; 主動(dòng)溫控?zé)嵩窗R達(dá)線圈發(fā)熱、 浮子右轉(zhuǎn)子線圈發(fā)熱以及殼體外加熱片的保溫加熱, 各部分的熱功率如表2 所示; 黏性損耗產(chǎn)生的熱功率只需在有限元仿真軟件中打開, 系統(tǒng)會(huì)自動(dòng)計(jì)算。

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

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

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

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

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

3.3 計(jì)算條件

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

4 流場(chǎng)與溫度場(chǎng)計(jì)算結(jié)果

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

表3 網(wǎng)格細(xì)化前后溫度對(duì)比Table 3 Temperature comparison before and after grid refinement

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

4.1 氣體與浮油流場(chǎng)

由式（5）可知, 固體與流體之間的熱傳導(dǎo)需要知道流體的流速, 因此有必要對(duì)陀螺儀流體的流場(chǎng)進(jìn)行仿真計(jì)算。 計(jì)算得到陀螺儀內(nèi)部流體的流場(chǎng)如圖5 和圖6 所示。

圖5 氣體流場(chǎng)流速圖Fig.5 Flow velocity diagram of gas flow field

圖6 浮油流場(chǎng)流速圖Fig.6 Flow velocity diagram of oil slick flow field

由圖5 可知, 高速旋轉(zhuǎn)的馬達(dá)轉(zhuǎn)子帶動(dòng)周圍的氣體一起運(yùn)動(dòng), 形成渦流。 氣體的最大流速為67.40m/s, 在轉(zhuǎn)子距軸線半徑最大處。

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

4.2 溫度場(chǎng)

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

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

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

5 結(jié)論

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

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

2）流體在固體接觸表面由于摩擦產(chǎn)生的熱量約占總熱量的13.90%, 馬達(dá)內(nèi)部氣體收縮膨脹產(chǎn)生的熱量約占總熱量的56%, 其余熱量為主動(dòng)熱源產(chǎn)生;

3）浮油溫度梯度較明顯, 最高溫與最低溫相差1.61K, 會(huì)造成浮油密度分布不均, 引起浮油流動(dòng), 進(jìn)而影響浮子外側(cè)的應(yīng)力分布和陀螺儀精度。