黃 濱，王 璐，陳 磊，劉錦濤，吳大轉(zhuǎn)

1.浙江大學(xué)海洋學(xué)院，舟山 316021 2.北京控制工程研究所，北京 100190 3.浙江大學(xué)能源工程學(xué)院，杭州 310027

0 引 言

隨著航天事業(yè)的迅速發(fā)展，對航天器的有效工作壽命與經(jīng)濟性有了更高的要求，因而推進劑在航天器中的填充量在不斷地增大.現(xiàn)階段的航天任務(wù)越來越復(fù)雜，對其姿態(tài)和機動性能要求越來越高，使得我們對航天器貯箱運行過程中推進劑在不同工況下的流動狀態(tài)越來越重視[1-3].板式表面張力貯箱是國際上航行器主流的推進劑貯箱,它利用表面張力進行液體輸送和氣液分離，為發(fā)動機或推力器提供不夾氣的推進劑，PMD作為板式表面張力貯箱的核心部件對貯箱為發(fā)動機或推進器提供不夾氣的推進劑有關(guān)鍵作用[4-7].由于貯箱內(nèi)推進劑的晃動會對PMD的液體管理能力產(chǎn)生重要影響，從而對航天器的正常運行產(chǎn)生重大影響，因此需要關(guān)注航天器飛行過程中不同姿態(tài)下的推進劑晃動問題[8-12].針對微重力環(huán)境下板式貯箱內(nèi)流體行為研究，由于無法長時間進行地面微重力驗證試驗，常采用數(shù)值仿真進行先驗性分析.通過數(shù)值仿真可以對航天器航行中推進劑行為進行預(yù)測，分析貯箱結(jié)構(gòu)對推進劑管理性能的影響，對PMD結(jié)構(gòu)設(shè)計與優(yōu)化起導(dǎo)向作用[13-14].本文根據(jù)衛(wèi)星可能的在軌飛行情況，針對重定位過程中推進劑有明顯晃動產(chǎn)生的東西位置保持和南北位置保持工況下的板式表面張力貯箱內(nèi)部流體運動過程進行數(shù)值仿真，并對貯箱內(nèi)推進劑晃動行為和PMD液體管理性能進行分析，為后續(xù)微重力落塔試驗及空間站試驗提供參考.

1 數(shù)值模擬模型

本文采用VOF兩相流動模型模擬側(cè)向加速度環(huán)境下板式表面張力貯箱內(nèi)流體分布規(guī)律，對貯箱內(nèi)推進劑晃動現(xiàn)象進行分析研究，驗證板式表面張力貯箱的推進劑管理能力.

1.1 VOF計算模型

VOF兩相流模型是一種求解氣液兩相流動和自由表面流動的算法，它最早是由美國Los Alamos科學(xué)實驗室開發(fā)用于處理含有自由面的二維粘性流體的數(shù)值計算方法[15-20].該模型通過預(yù)測氣液兩相界面的自由表面流動來監(jiān)測側(cè)向微重力 環(huán)境下貯箱內(nèi)流體運動過程，由于推進劑在貯箱內(nèi)的流動過程雷諾數(shù)較小，流動形態(tài)一般為層流，此種情況下VOF模型的基本方程式由物性方程、連續(xù)方程、動量方程等三類方程組成.

1)物性方程

流體的物性是由流體中不同相的體積組分確定的，物性方程就是要確定不同體積組分時流體物性的表達式，提供給其他輸運方程.貯箱內(nèi)部流動只有氣液兩相，假設(shè)第二相的體積組分被跟蹤，那么流體的密度屬性方程為

ρ=α1ρ1+α2ρ2

(1)

α1+α2=1

(2)

式中，ρ為混合流體的密度，α1、α2分別為第一相和第二相的體積組分；ρ1、ρ2分別為第一相和第二相的密度，它們?yōu)榻o定值.

2)連續(xù)方程

流體連續(xù)方程的基本表達式如下：

(3)

式中，ui為混合流體速度，xi為混合流體位置，R為源項.

3)動量方程

混合流體的動量方程為

(4)

式中，p為貯箱壓力，ui、uj分別為兩相混合流體的氣相和液相速度，xi、xj分別為兩相混合流體的氣相和液相位置，t為時間，gj為微重力加速度，μ為粘性系數(shù).由于表面張力的作用，有

(5)

式中，σ為表面張力系數(shù)，r為液帶的曲率半徑，將式(5)代入(4)得混合流體的動量方程表達式

(6)

這就是混合流體的動量方程.

1.2 貯箱模型

以容積為884 L，內(nèi)直徑為1128 mm的板式表面張力貯箱作為研究對象，貯箱的PMD部件主要由4個導(dǎo)流板和1個蓄液器構(gòu)成，結(jié)構(gòu)簡圖1所示.貯箱軸向為Y軸，導(dǎo)流板布置在+X、-X、+Z、-Z四個方向上起到傳輸推進劑的作用，蓄液器位于貯箱下方出口處起到蓄液的功能.

圖1 板式表面張力貯箱三維模型Fig.1 Model of the vane-type surface tension tank

在ICEM中采用六面體結(jié)構(gòu)化網(wǎng)格分別對貯箱和蓄液器進行計算域離散.為便于進行網(wǎng)格劃分，對于貯箱和導(dǎo)流板區(qū)域采用1/4周期性流域進行網(wǎng)格劃分；對于蓄液器區(qū)域采用1/24周期性流域進行網(wǎng)格劃分.將周期性流域在ICEM中陣列出全域模型如圖2所示，貯箱殼體和導(dǎo)流板流域總體網(wǎng)格數(shù)約為291萬，蓄液器流域總體網(wǎng)格數(shù)約為309萬，經(jīng)過網(wǎng)格無關(guān)性驗證，具有良好的收斂性，計算模型如圖2中(c)圖所示.

圖2 貯箱各部件網(wǎng)格Fig.2 Mesh of components of the vane-type surface tension tank

1.3 模擬求解方案

(1)邊界條件

計算時將各個壁面設(shè)為固壁邊界，即為無穿透、無滑移邊界.

(2)初始體積組分

計算域由液體和氣體兩部分組成，設(shè)定氣相為主相，液相為次相.液體介質(zhì)為NTO，密度為1444 kg/m3(20℃)，粘性系數(shù)為4.2×10-4Pa·S；氣體介質(zhì)為氦氣，密度為0.1625 kg/m3(20℃)，粘性系數(shù)為1.99×10-5Pa·S.NTO的表面張力系數(shù)為0.026 N/m.在進行計算時需設(shè)定計算域內(nèi)的初始體積組分分布，本文按照推進劑的填充比(容量占整個體積的比率)5%、25%、45%、60%，初始時推進劑處于貯箱底部進行設(shè)置.

(3)計算模型

利用 Fluent18.0進行一階隱式非穩(wěn)態(tài)數(shù)值計算，時間步長為 1×10-3s.壓力-速度耦合方法采用PISO，其他方程的離散格式設(shè)置如下：壓力離散：Body Force Weighted；動量離散：QUICK，體積分?jǐn)?shù)離散：Compressive.方程收斂標(biāo)準(zhǔn)：質(zhì)量守恒方程殘差 RRR<10-3，所有速度量殘差：RRR<10-3.能量殘差RRR<10-3.

2 流體行為數(shù)值仿真結(jié)果與分析

側(cè)向加速度環(huán)境模擬了航行器航行于東西位置保持和南北位置保持的工況下貯箱的受力情況，推進劑的初始液面處于貯箱底部，側(cè)向加速度大小為0.015 m/s2，方向為X軸方向.

2.1 側(cè)向微重力重定位數(shù)值仿真

通過數(shù)值仿真模擬側(cè)向微重力重定位過程，得到推進劑重定位后的分布情況，預(yù)測了板式表面張力貯箱內(nèi)推進劑在側(cè)向加速環(huán)境下重定位過程中的流體行為.為了解不同推進劑填充比對重定位過程中推進劑晃動行為的影響，設(shè)定初始推進劑填充比分別為5%、25%、45%、60%.各填充比下貯箱初始氣液云圖如圖3所示，重定位后貯箱氣液分布如圖4、圖5所示.

圖3 各填充比初始氣液云圖Fig.3 Initial gas-liquid contours of different fill ratio

圖4 各填充比重定位結(jié)束氣液云圖Fig.4 Gas-liquid contours after reorientation of different fill ratio

圖5 各填充比重定位結(jié)束氣液交界面Fig.5 Gas-liquid interface after reorientation of different fill ratio

重定位結(jié)束氣液云圖表明氣體和液體在重定位過程中無混合，氣液界面穩(wěn)定后呈現(xiàn)波浪狀，蓄液器內(nèi)部液體未流失.氣液界面呈現(xiàn)波浪狀說明液體內(nèi)部存在相互作用力使得液體內(nèi)部產(chǎn)生晃動；但液體的晃動對氣液分離過程并沒有造成破壞、蓄液器仍能保持液體儲蓄，表明該板式表面張力貯箱的PMD在該側(cè)向加速的環(huán)境下能夠?qū)崿F(xiàn)推進劑的管理與運輸.氣液交界面圖可以更清晰顯示液體表面的波浪狀，液體在側(cè)向微重力的作用下將貯箱底部、中部、頂部的側(cè)表面浸染.

2.2 重定位過程流量分析

為更清晰地了解板式貯箱側(cè)向加速度下重定位過程中推進劑體運動過程，對板式貯箱初始氣液交界面和豎直對稱面進行流量監(jiān)測，得到重定位過程中液體的晃動特性.流量監(jiān)測面設(shè)置示意圖如圖6所示，流量監(jiān)測曲線如圖7、圖8所示.

圖6 流量監(jiān)測面設(shè)置示意圖Fig.6 Schematic diagram of flow monitoring surface setting

通過液體晃動氣液云圖周期變化可以更直觀的對流量時間曲線的變化趨勢進行解釋，以填充比為25%時液體晃動周期圖為例如圖9、圖10所示，并將氣液云圖對應(yīng)的時刻點在相應(yīng)的流量時間曲線上進行了標(biāo)注如圖7(b)和圖8(b)所示.

圖7 各填充比交界面流量時間曲線Fig.7 Flow-time curve of interface flow of different fill ratio

圖8 各填充比垂直對稱面流量時間曲線Fig.8 Flow-time curves of vertical symmetry plane flow with different fill ratio

圖9 填充比25%時液體晃動周期氣液云圖Fig.9 Gas-liquid contour of sloshing cycle of fill ratio 25%

圖10 填充比25%時液體晃動周期氣液交界面Fig.10 Gas-liquid interface of sloshing cycle of fill ratio 25%

初始交界面和豎直對稱面的流量時間曲線顯示在重定位過程中流量的變化存在明顯的波動，通過液體晃動周期氣液云圖可以更直觀的看出液體的晃動行為.通過對比不同填充比的流量時間曲線發(fā)現(xiàn)：液體晃動程度與推進劑的填充比存在一定的聯(lián)系.隨著填充比的增加，重定位完成的時間與流量波動的最大幅值的變化呈現(xiàn)先增大后減小的趨勢.該現(xiàn)象表明隨填充比的增加液體的晃動劇烈程度呈拋物線狀的變化趨勢，且存在峰值狀態(tài)，為保證貯箱的正常運行和PMD的推進劑管理能力，應(yīng)在峰值區(qū)間內(nèi)進行密集試驗.

3 結(jié) 論

本文設(shè)置X軸方向0.015 m/s2的加速度環(huán)境模擬在軌飛行時東西位置保持和南北位置保持航行工況下的側(cè)向微重力環(huán)境.對不同填充比表面張力板式貯箱重定位過程進行數(shù)值模擬，得到側(cè)向微重力環(huán)境下重定位過程中推進劑的運動過程：氣液在重定位過程中無混合且液體能夠穩(wěn)定積蓄在貯箱的一側(cè)，蓄液器能夠有效的儲蓄推進劑，氣液交界面呈波浪狀.重定位結(jié)果表明該板式表面貯箱在該側(cè)向微重力環(huán)境下可以有效地實現(xiàn)氣液分離，雖然推進劑內(nèi)部存在明顯的晃動，但貯箱內(nèi)PMD能夠有效的管理貯箱內(nèi)部液體.流量時間曲線顯示貯箱內(nèi)液體的晃動與液體的填充比存在一定的關(guān)聯(lián)，在環(huán)境一定的情況下，重定位的時間和流量波動的最大幅值隨流量的增加呈先增加后減小的趨勢.本次模擬結(jié)果對板式表面張力貯箱的側(cè)向加速度下的運行有了一定的預(yù)測，對后續(xù)即將開展的微重力落塔試驗和空間站試驗提供了重要的參考依據(jù).