充液航天器變質量液體大幅晃動的SPH分析方法

于 強，王天舒

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

0 引 言

一直以來，充液航天器貯箱內液體燃料及氧化劑的晃動作用對于航天器的姿軌控制都有著不可忽視的影響，甚至可能導致任務失敗，對液體晃動問題進行準確分析具有重要的意義。許多學者都對該問題進行了研究。

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

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

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

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

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

1 控制方程

1.1 SPH方法基本理論

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

(1)

(2)

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

(3)

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

(4)

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

(5)

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

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

(6)

Γi+fout,i

(7)

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

(8)

式中：ν是液體的運動黏性系數。

1.3 壓力計算方法

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

(9)

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

2 邊界處理方法

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

2.1 固壁邊界條件

基于邊界虛粒子的設置，本文采用的固壁邊界處理方法詳見文獻[10]。考慮無滑移邊界條件，給邊界虛粒子設置了與鏡像虛粒子邊界處理方法意義類似的鏡像速度：

(10)

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

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

(11)

邊界虛粒子的密度根據式(11)求出的壓力，利用式(9)的反函數計算得到。邊界虛粒子參與流體域實粒子的SPH計算，可以在沒有強沖擊的情況下有效防止流體實粒子的非物理穿透，這種處理方法對于充液航天器貯箱內液體晃動仿真非常適用。

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

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

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

(12)

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

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

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

vi=si0/Δtd

(13)

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

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

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

(14)

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

3 變質量液體晃動作用力及力矩的求解方法

在實際的應用當中，航天器的姿態信息變化是比較復雜的，動力學閉環仿真中的晃動仿真模塊需要通過航天器的姿態信息和充液比，反饋該時刻液體晃動對航天器產生的作用力和作用力矩。本文采用非慣性系下的SPH方法[16]，根據航天器姿態信息對每一個流體實粒子施加慣性力Finertial,i：

Finertial,i=-mi·

(15)

式中:V0是航天器質心的平動速度，dV0/dt為航天器質心的絕對加速度，ω和β分別為航天器質心本體坐標系的角速度和角加速度的矢量形式，ri為流體實粒子在航天器質心本體坐標系下的位置矢量，vi為流體實粒子相對于航天器質心本體坐標系的相對速度矢量。

至于液體晃動作用力及作用力矩的求解，本文采用了變質量質點系的動量定理及動量矩定理進行計算，不計算緩沖層流體實粒子的動量和動量矩，在每一個航天器動力學時間步Δtd中對累積流出貯箱進入緩沖層的粒子動量和角動量求和，分別得到ΔPout和ΔLout，則該時間步內液體晃動對于航天器產生的作用力Fslosh及作用力矩Mslosh分別為：

(16)

(17)

式中：N為該動力學時間步計算之后的貯箱內流體實粒子數目，Δ后加括號表示括號內物理量在該時間步的變化量，finertial,i和fout,i分別為第i個流體實粒子所受慣性力和重力導致的加速度。

SPH仿真液體晃動的計算模塊可以作為一個反饋模塊反映液體晃動產生的影響，與航天器系統動力學仿真形成一個閉環，從而實現兩者的聯合仿真。

4 算例分析和討論

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

4.1 常重環境球形貯箱短時間仿真

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

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

仿真時間為10 s，充液比從15%勻速減少到5%。考慮環境為常重環境，重力加速度g=9.8 m/s2，垂直初始水平液面向下(設為z軸負方向)，然后在水平x和y方向(與z方向構成了貯箱的隨體坐標系)分別施加強制的簡諧振動，振動位移形式分別如下所示：

x=0.15sinπt

(18)

y=0.1sin0.25πt

(19)

預平衡仿真時間2 s(實際仿真前僅施加重力進行的初始狀態計算)，動力學時間步長為0.01 s，SPH仿真時間步長Δt由CFL(Courant-Friedrichs-Lewy)穩定性條件[22]決定：

(20)

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

仿真實際用時52 s，效率上能夠達到嵌入航天器動力學仿真的要求。表1展示了部分時間點， SPH計算得到的液體粒子分布與Flow-3d計算結果的對比圖吻合程度較好，SPH方法可以比較好地刻畫自由液面的非線性現象。晃動作用力及力矩在貯箱隨體坐標系下投影的對比曲線如圖3和圖4所示。

表1 常重情況下球形貯箱內變質量液體晃動：SPH方法與Flow-3d液面對比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可知，此時SPH方法求解得到的晃動作用力及力矩曲線與Flow-3d仿真結果吻合良好，計算精度可以得到保證。相比于Flow-3d的計算結果，SPH求解得到的三軸晃動作用力峰值的相對偏差分別為0.94%，15.74%，9.72%；三軸晃動作用力矩峰值的相對偏差分別為8.64%，9.05%，2.12%。其中相對偏差最大的是水平y方向上的作用力，這主要是由于y方向上的激勵加速度峰值是事實上，包括SPH方法在內的粒子類方法會有比網格類方法更強的數值振蕩，即便采用數值耗散的方法，也會使計算得到的晃動作用力及力矩存在小幅高頻的非物理振蕩。但是，該振蕩沒有改變作用力變化趨勢，且不會累積沖量，在晃動幅度較大方向上衰減明顯，對整體計算結果不會產生顯著影響。總體而言，對于常重環境下的變質量液體晃動問題，本文提出的方法能夠得到合理的計算結果。

圖3 常重環境SPH方法與Flow-3d仿真球形貯箱內變質量液體晃動的作用力對比曲線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 常重環境SPH方法與Flow-3d仿真球形貯箱內變質量液體晃動的作用力矩對比曲線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，主要的晃動現象會呈現在x這一個方向上(見表1)。從曲線可以看出實際計算的y方向上晃動作用力量級大約是x方向上的1/10，其大小受離散粒子產生的數值振蕩影響相對較大。

4.2 低重環境卡西尼貯箱長時間仿真

卡西尼(Cassini)貯箱的中間段為圓柱，兩端為半球，在航天器中有廣泛的應用。本節針對卡西尼貯箱在低重環境下的復雜激勵情況，進行了長時間的變質量液體晃動仿真。

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

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

實際航天器在進行變軌或軟著陸時除了用推進的方式，還會使用動量輪調整姿態。算例采用了復雜的運動形式作為激勵輸入，航天器的質心平動加速度以及姿態變化規律在貯箱隨體坐標系下的投影如圖6和圖7所示(垂直于初始液面向上為z軸正方向)。

圖6 航天器質心的平動加速度Fig.6 Translational acceleration of the spacecraft mass center

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

SPH仿真實際用時2795 s，約為仿真時間的9倍，效率能夠滿足閉環仿真。計算得到的晃動作用力及力矩在貯箱隨體坐標系下的投影與Flow-3d計算結果對比曲線如圖8和圖9所示。

圖8 低重環境SPH方法與Flow-3d仿真卡西尼貯箱內變質量液體晃動的作用力對比曲線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 低重環境SPH方法與Flow-3d仿真卡西尼貯箱內變質量液體晃動的作用力矩的對比曲線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計算得到的結果與Flow-3d計算結果吻合較好；相比于Flow-3d計算結果，SPH方法求解得到的三軸晃動作用力峰值相對偏差分別為1.64%，10.50%，5.80%；三軸晃動作用力矩峰值的相對偏差分別為6.93%，2.12%，1.35%。

為了更清晰地區分兩種仿真結果，對圖8進行了局部放大，選取100～150 s區間內的結果如圖10所示。

圖10 低重環境SPH方法與Flow-3d仿真卡西尼貯箱內變質量液體晃動的作用力局部對比曲線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

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

綜上所述，本文提出的算法能夠解決低重環境復雜激勵下長時間變質量液體晃動的動力學問題，對充液航天器貯箱內變質量液體晃動的仿真分析十分有效。

5 結 論

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

最后通過與商業軟件Flow-3d計算結果的對比，該算法對于解決充液航天器貯箱內變質量液體大幅晃動問題的適用性和有效性得到了驗證。

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