賀志遠(yuǎn) 胡文林 馬貴軍

（1.海軍航空大學(xué) 煙臺(tái) 264001）（2.中國(guó)人民解放軍91049部隊(duì) 青島 266001）

1 引言

FLUENT軟件是目前世界上最先進(jìn)的流體動(dòng)力仿真軟件，可以解決大部分流體問(wèn)題。但是在液面晃動(dòng)問(wèn)題中的應(yīng)用相對(duì)較少，主要有兩個(gè)原因：一是液面的晃動(dòng)問(wèn)題屬于非穩(wěn)態(tài)流動(dòng)，在非慣性坐標(biāo)系下要遵循物理守恒定律（質(zhì)量守恒定律和動(dòng)量守恒定律），而動(dòng)量守恒定律是在慣性系下建立的，在非慣性系下的形式是不同的［1～3］，所以無(wú)法在軟件中直接利用Navier-Stokes方程來(lái)求解；二是容器內(nèi)部液體受到質(zhì)量力和重力共同作用，如果容器加速度隨時(shí)間變化，它所受的質(zhì)量力也將不停變換，無(wú)法選定合適的參考系，導(dǎo)致應(yīng)用起來(lái)比較困難［4～6］。由于軟件在應(yīng)用中存在這些問(wèn)題，許多專家、學(xué)者利用不同的方法和理論對(duì)該問(wèn)題進(jìn)行了廣泛的研究。其中尚春雨，趙金城［7］利用等加速度直線運(yùn)動(dòng)水體、矩形水箱基本自振等算例，并結(jié)合Tait［8］等的研究成果，驗(yàn)證了該軟件在解決剛性容器內(nèi)頁(yè)面晃動(dòng)問(wèn)題的可行性。本文針對(duì)導(dǎo)彈推進(jìn)劑貯箱在使用過(guò)程中晃動(dòng)的問(wèn)題，使用FLUENT進(jìn)行仿真分析，根據(jù)模擬計(jì)算結(jié)果，給出推進(jìn)劑貯箱使用的風(fēng)險(xiǎn)評(píng)估建議。

2 FLUENT軟件簡(jiǎn)介

2.1 計(jì)算流體力學(xué)基礎(chǔ)

計(jì)算流體力學(xué)分析的基本思想是選用有限個(gè)離散點(diǎn)上的物理場(chǎng)變量值的集合，通過(guò)一定的法和原則，來(lái)代替流體在時(shí)間系及空間系上原本連續(xù)的物理量場(chǎng)（壓力場(chǎng)和速度場(chǎng)等），從而建立起關(guān)于這些有限個(gè)離散點(diǎn)上場(chǎng)變量之間關(guān)系的代數(shù)方程組，然后運(yùn)用有限元的基本方法求解代數(shù)方程組，進(jìn)而獲得流體在時(shí)間系和空間系上個(gè)物理場(chǎng)變量的相對(duì)近似值，最終得出具體結(jié)論［9］。

2.2 湍流模型

1972 年，Launder和 Spalding［11］在雷諾方程的基礎(chǔ)上，提出了標(biāo)準(zhǔn)k-ε湍流模型。其中湍流耗散方程是由經(jīng)驗(yàn)公式推導(dǎo)出來(lái)的方程，湍流動(dòng)能方程是精確方程，因此這是個(gè)半經(jīng)驗(yàn)公式。同時(shí)，k-ε模型的使用需要一定的前提條件，需要忽略流體分子之間的粘性特性，假定流場(chǎng)是完全發(fā)展的湍流場(chǎng)，因此標(biāo)準(zhǔn)k-ε模型主要針對(duì)完全湍流的流場(chǎng)［10］。

3 仿真推進(jìn)劑貯箱晃動(dòng)問(wèn)題

3.1 推進(jìn)劑貯箱仿真模型的基本假設(shè)

對(duì)于貯箱內(nèi)推進(jìn)劑晃動(dòng)問(wèn)題的完全模擬是很難實(shí)現(xiàn)的，為方便研究，現(xiàn)做出如下假設(shè)和簡(jiǎn)化。

1）貯箱為剛性容器，運(yùn)動(dòng)過(guò)程中不發(fā)生形變，不與內(nèi)部推進(jìn)劑發(fā)生反應(yīng)；

2）貯箱內(nèi)推進(jìn)劑運(yùn)動(dòng)過(guò)程中不發(fā)生相變變化，不發(fā)生化學(xué)變化；

3）貯箱內(nèi)壓強(qiáng)保持不變；

4）貯箱簡(jiǎn)化為球體，忽略不規(guī)則部分；

5）貯箱運(yùn)動(dòng)方式為簡(jiǎn)化后的運(yùn)動(dòng)；

6）忽略摩擦及熱量交換。

3.2 液體推進(jìn)劑非靜態(tài)過(guò)程的基本方程

液體推進(jìn)劑非靜態(tài)過(guò)程主要是流體的流動(dòng)，流體流動(dòng)的基本守恒定律包括：質(zhì)量守恒定律、動(dòng)量守恒定律、能量守恒定律，本問(wèn)題中由于涉及不同組分間的相互作用，還需要加上組分守恒定律和湍流輸運(yùn)方程，具體方程如下。

質(zhì)量守恒方程：

動(dòng)量守恒方程：

式中：p為靜壓；τij為應(yīng)力張量；Fi和gi分別為i方向上的外部體積力和重力體積力。

能量守恒方程：

式中：cp為比熱容；T為溫度；k為流體傳熱系數(shù)；ST為流體的粘性耗散項(xiàng)。

湍流輸運(yùn)方程：

式中：為組分的擴(kuò)散通量；Ri為第i種物質(zhì)的化學(xué)反應(yīng)凈剩量；Si為用戶自定義源項(xiàng)和離散項(xiàng)的產(chǎn)生率。

3.3 模型的選擇

在該仿真軟件中主要提供了三種多相流模型：VOF模型［12］（Volume of Fluid Model）、混合模型（混合Model）、歐拉模型（歐拉Model），本文研究的問(wèn)題采用VOF模型。

3.4 液體推進(jìn)劑非靜態(tài)過(guò)程的數(shù)值模擬

根據(jù)前文所做分析，在使用過(guò)程中，導(dǎo)彈會(huì)受到不同程度的外部沖擊，進(jìn)而帶動(dòng)推進(jìn)劑作各種無(wú)規(guī)則運(yùn)動(dòng)。在本節(jié)中我們將對(duì)貯箱無(wú)規(guī)則運(yùn)動(dòng)進(jìn)行數(shù)值模擬。

1）幾何模型的建立

利用FLUENT的ICEM前處理軟件，對(duì)液體推進(jìn)劑兩種非靜態(tài)過(guò)程建立幾何模型。

2）計(jì)算網(wǎng)格的劃分

由于實(shí)際工程計(jì)算中存在較多復(fù)雜且不規(guī)則的區(qū)域，要計(jì)算這些區(qū)域的流體力學(xué)，就需要?jiǎng)澐志W(wǎng)格，計(jì)算的效率和計(jì)算結(jié)果的精確度往往取決于網(wǎng)格劃分的質(zhì)量。

本文主要模擬推進(jìn)劑貯箱不規(guī)則運(yùn)動(dòng)過(guò)程中液體推進(jìn)劑的流動(dòng)，要關(guān)注推進(jìn)劑運(yùn)動(dòng)過(guò)程中對(duì)貯箱壁的作用力，同時(shí)，加注口與貯箱連接部分的特殊性，不可忽視，因此采取非均勻網(wǎng)格劃分，在貯箱和加注口連接的地方需要加密網(wǎng)格，從而得到更加精確的結(jié)果。

我們初步設(shè)定網(wǎng)格單元尺寸設(shè)定為10mm，同時(shí)考慮到加注口較窄且連接部分需要重點(diǎn)研究，該部分進(jìn)行網(wǎng)格加密，所選網(wǎng)格單元尺寸為1mm，經(jīng)過(guò)網(wǎng)格劃分軟件的劃分，兩種模型最終的網(wǎng)格數(shù)量為388204個(gè)。

3）初始及邊界條件的設(shè)定

計(jì)算中，液相部分分別為綠色四氧化二氮和偏二甲肼，氣相為氮?dú)?，工作溫度選擇室溫（20℃），工質(zhì)具體數(shù)據(jù)如表1。

4）仿真計(jì)算條件和方法

（1）求解器的選擇

在FLUENT軟件的仿真計(jì)算中，解決問(wèn)題的關(guān)鍵是選擇合適的求解器，合理選擇求解器將使計(jì)算過(guò)程更高效、計(jì)算結(jié)果更準(zhǔn)確。根據(jù)本文解決問(wèn)題的特點(diǎn)，采用基于壓力的求解器。

表1 工質(zhì)參數(shù)表

（2）時(shí)間步長(zhǎng)的設(shè)定

時(shí)間步長(zhǎng)（Time Step Size）的設(shè)定會(huì)直接影響求解過(guò)程的穩(wěn)定性，如果時(shí)間步長(zhǎng)選擇不合理，與網(wǎng)格尺寸不相匹配，則會(huì)導(dǎo)致計(jì)算結(jié)果發(fā)散，難以收斂。本節(jié)根據(jù)Von Neumann法則［13］來(lái)對(duì)時(shí)間步長(zhǎng)進(jìn)行設(shè)置，該法則的基本原則是：計(jì)算過(guò)程中中引入誤差，如果其強(qiáng)度（或振幅）隨時(shí)間的推移而衰減或保持不變，則這種格式是穩(wěn)定的，也就意味著時(shí)間步長(zhǎng)的設(shè)定是合理的，反之是不合理的。

由于本文要研究一個(gè)長(zhǎng)期的過(guò)程，時(shí)間步長(zhǎng)采用分段設(shè)計(jì)的方式。結(jié)合這一法則，在提高精確度又不增加巨大的運(yùn)算量的前提下，在前100天的仿真中，時(shí)間步長(zhǎng)選擇60s；100天以后至15年的過(guò)程中，時(shí)間步長(zhǎng)選擇600000s（近似7天），時(shí)間步長(zhǎng)初步設(shè)定，在仿真過(guò)程中還可能適當(dāng)修改。

（3）松弛因子的選擇

在FLUENT軟件的仿真計(jì)算過(guò)程中，松弛因子通過(guò)影響迭代次數(shù)進(jìn)而影響計(jì)算的穩(wěn)定性、效率和求解速度。在本文研究的問(wèn)題中，綠色四氧化二氮貯箱的運(yùn)動(dòng)方式較為簡(jiǎn)單，計(jì)算過(guò)程中不需要修改松弛因子，可以保證計(jì)算很好的收斂。而偏二甲肼貯箱運(yùn)動(dòng)方式復(fù)雜，內(nèi)部推進(jìn)劑晃動(dòng)較大，經(jīng)過(guò)初步驗(yàn)證，默認(rèn)的松弛因子很難使計(jì)算結(jié)果正確收斂，結(jié)合參考文獻(xiàn)，修改壓強(qiáng)、動(dòng)量、湍流動(dòng)能和湍流耗散率松弛因子的默認(rèn)值，分別修改為為0.2、0.8、0.6和0.6。

4 液體推進(jìn)劑非靜態(tài)過(guò)程仿真結(jié)果與分析

根據(jù)本文的設(shè)計(jì)思路，對(duì)于綠色四氧化二氮和偏二甲肼貯箱模型的運(yùn)動(dòng)仿真，分別選取時(shí)間為7天、28天、100天、300天、3年和15年這幾個(gè)時(shí)間節(jié)點(diǎn)的仿真數(shù)據(jù)，使用Fluent中的CFD-post軟件轉(zhuǎn)換為壓力云圖。其中，15年時(shí)間的壓力云圖如圖2～3所示。

選取貯箱與加注口的連接部分作為采樣點(diǎn)，采集仿真時(shí)間點(diǎn)7天、28天、100天、300天、3年和15年的壓力變化情況，如圖3所示。

圖1 四氧化二氮貯箱壓力云圖

圖2 偏二甲肼貯箱壓力云圖

圖3 貯箱連接部分壓力變化

從上面的壓力云圖可以看出，藍(lán)色部分為氮?dú)?，黃色和紅色部分為模擬的液體推進(jìn)劑，紅色越深代表壓力越大。因?yàn)檫\(yùn)動(dòng)的不規(guī)則性，每幅圖中都有不同的紅色區(qū)域，此區(qū)域表示壓力集中區(qū)，由具體數(shù)值可以發(fā)現(xiàn)，壓力較大區(qū)域超過(guò)正常區(qū)域兩倍以上，焊點(diǎn)附近壓力要比周圍大，所承受的沖擊力也大。從圖3中可以看出，連接口的壓力隨著仿真時(shí)間增長(zhǎng)而上升，在導(dǎo)彈運(yùn)輸過(guò)程中的顛簸和航貯過(guò)程中洋流的影響下，這種情況十分常見。在真實(shí)的推進(jìn)劑貯箱中，貯箱與加注口是采用焊接連接的，在連接部分存在大量焊點(diǎn)，與整體結(jié)構(gòu)不同，這些焊點(diǎn)部位強(qiáng)度不大，當(dāng)成為應(yīng)力集中點(diǎn)后，長(zhǎng)期作用下，很容易發(fā)生破壞，進(jìn)而導(dǎo)致液體推進(jìn)劑泄露。

5 結(jié)語(yǔ)

本文采用FLUENT仿真軟件模擬了導(dǎo)彈推進(jìn)劑貯箱模型的晃動(dòng)問(wèn)題，模擬液體推進(jìn)劑的非靜態(tài)狀態(tài)，與導(dǎo)彈使用過(guò)程中，實(shí)際的操作使用情況也是一一相對(duì)應(yīng)的。仿真結(jié)果表明，不同運(yùn)動(dòng)方式下，推進(jìn)劑對(duì)貯箱有不同的作用力，局部作用力明顯使得貯箱存在被破壞的風(fēng)險(xiǎn)，也就意味著液體推進(jìn)劑存在著泄漏的風(fēng)險(xiǎn)，對(duì)于導(dǎo)彈裝備維護(hù)保養(yǎng)有一定的指導(dǎo)意義。但是，由于仿真的時(shí)間跨度較大，過(guò)程中存在誤差積累現(xiàn)象，對(duì)于造成誤差的原因以及誤差修正還需要進(jìn)一步研究。