于 強(qiáng)，王天舒

(清華大學(xué)航天航空學(xué)院，北京 100084)

0 引 言

一直以來，充液航天器貯箱內(nèi)液體燃料及氧化劑的晃動(dòng)作用對(duì)于航天器的姿軌控制都有著不可忽視的影響，甚至可能導(dǎo)致任務(wù)失敗，對(duì)液體晃動(dòng)問題進(jìn)行準(zhǔn)確分析具有重要的意義。許多學(xué)者都對(duì)該問題進(jìn)行了研究。

等效力學(xué)模型方法是目前廣泛應(yīng)用于航天領(lǐng)域液體小幅晃動(dòng)分析的研究方法[1-3]，目前也拓展到了液面非線性不太強(qiáng)的大幅晃動(dòng)問題中[4-5]。然而，當(dāng)航天器在進(jìn)行變軌或軟著陸等較大幅度機(jī)動(dòng)的時(shí)候，可能會(huì)激發(fā)貯箱內(nèi)液體的強(qiáng)非線性晃動(dòng)，甚至產(chǎn)生液面翻卷破碎的現(xiàn)象。此時(shí)等效力學(xué)模型不再適用，合理快速的數(shù)值仿真方法研究有其必要性。此外，航天器在進(jìn)行大幅度機(jī)動(dòng)的時(shí)候往往伴隨著燃料的消耗，貯箱內(nèi)的液體質(zhì)量在不斷減少，所以對(duì)液體晃動(dòng)的數(shù)值仿真方法還要考慮變質(zhì)量的因素。

基于網(wǎng)格的計(jì)算流體動(dòng)力學(xué)方法(CFD)是處理液體晃動(dòng)問題的一類常用方法，如任意拉格朗日-歐拉(ALE)方法[6]。但網(wǎng)格類方法處理自由液面非線性很強(qiáng)的情況，需要很精細(xì)的網(wǎng)格劃分才能比較好地刻畫自由液面，會(huì)降低計(jì)算效率，在實(shí)際工程應(yīng)用時(shí)有一定的局限性。光滑粒子流體動(dòng)力學(xué)(SPH)方法是一種基于拉格朗日描述的無網(wǎng)格粒子類方法，適合自由表面流動(dòng)問題的仿真[7]。其特點(diǎn)是將流體離散成一系列任意分布的粒子，并利用空間核函數(shù)近似以及粒子求和近似估算場(chǎng)函數(shù)及其空間導(dǎo)數(shù)。粒子之間不存在物理上的連接，在計(jì)算大變形自由液面問題時(shí)存在優(yōu)勢(shì)。

變質(zhì)量晃動(dòng)問題反映在SPH方法中是流出流量可控的開口邊界問題。邊界條件處理一直都是SPH方法的重要研究方向，大致可分為固壁和開口兩類。目前為止，SPH處理固壁邊界的方法主要可以分為兩大類，第一類是設(shè)置虛粒子的方法，主要有三種實(shí)現(xiàn)方式：一是邊界虛粒子對(duì)流體粒子直接建立斥力模型[8]，阻止粒子穿透邊界；二是認(rèn)為相對(duì)靜止的邊界虛粒子具有流體粒子的性質(zhì)[9-10]，通過計(jì)算合理的邊界虛粒子壓力阻止流體粒子流出；三是通過設(shè)置與流體粒子相對(duì)于邊界對(duì)稱的鏡像虛粒子，使得流體粒子在固壁邊界處形成近似反彈的效果以達(dá)到固壁邊界條件的實(shí)現(xiàn)[11]。第二類方法是一致化半解析邊界條件(Unified semi-analytical wall boundary conditions, USAW)處理方法[12]。在該方法中，邊界是按照網(wǎng)格進(jìn)行劃分，核函數(shù)近似解在邊界被截?cái)嗟牟糠蛛x散到邊界網(wǎng)格節(jié)點(diǎn)上，進(jìn)而推導(dǎo)出邊界上的控制方程。相比于虛粒子設(shè)置方法，USAW方法具備一定的數(shù)學(xué)理論基礎(chǔ)，可以比較準(zhǔn)確地反映固壁邊界附近的流動(dòng)，并描述任意紐曼(Neumann)或狄利克雷(Dirichlet)邊界條件，但該方法在計(jì)算效率方面相對(duì)于虛粒子方法有明顯不足。Mayrhofer等[13]將USAW方法應(yīng)用到了三維情況，拓展了其應(yīng)用范圍。駱釗等[14]在USAW方法基礎(chǔ)上提出了無質(zhì)量邊界粒子概念，避免了邊界零粒子層現(xiàn)象。

對(duì)應(yīng)固壁邊界條件，SPH方法對(duì)于開口邊界條件的處理方法也主要分為虛粒子類方法和半解析類方法兩類。前者主要采用設(shè)置緩沖層并填滿虛粒子的處理方式[15]，借鑒網(wǎng)格類處理開口邊界的思想，為緩沖層虛粒子設(shè)置滿足邊界物理?xiàng)l件的物理量，緩沖區(qū)粒子參與SPH計(jì)算。該方法實(shí)現(xiàn)簡(jiǎn)單，但由于粒子的離散特性，難以準(zhǔn)確描述邊界處的流量。后者則利用USAW方法的處理思想控制開口邊界處的流量以實(shí)現(xiàn)開口邊界條件。Leroy等[16]結(jié)合ISPH方法和USAW條件，實(shí)現(xiàn)了二維及三維情況下定流速開口邊界工況的計(jì)算，并得到了很好的效果。但是，半解析的開口邊界處理方法在實(shí)際應(yīng)用時(shí)，尤其對(duì)于三維情況，龐大的計(jì)算量限制了其在實(shí)際工程領(lǐng)域上的應(yīng)用。

針對(duì)充液航天器燃料消耗過程中大機(jī)動(dòng)可能引起的液體晃動(dòng)問題，本文出于對(duì)計(jì)算效率以及與航天器動(dòng)力學(xué)程序聯(lián)合仿真的需求，在非慣性系下SPH(Non-inertial SPH, NI-SPH)方法[17]的基礎(chǔ)上，對(duì)虛粒子類開口邊界處理方法進(jìn)行了改進(jìn)，使其能夠準(zhǔn)確控制出口邊界處的流量，從而實(shí)現(xiàn)變質(zhì)量液體晃動(dòng)的SPH仿真。此外，本文還利用變質(zhì)量質(zhì)點(diǎn)系的動(dòng)量和動(dòng)量矩定理計(jì)算了液體晃動(dòng)的作用力及作用力矩，與航天器系統(tǒng)動(dòng)力學(xué)仿真程序形成了閉環(huán)。最后，通過將仿真結(jié)果與CFD商業(yè)軟件Flow-3d結(jié)果進(jìn)行對(duì)比，驗(yàn)證了算法的有效性。

1 控制方程

1.1 SPH方法基本理論

標(biāo)準(zhǔn)SPH方法是將液體部分離散成一系列具有質(zhì)量的拉格朗日粒子，不存在碰撞體積，然后利用核函數(shù)近似和粒子求和近似的方法用粒子所攜帶的物理量對(duì)宏觀物理量f(xi)(如密度、壓強(qiáng)等)及其空間導(dǎo)數(shù)進(jìn)行近似，形式如下：

(1)

(2)

其中，下標(biāo)代表粒子編號(hào)，Wij=W(xi-xj,h)為核函數(shù)；x，ρ，m分別表示相應(yīng)粒子的位置向量、密度和質(zhì)量；h表示核函數(shù)的光滑長(zhǎng)度，根據(jù)核函數(shù)的不同，其作用范圍，即支持域的半徑為κh(κ取2或3)，h一般可以取粒子初始間距r0的倍數(shù)，即h=kr0；N為粒子i支持域內(nèi)的粒子數(shù)目；此外，有

(3)

本文的核函數(shù)采用的是目前最為廣泛的B-樣條函數(shù)(取κ=2,k=1.3)，如下所示：

(4)

其中,R=r/h，r為粒子間距；對(duì)于一維、二維和三維問題分別有

(5)

1.2 N-S方程的SPH離散形式

利用SPH近似方法對(duì)拉格朗日描述下Navier-Stokes方程進(jìn)行離散，可以得到N-S方程的SPH離散形式[18]：

(6)

Γi+fout,i

(7)

其中,v和fout分別代表速度矢量和重力產(chǎn)生的加速度，物理量的下標(biāo)ij代表粒子i和j的相應(yīng)物理量之差。Γi為簡(jiǎn)化的物理黏性耗散項(xiàng)，采用由Lo和Shao[19]提出的公式：

(8)

式中：ν是液體的運(yùn)動(dòng)黏性系數(shù)。

1.3 壓力計(jì)算方法

SPH方法在計(jì)算不可壓縮流體問題時(shí)，通過引進(jìn)人工壓縮率將不可壓縮流近似看作弱可壓流體。對(duì)于一般液體的晃動(dòng)問題，本文采用下面的狀態(tài)方程由密度變化顯式計(jì)算壓力[20]：

(9)

式中:ρ0是液體初始密度；γ為常數(shù)，一般γ=7；pb為背景壓力，一般可置為0，式(9)計(jì)算的便是流體的相對(duì)壓力場(chǎng)；cs為人工聲速，本文取cs=15。γ和cs的選取決定了液體的剛性，限制液體粒子的密度變化率在1%之內(nèi)，能夠保證流體的弱可壓性質(zhì)。

2 邊界處理方法

由于SPH離散粒子的拉格朗日性質(zhì)，SPH方法在表面張力作用不明顯的情況下可以自然滿足自由邊界條件，而對(duì)于固壁邊界條件和流量可控的出口邊界條件，本文均基于虛粒子類方法進(jìn)行了實(shí)現(xiàn)。

2.1 固壁邊界條件

基于邊界虛粒子的設(shè)置，本文采用的固壁邊界處理方法詳見文獻(xiàn)[10]?？紤]無滑移邊界條件，給邊界虛粒子設(shè)置了與鏡像虛粒子邊界處理方法意義類似的鏡像速度：

(10)

其中，w代表邊界虛粒子。該速度用于方程(6)和(8)的計(jì)算，但不用于虛粒子相對(duì)位置更新，邊界虛粒子在非慣性系下的位置始終保持不變。

邊界虛粒子的壓力則需要較好地描述固壁邊界附近的液體壓力梯度，綜合考慮非慣性系下的慣性力和重力，邊界虛粒子的壓力由SPH方法離散得到：

(11)

邊界虛粒子的密度根據(jù)式(11)求出的壓力，利用式(9)的反函數(shù)計(jì)算得到。邊界虛粒子參與流體域?qū)嵙Ｗ拥腟PH計(jì)算，可以在沒有強(qiáng)沖擊的情況下有效防止流體實(shí)粒子的非物理穿透，這種處理方法對(duì)于充液航天器貯箱內(nèi)液體晃動(dòng)仿真非常適用。

2.2 流量可控的出口邊界條件

由于SPH粒子的離散性質(zhì)，SPH方法中的物理量分布是離散的，不是連續(xù)變化的，難以準(zhǔn)確刻畫一個(gè)截面上物理量的變化率，因此基于虛粒子設(shè)置的緩沖層開口邊界處理方法在處理控制邊界流量的開口邊界問題有著天然的劣勢(shì)，無法像網(wǎng)格類方法簡(jiǎn)單地在邊界結(jié)點(diǎn)處設(shè)置固定速度。

針對(duì)這一問題，本文將貯箱及緩沖層區(qū)域(可認(rèn)為是貯箱出口管)統(tǒng)一用固壁邊界虛粒子封閉，正常情況下均按照SPH控制方程進(jìn)行計(jì)算。在每一個(gè)航天器動(dòng)力學(xué)時(shí)間步Δtd(遠(yuǎn)大于SPH仿真時(shí)間步)中，根據(jù)所給的當(dāng)前時(shí)刻充液比σ，計(jì)算得到該動(dòng)力學(xué)時(shí)間步要流出流體區(qū)域的粒子數(shù)目：

(12)

其中，N和N0分別是上一時(shí)刻剩余的流體粒子數(shù)和初始流體粒子數(shù)；σ0是初始的充液比；方括號(hào)代表高斯取整函數(shù)。根據(jù)需要流出粒子的數(shù)目ΔN，在緩沖層流出邊界的鄰近網(wǎng)格(鏈表搜索法搜索鄰近粒子對(duì)設(shè)置的網(wǎng)格)中搜索最接近邊界的ΔN個(gè)粒子設(shè)置為邊緣流出粒子，如圖1所示，記錄此時(shí)這些粒子與邊界的垂直距離分別為si0,i=1,2,…,ΔN。這些粒子具有如下性質(zhì)：

圖1 流出邊界附近的粒子分布示意圖Fig.1 Particles distribution near the outlet boundary

1)解除邊界虛粒子對(duì)邊緣流出粒子的影響，該類粒子具備垂直于流出邊界向外的速度，其大小為

vi=si0/Δtd

(13)

2)邊緣流出粒子的壓力和密度分別按照固壁邊界條件的方法進(jìn)行計(jì)算。

3)邊緣流出粒子對(duì)于其他流體粒子的作用均乘以一個(gè)衰減系數(shù)ξi，該系數(shù)由粒子與流出邊界的實(shí)時(shí)距離si計(jì)算得到，其表達(dá)式為：

ξi=(esi/si0-1)/(e-1)

(14)

第一條性質(zhì)使得邊緣流出粒子能夠在該動(dòng)力學(xué)時(shí)間步的計(jì)算后運(yùn)動(dòng)到邊界處，進(jìn)而被設(shè)置成場(chǎng)外粒子，不再參與計(jì)算。而后兩條性質(zhì)則是在保持邊緣流出粒子對(duì)于緩沖層粒子作用的同時(shí)，保證邊界虛粒子對(duì)緩沖層其他流體實(shí)粒子的阻擋作用依然合理。該方法可以有效控制流出邊界的流量，并讓出口邊界附近的流場(chǎng)比較合理，減少數(shù)值震蕩。這種離散的邊界處理方法是基于SPH粒子進(jìn)行人工操作的，數(shù)學(xué)上不能嚴(yán)格對(duì)應(yīng)流體在邊界處要滿足的方程，但由于緩沖層的存在，這種人工處理方法對(duì)于貯箱內(nèi)流體實(shí)粒子的擾動(dòng)很小，可以適用于充液航天器內(nèi)的變質(zhì)量液體晃動(dòng)問題。此外，方法中衰減系數(shù)的表達(dá)式也可以采用其他類型函數(shù)，滿足初始為1，粒子到達(dá)邊界時(shí)為0即可。

3 變質(zhì)量液體晃動(dòng)作用力及力矩的求解方法

在實(shí)際的應(yīng)用當(dāng)中，航天器的姿態(tài)信息變化是比較復(fù)雜的，動(dòng)力學(xué)閉環(huán)仿真中的晃動(dòng)仿真模塊需要通過航天器的姿態(tài)信息和充液比，反饋該時(shí)刻液體晃動(dòng)對(duì)航天器產(chǎn)生的作用力和作用力矩。本文采用非慣性系下的SPH方法[16]，根據(jù)航天器姿態(tài)信息對(duì)每一個(gè)流體實(shí)粒子施加慣性力Finertial,i：

Finertial,i=-mi·

(15)

式中:V0是航天器質(zhì)心的平動(dòng)速度，dV0/dt為航天器質(zhì)心的絕對(duì)加速度，ω和β分別為航天器質(zhì)心本體坐標(biāo)系的角速度和角加速度的矢量形式，ri為流體實(shí)粒子在航天器質(zhì)心本體坐標(biāo)系下的位置矢量，vi為流體實(shí)粒子相對(duì)于航天器質(zhì)心本體坐標(biāo)系的相對(duì)速度矢量。

至于液體晃動(dòng)作用力及作用力矩的求解，本文采用了變質(zhì)量質(zhì)點(diǎn)系的動(dòng)量定理及動(dòng)量矩定理進(jìn)行計(jì)算，不計(jì)算緩沖層流體實(shí)粒子的動(dòng)量和動(dòng)量矩，在每一個(gè)航天器動(dòng)力學(xué)時(shí)間步Δtd中對(duì)累積流出貯箱進(jìn)入緩沖層的粒子動(dòng)量和角動(dòng)量求和，分別得到ΔPout和ΔLout，則該時(shí)間步內(nèi)液體晃動(dòng)對(duì)于航天器產(chǎn)生的作用力Fslosh及作用力矩Mslosh分別為：

(16)

(17)

式中：N為該動(dòng)力學(xué)時(shí)間步計(jì)算之后的貯箱內(nèi)流體實(shí)粒子數(shù)目，Δ后加括號(hào)表示括號(hào)內(nèi)物理量在該時(shí)間步的變化量，finertial,i和fout,i分別為第i個(gè)流體實(shí)粒子所受慣性力和重力導(dǎo)致的加速度。

SPH仿真液體晃動(dòng)的計(jì)算模塊可以作為一個(gè)反饋模塊反映液體晃動(dòng)產(chǎn)生的影響，與航天器系統(tǒng)動(dòng)力學(xué)仿真形成一個(gè)閉環(huán)，從而實(shí)現(xiàn)兩者的聯(lián)合仿真。

4 算例分析和討論

本文所有算例測(cè)試均在一臺(tái)配置4核8線程4.20 GHz主頻的Intel i7處理器和16 GB內(nèi)存的個(gè)人電腦上完成。根據(jù)文獻(xiàn)[21]中的驗(yàn)證，非慣性系下SPH方法可以在粒子初始間距較大，即粒子數(shù)較少的情況下得到合理的結(jié)果。為了體現(xiàn)算法的效率，本文算例采取了粒子初始間距較大的情況(將粒子初始間距選為貯箱半徑的1/10左右)。此外，算例以CFD商業(yè)軟件Flow-3d的計(jì)算結(jié)果作為參照，該軟件采用了VOF方法計(jì)算自由表面流，在仿真貯箱內(nèi)液體晃動(dòng)問題上有穩(wěn)定良好的表現(xiàn)。

4.1 常重環(huán)境球形貯箱短時(shí)間仿真

球形貯箱的半徑為0.62 m，貯箱中心位置為(0.62，0.62，0.62)m，液體為水，其密度為1000 kg/m3，動(dòng)力黏性系數(shù)為9.29×10-4kg/(m·s)。粒子初始間距為0.05 m，初始充液比為15%，設(shè)置了1322個(gè)實(shí)粒子(包括緩沖區(qū))和2284個(gè)邊界虛粒子。初始狀態(tài)下的液體分布及相應(yīng)SPH粒子離散如圖2所示。

圖2 球形貯箱初始液體及SPH粒子分布圖Fig.2 Initial liquid and SPH particles distribution in the sphere tank

仿真時(shí)間為10 s，充液比從15%勻速減少到5%?？紤]環(huán)境為常重環(huán)境，重力加速度g=9.8 m/s2，垂直初始水平液面向下(設(shè)為z軸負(fù)方向)，然后在水平x和y方向(與z方向構(gòu)成了貯箱的隨體坐標(biāo)系)分別施加強(qiáng)制的簡(jiǎn)諧振動(dòng)，振動(dòng)位移形式分別如下所示：

x=0.15sinπt

(18)

y=0.1sin0.25πt

(19)

預(yù)平衡仿真時(shí)間2 s(實(shí)際仿真前僅施加重力進(jìn)行的初始狀態(tài)計(jì)算)，動(dòng)力學(xué)時(shí)間步長(zhǎng)為0.01 s，SPH仿真時(shí)間步長(zhǎng)Δt由CFL(Courant-Friedrichs-Lewy)穩(wěn)定性條件[22]決定：

(20)

式中:h為光滑長(zhǎng)度，cs為人工聲速,vi為粒子i的速度大小。

仿真實(shí)際用時(shí)52 s，效率上能夠達(dá)到嵌入航天器動(dòng)力學(xué)仿真的要求。表1展示了部分時(shí)間點(diǎn)， SPH計(jì)算得到的液體粒子分布與Flow-3d計(jì)算結(jié)果的對(duì)比圖吻合程度較好，SPH方法可以比較好地刻畫自由液面的非線性現(xiàn)象?；蝿?dòng)作用力及力矩在貯箱隨體坐標(biāo)系下投影的對(duì)比曲線如圖3和圖4所示。

表1 常重情況下球形貯箱內(nèi)變質(zhì)量液體晃動(dòng)：SPH方法與Flow-3d液面對(duì)比Table 1 Variable-mass liquid sloshing in the sphere tank under normal gravity: SPH results of flow patterns in comparison with those from Flow-3d

由圖3～圖4可知，此時(shí)SPH方法求解得到的晃動(dòng)作用力及力矩曲線與Flow-3d仿真結(jié)果吻合良好，計(jì)算精度可以得到保證。相比于Flow-3d的計(jì)算結(jié)果，SPH求解得到的三軸晃動(dòng)作用力峰值的相對(duì)偏差分別為0.94%，15.74%，9.72%；三軸晃動(dòng)作用力矩峰值的相對(duì)偏差分別為8.64%，9.05%，2.12%。其中相對(duì)偏差最大的是水平y(tǒng)方向上的作用力，這主要是由于y方向上的激勵(lì)加速度峰值是事實(shí)上，包括SPH方法在內(nèi)的粒子類方法會(huì)有比網(wǎng)格類方法更強(qiáng)的數(shù)值振蕩，即便采用數(shù)值耗散的方法，也會(huì)使計(jì)算得到的晃動(dòng)作用力及力矩存在小幅高頻的非物理振蕩。但是，該振蕩沒有改變作用力變化趨勢(shì)，且不會(huì)累積沖量，在晃動(dòng)幅度較大方向上衰減明顯，對(duì)整體計(jì)算結(jié)果不會(huì)產(chǎn)生顯著影響?？傮w而言，對(duì)于常重環(huán)境下的變質(zhì)量液體晃動(dòng)問題，本文提出的方法能夠得到合理的計(jì)算結(jié)果。

圖3 常重環(huán)境SPH方法與Flow-3d仿真球形貯箱內(nèi)變質(zhì)量液體晃動(dòng)的作用力對(duì)比曲線Fig.3 Comparison of variable-mass liquid sloshing force histories predicted by the SPH and those simulated by Flow-3d in the sphere tank under normal gravity

圖4 常重環(huán)境SPH方法與Flow-3d仿真球形貯箱內(nèi)變質(zhì)量液體晃動(dòng)的作用力矩對(duì)比曲線Fig.4 Comparison of variable-mass liquid sloshing moment histories predicted by the SPH and those simulated by Flow-3d in the sphere tank under normal gravity

x方向上的1/24，主要的晃動(dòng)現(xiàn)象會(huì)呈現(xiàn)在x這一個(gè)方向上(見表1)。從曲線可以看出實(shí)際計(jì)算的y方向上晃動(dòng)作用力量級(jí)大約是x方向上的1/10，其大小受離散粒子產(chǎn)生的數(shù)值振蕩影響相對(duì)較大。

4.2 低重環(huán)境卡西尼貯箱長(zhǎng)時(shí)間仿真

卡西尼(Cassini)貯箱的中間段為圓柱，兩端為半球，在航天器中有廣泛的應(yīng)用。本節(jié)針對(duì)卡西尼貯箱在低重環(huán)境下的復(fù)雜激勵(lì)情況，進(jìn)行了長(zhǎng)時(shí)間的變質(zhì)量液體晃動(dòng)仿真。

設(shè)置卡西尼貯箱的半球及圓柱段半徑為0.492 m，圓柱段的高度為0.1 m，貯箱幾何中心在航天器質(zhì)心本體坐標(biāo)系中的位置坐標(biāo)為(0.0，-0.875，-0.681)m。液體為航天器中常用的燃燒劑甲基肼(MMH)，密度為874.4 kg/m3，動(dòng)力黏性系數(shù)為8.50×10-4kg/(m·s)。初始充液比為60%，仿真時(shí)間為300 s，勻速減少到35%。粒子初始間距為0.045 m，初始共設(shè)置了3324個(gè)流體實(shí)粒子(含緩沖層)和1982個(gè)虛粒子，初始狀態(tài)下的液體分布及相應(yīng)SPH粒子離散如圖5所示。

圖5 卡西尼貯箱初始液體及粒子分布圖Fig.5 Initial liquid and particles distribution in the Cassini tank

實(shí)際航天器在進(jìn)行變軌或軟著陸時(shí)除了用推進(jìn)的方式，還會(huì)使用動(dòng)量輪調(diào)整姿態(tài)。算例采用了復(fù)雜的運(yùn)動(dòng)形式作為激勵(lì)輸入，航天器的質(zhì)心平動(dòng)加速度以及姿態(tài)變化規(guī)律在貯箱隨體坐標(biāo)系下的投影如圖6和圖7所示(垂直于初始液面向上為z軸正方向)。

圖6 航天器質(zhì)心的平動(dòng)加速度Fig.6 Translational acceleration of the spacecraft mass center

圖7 航天器角速度Fig.7 Angular velocity of the spacecraft

SPH仿真實(shí)際用時(shí)2795 s，約為仿真時(shí)間的9倍，效率能夠滿足閉環(huán)仿真。計(jì)算得到的晃動(dòng)作用力及力矩在貯箱隨體坐標(biāo)系下的投影與Flow-3d計(jì)算結(jié)果對(duì)比曲線如圖8和圖9所示。

圖8 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動(dòng)的作用力對(duì)比曲線Fig.8 Comparison of variable-mass liquid sloshing force histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

圖9 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動(dòng)的作用力矩的對(duì)比曲線Fig.9 Comparison of variable-mass liquid sloshing moment histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

從圖8～9可以看出，整體上SPH計(jì)算得到的結(jié)果與Flow-3d計(jì)算結(jié)果吻合較好；相比于Flow-3d計(jì)算結(jié)果，SPH方法求解得到的三軸晃動(dòng)作用力峰值相對(duì)偏差分別為1.64%，10.50%，5.80%；三軸晃動(dòng)作用力矩峰值的相對(duì)偏差分別為6.93%，2.12%，1.35%。

為了更清晰地區(qū)分兩種仿真結(jié)果，對(duì)圖8進(jìn)行了局部放大，選取100～150 s區(qū)間內(nèi)的結(jié)果如圖10所示。

圖10 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動(dòng)的作用力局部對(duì)比曲線Fig.10 Comparison of variable-mass liquid sloshing local force histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

兩種仿真結(jié)果變化趨勢(shì)吻合良好，但通過進(jìn)一步觀察，可以發(fā)現(xiàn)隨著時(shí)間的推移，兩種仿真結(jié)果差異有明顯的增大趨勢(shì)，尤其是z方向上的晃動(dòng)作用力，在仿真最后已經(jīng)可以觀察出有明顯的偏移。經(jīng)過檢查，發(fā)現(xiàn)Flow-3d在計(jì)算長(zhǎng)時(shí)間自由表面流動(dòng)的時(shí)候，由于自由液面處網(wǎng)格之間的流量計(jì)算存在誤差，會(huì)出現(xiàn)比較明顯的體積計(jì)算誤差。針對(duì)本算例，F(xiàn)low-3d中液體體積在300 s的仿真過程中大約增加了7%，后期的實(shí)際液體體積大于既定值，最后的晃動(dòng)作用力及力矩也因此相對(duì)偏大。當(dāng)仿真時(shí)間進(jìn)一步加長(zhǎng)時(shí)，F(xiàn)low-3d的體積累積誤差會(huì)進(jìn)一步增大。此外，根據(jù)實(shí)際測(cè)試，當(dāng)自由液面面積更大，變化更加劇烈的時(shí)候，液體體積的計(jì)算誤差累積速度也會(huì)變得更快。因此，考慮到SPH粒子系的質(zhì)量守恒性，SPH方法在長(zhǎng)時(shí)間仿真時(shí)具有更好的穩(wěn)定性。

綜上所述，本文提出的算法能夠解決低重環(huán)境復(fù)雜激勵(lì)下長(zhǎng)時(shí)間變質(zhì)量液體晃動(dòng)的動(dòng)力學(xué)問題，對(duì)充液航天器貯箱內(nèi)變質(zhì)量液體晃動(dòng)的仿真分析十分有效。

5 結(jié) 論

本文基于非慣性系下SPH方法的基本理論，針對(duì)存在燃料消耗的充液航天器大幅機(jī)動(dòng)情況，提出了一種可以控制流量的虛粒子類出口邊界處理方法。該方法能夠根據(jù)輸入的航天器姿態(tài)信息和充液比，完成單動(dòng)力學(xué)時(shí)間步的SPH更新，并采用變質(zhì)量質(zhì)點(diǎn)系動(dòng)量定理和動(dòng)量矩定理計(jì)算并輸出晃動(dòng)作用力及力矩，與航天器系統(tǒng)動(dòng)力學(xué)仿真形成閉環(huán)。

最后通過與商業(yè)軟件Flow-3d計(jì)算結(jié)果的對(duì)比，該算法對(duì)于解決充液航天器貯箱內(nèi)變質(zhì)量液體大幅晃動(dòng)問題的適用性和有效性得到了驗(yàn)證。

SPH方法能夠較好地解決常重及低重環(huán)境下的液體晃動(dòng)問題，但對(duì)于微重環(huán)境表面張力作為主要驅(qū)動(dòng)力的情況還未有簡(jiǎn)潔有效的處理方法，這也是今后需要研究的重要方向。