李青，黃俊，繆遠(yuǎn)明，于杭健

北京空間飛行器總體設(shè)計(jì)部，北京 100094

1 引言

運(yùn)載火箭、人造地球衛(wèi)星、深空探測(cè)器、載人空間站等現(xiàn)代航天器通常都是充液系統(tǒng)，攜帶著大量的液體推進(jìn)劑以及其他液體工質(zhì)，其中的液體質(zhì)量均占總質(zhì)量的一半以上。充液系統(tǒng)可分為全充液系統(tǒng)和部分充液系統(tǒng)。全充液系統(tǒng)中的液體無(wú)自由液面，例如加壓狀態(tài)的囊式貯箱、金屬膜片貯箱或滿箱狀態(tài)的表面張力貯箱等；而部分充液系統(tǒng)中的液體由于自由液面的運(yùn)動(dòng)，還存在晃動(dòng)的問題，例如充液比(即充液體積與貯箱容積之比)小于1的表面張力貯箱或未加壓的金屬膜片貯箱等。

對(duì)于充液系統(tǒng)中液體晃動(dòng)問題的研究始于20世紀(jì)60年代，研究方法主要有解析法、數(shù)值法和試驗(yàn)法。Abramson[1]系統(tǒng)介紹了幾種特殊形狀(包括矩形、圓柱形、圓環(huán)形等)容器內(nèi)液體橫向小幅晃動(dòng)的基本理論和解析解。Levin[2]給出了直角圓錐容器內(nèi)液體橫向晃動(dòng)的解析解。Hasheminejad等[3]分析了橫放圓柱形容器內(nèi)液體橫向晃動(dòng)問題。實(shí)際上，只有少數(shù)幾種形狀容器內(nèi)的液體橫向晃動(dòng)問題可進(jìn)行解析求解。Aslam[4]采用有限元法研究了地震導(dǎo)致的軸對(duì)稱容器內(nèi)的液體晃動(dòng)；Wang等[5]在對(duì)液體晃動(dòng)模態(tài)進(jìn)行有限元數(shù)值分析的基礎(chǔ)上，研究了晃動(dòng)阻尼的計(jì)算方法；Gedikli等[6]采用邊界元法進(jìn)行液體晃動(dòng)模態(tài)分析，并應(yīng)用該方法研究了帶防晃板貯箱內(nèi)的液體晃動(dòng)問題。計(jì)算流體動(dòng)力學(xué)(CFD)方法是一類求解Navier-Stokes方程組的數(shù)值方法，包括有限差分法[7]、有限體積法[8-9]等。試驗(yàn)法[10-13]是驗(yàn)證數(shù)值計(jì)算正確性的有效途徑。歐空局于2005年發(fā)射了一顆液體晃動(dòng)試驗(yàn)衛(wèi)星[14]——Sloshsat FLEVO，用于驗(yàn)證其基于有限體積法開發(fā)的CFD數(shù)值計(jì)算程序Comflo的有效性[15]。He等[16]采用有限體積法研究了油罐車內(nèi)自由油面在剎車、加速、減速等工況下產(chǎn)生的晃動(dòng)壓力，并通過了試驗(yàn)驗(yàn)證。在航天器動(dòng)力學(xué)分析和控制系統(tǒng)工程設(shè)計(jì)中，通常采用簡(jiǎn)單高效的等效力學(xué)模型計(jì)算來(lái)替代復(fù)雜費(fèi)時(shí)的晃動(dòng)流場(chǎng)計(jì)算，如單擺模型[17-18]、彈簧-質(zhì)量模型[19]、MPB(運(yùn)動(dòng)脈動(dòng)球)模型[20-21]等，從而可以比較容易地將液體晃動(dòng)動(dòng)力學(xué)方程整合到航天器系統(tǒng)動(dòng)力學(xué)方程中。

無(wú)論是全充液系統(tǒng)還是部分充液系統(tǒng)，當(dāng)研究其轉(zhuǎn)動(dòng)或姿態(tài)動(dòng)力學(xué)時(shí)，必須對(duì)其慣性張量進(jìn)行計(jì)算。Dodge[22]給出了矩形貯箱內(nèi)液體轉(zhuǎn)動(dòng)慣量的解析解以及圓柱形貯箱內(nèi)液體轉(zhuǎn)動(dòng)慣量的數(shù)值解曲線；馬斌捷等[23-24]分析整理了航天器圓環(huán)筒形平底貯箱內(nèi)液體有效轉(zhuǎn)動(dòng)慣量計(jì)算方法，研究了縱向擋板和橫向擋板對(duì)液體轉(zhuǎn)動(dòng)慣量的影響。Lee[25]采用半解析方法計(jì)算了橢圓形、矩形、六邊形和八邊形全充液貯箱內(nèi)液體的轉(zhuǎn)動(dòng)慣量。

綜上，關(guān)于充液系統(tǒng)的已有研究文獻(xiàn)主要側(cè)重于液體晃動(dòng)問題，對(duì)于液體質(zhì)量特性的計(jì)算，特別是慣性張量的計(jì)算研究較少，且主要局限于若干特殊形狀的貯箱，計(jì)算方法沒有通用性。實(shí)際上，充液航天器質(zhì)量特性的計(jì)算在工程中還沒有得到足夠的重視。在航天器設(shè)計(jì)中，通常將液體當(dāng)作等量等分布的固體計(jì)算慣性張量，這種固化液體的簡(jiǎn)化計(jì)算方法顯然是不符合實(shí)際情況的，會(huì)加大質(zhì)量特性參數(shù)的偏差范圍，從而增加動(dòng)力學(xué)模型的誤差和控制系統(tǒng)設(shè)計(jì)的難度。本工作采用勢(shì)流理論推導(dǎo)了任意形狀貯箱內(nèi)的液體質(zhì)量特性，采用有限元方法建立了計(jì)算液體質(zhì)量特性的數(shù)值算法，通過一個(gè)矩形貯箱的算例驗(yàn)證了算法的正確性，將該方法應(yīng)用于充液航天器質(zhì)量特性的計(jì)算，研究了不同的液體處理方法對(duì)計(jì)算結(jié)果的影響，以期為充液航天器質(zhì)量特性計(jì)算方法的選用提供參考。

2 理論推導(dǎo)

如圖 1所示，O0-XYZ為慣性坐標(biāo)系，O-xyz為充液系統(tǒng)本體坐標(biāo)系，Ω為液體內(nèi)部，Sf為自由液面，Sw為貯箱壁面，g為重力加速度。

圖 1 充液系統(tǒng)示意Fig.1 Illustration for liquid-filled system

當(dāng)充液系統(tǒng)運(yùn)動(dòng)時(shí)，相對(duì)O點(diǎn)矢徑為r的液體質(zhì)點(diǎn)的絕對(duì)速度v等于該點(diǎn)牽連速度與相對(duì)速度u之和，即

v=vO+ω×r+u

(1)

式中：vO和ω分別為O點(diǎn)速度和系統(tǒng)的剛體角速度。假設(shè)貯箱內(nèi)的液體是無(wú)粘、無(wú)旋、不可壓縮的，則根據(jù)勢(shì)流理論存在速度勢(shì)Φ，使得v=?Φ；在重力或慣性力加速度占主導(dǎo)、液體晃動(dòng)幅度較小、液體表面張力可忽略時(shí)，Φ在Ω內(nèi)滿足

ΔΦ=0

(2)

Φ在Sw上滿足

(3)

式中：n為單位外法線矢量。Φ在Sf上滿足

(4)

式中：t為時(shí)間。Φ具有如下形式：

(5)

式中：ψ為Stokes-Zhukovskiy矢勢(shì)；φi為第i階液體晃動(dòng)模態(tài)函數(shù)；qi為第i階液體晃動(dòng)模態(tài)廣義坐標(biāo)。

Δφi=0

(6)

φi在Sw上滿足

(7)

φi在Sf上滿足

(8)

由方程(6)～(8)可解出φi及相應(yīng)的固有頻率ωi。

ψ在Ω內(nèi)滿足

Δψ=0

(9)

ψ在Sw∪Sf上滿足

(10)

假設(shè)液體的密度為ρ，液體的動(dòng)能可以寫成

(11)

將式(1)代入式(11)，并根據(jù)拉格朗日方程推導(dǎo)出

(12)

(13)

(14)

式中：m為液體總質(zhì)量；rc為液體質(zhì)心位置矢量；CTi為第i階液體晃動(dòng)與貯箱的平動(dòng)耦合系數(shù)矢量；Fd為液體對(duì)貯箱的動(dòng)作用力矢量；J為液體等效剛體關(guān)于O點(diǎn)的慣性張量；CRi為第i階液體晃動(dòng)與貯箱的轉(zhuǎn)動(dòng)耦合系數(shù)矢量；Md為液體對(duì)貯箱關(guān)于O點(diǎn)的動(dòng)作用力矩矢量；μi為第i階液體晃動(dòng)廣義質(zhì)量。它們的表達(dá)式如下：

(15)

(16)

(17)

(18)

(19)

(20)

對(duì)于全充液系統(tǒng)，CTi、CRi和μi均為零，J即為液體慣性張量的理論表達(dá)式，由于沒有自由液面，式(17)中略去Sf項(xiàng)。式(15)～(17)給出了全充液系統(tǒng)內(nèi)液體質(zhì)量特性的表達(dá)式。

(21)

因此，部分充液系統(tǒng)內(nèi)液體慣性張量的理論表達(dá)式為

(22)

式(15)(16)(22)給出了部分充液系統(tǒng)內(nèi)液體質(zhì)量特性的表達(dá)式。

3 數(shù)值算法

采用三維有限元?jiǎng)澐忠后w域網(wǎng)格，由第2節(jié)的理論推導(dǎo)可知，充液系統(tǒng)內(nèi)液體的質(zhì)量和質(zhì)心位置可根據(jù)式(15)(16)進(jìn)行數(shù)值積分得到，而慣性張量與ωi、φi和ψ有關(guān)，需要首先采用有限元方法得到它們的數(shù)值解。

方程(6)～(8)和方程(9)(10)可分別寫成如下變分形式：

(23)

(24)

采用三維單元?jiǎng)澐忠后w網(wǎng)格，在每個(gè)單元內(nèi)場(chǎng)函數(shù)可以表示成

(25)

(26)

(27)

(28)

對(duì)式(27)(28)分別應(yīng)用變分原理，得到

(29)

(30)

求解方程(29)的廣義特征值問題得到ωi和φi的數(shù)值解。求解方程(30)的線性代數(shù)方程組得到ψ的數(shù)值解。然后再根據(jù)式(17)(22)進(jìn)行數(shù)值積分，分別得到全充液系統(tǒng)和部分充液系統(tǒng)的慣性張量。

4 驗(yàn)證算例

如圖2所示，對(duì)于長(zhǎng)為a、寬為b、液高為h的矩形貯箱液體，本體坐標(biāo)系原點(diǎn)O位于液體質(zhì)心，根據(jù)文獻(xiàn)[22]的推導(dǎo)，其液體轉(zhuǎn)動(dòng)慣量的解析解為

圖 2 矩形貯箱示意Fig.2 Illustration for a rectangular tank

(31)

(32)

(33)

若存在自由液面，則液體轉(zhuǎn)動(dòng)慣量的解析解為

(34)

(35)

Iz=Jz

(36)

在本驗(yàn)證算例中，取a=2.11 m、b=1.33 m、h=1.1 m，液體為水。如圖 3所示，液體有限元模型劃分為20 400個(gè)十結(jié)點(diǎn)四面體單元，總結(jié)點(diǎn)數(shù)為29 842。對(duì)于部分充液狀態(tài)，在計(jì)算Iy時(shí)截取了如圖 4所示的前3階x向液體晃動(dòng)模態(tài)，在計(jì)算Ix時(shí)截取了如圖 5所示的前3階y向液體晃動(dòng)模態(tài)。

圖3 矩形貯箱內(nèi)液體的有限元模型Fig.3 Finite element model for the liquid in a rectangular tank

圖4 沿x軸方向的前3階液體晃動(dòng)模態(tài)振型Fig.4 First three sloshing mode shapes along x-axis

圖5 沿y軸方向的前三階液體晃動(dòng)模態(tài)振型Fig.5 First three sloshing mode shapes along y-axis

表 1對(duì)比了液體質(zhì)量特性的數(shù)值解與解析解的計(jì)算結(jié)果，可見數(shù)值解與解析解的相對(duì)誤差非常小，從而驗(yàn)證了本文數(shù)值計(jì)算方法的正確性。由于矩形貯箱內(nèi)液體關(guān)于其質(zhì)心的慣性主軸始終為x軸、y軸和z軸，故所有的慣性積計(jì)算結(jié)果均為0且未在表1中列出。

表1 液體質(zhì)量特性的數(shù)值解與解析解對(duì)比

5 應(yīng)用研究

某充液航天器及其坐標(biāo)系定義如圖 6所示，OP-XPYPZP為機(jī)械坐標(biāo)系，OPC-XPCYPCZPC為質(zhì)心坐標(biāo)系(由機(jī)械坐標(biāo)系繞XP軸旋轉(zhuǎn)45°并平移得到)，OPC為航天器質(zhì)心；推進(jìn)劑貯箱為球形金屬膜片貯箱，包括2個(gè)氧化劑貯箱和2個(gè)燃燒劑貯箱，兩兩對(duì)稱布置。根據(jù)質(zhì)量特性測(cè)量結(jié)果，該航天器干重狀態(tài)下的質(zhì)量特性如表 2所示。

圖6 某充液航天器及其坐標(biāo)系定義Fig.6 Illustration for a liquid-filled spacecraft and its coordinate systems

表2 航天器干重狀態(tài)下的質(zhì)量特性

氧化劑密度為1445.1kg/m3，加注質(zhì)量為249.24 kg；燃燒劑密度為874.1 kg/m3，加注質(zhì)量為150.76 kg。建立推進(jìn)劑貯箱內(nèi)液體的有限元模型，如圖 7所示。根據(jù)表 2數(shù)據(jù)和本文方法可以計(jì)算出該航天器充液狀態(tài)下的質(zhì)量特性，如表 3所示。其中，金屬膜片貯箱加壓狀態(tài)對(duì)應(yīng)全充液工況，未加壓狀態(tài)對(duì)應(yīng)部分充液工況。與傳統(tǒng)的“固化液體”處理方法計(jì)算結(jié)果相比，由本文方法計(jì)算出來(lái)的轉(zhuǎn)動(dòng)慣量明顯減小，全充液工況最多減小4%，部分充液工況最多減小16%；而不同方法給出的質(zhì)量和質(zhì)心位置計(jì)算結(jié)果都是相同的。

圖7 充液航天器推進(jìn)劑貯箱內(nèi)液體的有限元模型Fig.7 Finite element model for the liquid in a propellant tank for a liquid-filled spacecraft

表3 航天器充液狀態(tài)下的質(zhì)量特性

6 結(jié)論

本文基于勢(shì)流理論推導(dǎo)了任意形狀貯箱在全充液和部分充液兩種狀態(tài)下液體質(zhì)量特性的一般表達(dá)式，建立了計(jì)算液體質(zhì)量特性的有限元數(shù)值方法；通過矩形貯箱液體質(zhì)量特性數(shù)值解與解析解的對(duì)比，驗(yàn)證了計(jì)算方法的正確性；對(duì)于布置多個(gè)推進(jìn)劑貯箱的充液航天器，分別采用“固化液體”方法和本文方法計(jì)算了它的質(zhì)量特性。研究表明，部分充液狀態(tài)下的轉(zhuǎn)動(dòng)慣量低于全充液狀態(tài)下的轉(zhuǎn)動(dòng)慣量，而兩者均低于固化液體狀態(tài)下的轉(zhuǎn)動(dòng)慣量。本文方法有助于提高充液航天器質(zhì)量特性計(jì)算的準(zhǔn)確性。除航天器外，所提出的方法亦可用于航空飛行器、油罐車、液貨船等其他充液系統(tǒng)質(zhì)量特性的計(jì)算和研究。