液體大幅晃動情形的航天器剛液耦合動響應仿真分析

岳寶增，馬伯樂，唐 勇，劉 峰

(1. 北京理工大學宇航學院，北京 100081； 2. 深圳大學電子與信息工程學院，深圳 518000；3.中國空間技術研究院，北京 100094)

0 引 言

現(xiàn)代航天器的液體燃料質量占比大，因此當其發(fā)生大幅晃動時，將影響航天器整體質量的分布。這將會導致航天器整體的慣性特性不斷劇烈變化，外加液體黏性引起的能量耗散等因素，最終可能造成充液航天器姿態(tài)失穩(wěn)。因此定量分析液體晃動與航天器的剛體運動相互耦合引發(fā)的復雜剛-液耦合動力學問題，是本文的重點研究方向。

剛-液耦合動力學的重難點是液體晃動動力學。特別在微重力環(huán)境下，液體的表面張力和接觸角等影響不可忽略。目前求得的解析解只限于如圓柱形和矩形等規(guī)則幾何形狀的貯箱。基于此學者們建立了航天器剛-液耦合動力學模型，并采用攝動法研究耦合系統(tǒng)的非線性特性或探究模塊化的求解方法等。但現(xiàn)代航天器通常采用球形或卡西尼(Cassini)貯箱，則難以采用解析方法處理液體晃動進而對剛-液耦合系統(tǒng)建模。目前等效力學模型(EMMs)和計算流體力學(CFD)方法在這方面取得了良好的效果。諸多學者基于等效球擺、質心約束面、運動脈動球等模型處理液體晃動問題，并建立了航天器的剛-液耦合動力學模型。從研究結果來看：質心約束面和運動脈動球模型，特別是液體發(fā)生大幅晃動時能更有效地預測航天器與液體的相互作用，這將是預測液體小幅晃動的彈簧質量阻尼器模型和球擺模型的有效補充。然而，基于等效力學模型建立的剛-液模型雖較為簡潔，但從精度上看還存在較大的缺陷，而且也需要通過CFD方法等進行相關參數(shù)的修正。CFD方法在無網(wǎng)格方法方面以SPH(Smoothed particle hydrodynamics)方法與MPS(Moving particle simulation)方法發(fā)展的最為成熟，應用最為廣泛。與網(wǎng)格類方法相比，無網(wǎng)格方法具有更大的靈活性，然而仍然存在一些需要解決的問題，比如不易施加本質邊界條件、計算量大、不易準確處理表面張力等。在網(wǎng)格方法方面，傳統(tǒng)的拉格朗日法在研究液體大幅晃動時會出現(xiàn)復雜曲面的網(wǎng)格畸變，因此很難被有效應用。歐拉方法中的VOF(Volume of Fluid)方法是目前商用軟件(如Flow3D，Ansys等)處理液體晃動采用最多的方法，在耦合航天器實際工程上實現(xiàn)了有效的應用，但在VOF方法中的流體體積分數(shù)函數(shù)是不連續(xù)的分段函數(shù)，因此無法模擬一些真實的液體晃動現(xiàn)象。Level-Set方法比VOF方法更容易重構自由面，但目前的研究多基于矩形貯箱，對于復雜貯箱的研究有待進一步加強。

ALE(Arbitrary Lagrange-Euler)方法旨在結合拉格朗日方法和歐拉方法各自的優(yōu)點，在研究液體晃動方面具有較大的發(fā)展空間。盡管ALE方法被廣泛研究，但多數(shù)只能解決規(guī)則幾何貯箱內(nèi)的液體晃動情況。如：Souli等和Alipour等分別采用ALE方法研究了橫向激勵下二維矩形貯箱內(nèi)的液體大幅晃動和三維矩形貯箱內(nèi)的液體大幅晃動。Okamoto等提出了一種新的分析多邊形水槽三維大幅晃動問題的ALE方法。本人編寫了一本關于非線性大振幅晃動動力學的綜合性書籍，詳細介紹了ALE有限元方法。為了解決復雜貯箱內(nèi)的液體晃動問題，設計合理有效的網(wǎng)格移動策略是根本措施。目前改進的ALE有限元方法有以下優(yōu)點:

1) 適用于球形貯箱內(nèi)三維大幅液體晃動(提出新的網(wǎng)格運動策略);

2) 適用于微重力環(huán)境下的大幅液體晃動(提出“無窮小接觸面元法”來施加接觸角邊界條件)。

本文基于改進的ALE有限元方法處理液體三維大幅晃動的問題。在耦合系統(tǒng)中，液體晃動力和晃動力矩作用于貯箱，影響貯箱的運動狀態(tài)；貯箱的運動同時反作用于液體，影響液體的運動。對于耦合系統(tǒng)中的液體部分，仍然考慮其相對貯箱的運動，貯箱運動對液體的影響以慣性力的形式來體現(xiàn)。液體在慣性參考系中的運動可以通過運動合成得到。貯箱在液體晃動力、晃動力矩以及外力、外力矩的作用下運動，不考慮貯箱的變形，并且貯箱與航天器平臺固連，被視為航天器主剛體。本文將給出剛-液耦合系統(tǒng)的動力學方程，并導出其具有數(shù)值穩(wěn)定性的形式，然后進行算例仿真與分析。

1 ALE方法基本理論

1.1 控制方程和邊界條件

對于黏性不可壓縮牛頓流體，ALE描述下的Navier-Stokes方程為可以寫成:

(1)

(2)

邊界條件如下，濕壁面包含無穿透邊界條件和Navier-滑移邊界條件，分別為：

·=0

(3)

··(-)=·(-)

(4)

自由液面可表示為：

-+··=-+2

(5)

式中：表示邊界處的單位外法向量;表示單位張量;為環(huán)境壓力;和分別為液體表面張力系數(shù)和平均曲率;稱為Navier滑移系數(shù)。

1.2 離散方法

為求解具有邊界條件的控制方程，采用特征線的時間分裂(CBS)算法進行時間離散，并采用標準伽遼金有限元法進行空間離散。同時，將壓力梯度投影的壓力穩(wěn)定方法應用于速度和壓力的求解，避免了Ladyzhenskaya-Brezzi-Babuska (LBB)條件的局限性。最后，將這三個步驟總結如下：

1)求解中間速度：

(6)

2)求解壓力:

(7)

3)速度修正，得到最終速度：

(8)

(9)

1.3 網(wǎng)格運動方法

球形貯箱的網(wǎng)格更新過程可分為接觸線、自由表面、濕容器壁和內(nèi)部流體域4個部分。更詳細的網(wǎng)格更新方法已在文獻[25]中提出。圖1～3分別為接觸線、自由表面和濕壁上節(jié)點的移動示意圖。

圖1 接觸線網(wǎng)格點(記為A1)的運動Fig.1 Movement of an arbitrary node on the contact line

圖2 自由液面網(wǎng)格點(記為A2)的運動Fig.2 Movement of an arbitrary node on the free surface

圖3 濕壁面網(wǎng)格點(記為A3)的運動Fig.3 Movement of an arbitrary node on the wet wall

自由表面和接觸線上的節(jié)點沿著規(guī)定的軌道運動,容器壁上的節(jié)點通過代數(shù)網(wǎng)格移動算法進行移動，節(jié)點內(nèi)部流體域通過求解拉普拉斯方程所得到。

1.4 接觸角邊界條件

圖4(a)顯示了部分原始自由曲面網(wǎng)格和接觸線。在該方法中，在接觸線與原始自由液面面單元之間添加。如圖4(b)所示，較暗的網(wǎng)格面是添加的無窮小接觸面元的一部分。

圖4 無窮小接觸面元法Fig.4 The infinitely small contact free-surface mesh method

更詳細的接觸角邊界條件已在文獻[26]中提出。將該方法嵌入ALE有限元求解框架中，并結合相應的網(wǎng)格運動算法，可以進行微重力環(huán)境下液體晃動的數(shù)值模擬。

2 航天器剛-液耦合動力學建模

2.1 坐標系

描述航天器剛-液耦合系統(tǒng)的各坐標系如圖5所示，′′′′為慣性參考系，為動參考系。將動參考系與主剛體固連，即參考系為剛體的隨體系，考慮隨體系下的液體運動，剛體運動對液體的影響以慣性力的形式來體現(xiàn)。

圖5 剛-液耦合系統(tǒng)示意圖Fig.5 Rigid-liquid coupling system

2.1 貯箱內(nèi)流體的晃動力與晃動力矩

液體對剛體的力為：

(10)

式中：表示液體區(qū)域,表示是貯箱內(nèi)流體對剛體的總力。

根據(jù)動量方程和散度定理，式(10)表達為：

(11)

液體對剛體產(chǎn)生的力矩為：

(12)

與力的推導類似，最后得到：

(13)

式中：是液體相對于點的慣性張量;表示流體質心相對于原點的距離。

2.2 考慮液體影響的剛體動力學方程

動參考系與剛體固連，所以剛體質點的相對速度與相對加速度恒為零，對剛體內(nèi)任一質點，有運動學關系：

(14)

由牛頓第二定律得：

(15)

將上式對剛體內(nèi)所有質點求和，可得：

(16)

剛體質點間內(nèi)力總是成對出現(xiàn)的，所以剛體內(nèi)力之和為零,再根據(jù)質心公式得到：

(17)

式中：為剛體質量；為剛體質心在參考系中的位置矢量；為剛體所受外力。

與力的推導類似，剛體的動量方程可以表示為：

(18)

式中：是剛體相對于點的慣性張量。

為了避免液體質量較大時的數(shù)值不穩(wěn)定，對剛-液耦合動力學方程進行變形，得到數(shù)值穩(wěn)定形式。最后，給出了系統(tǒng)的剛體動力學控制方程:

(19)

(20)

式中：為系統(tǒng)的總質量等于和之和；=(+)為系統(tǒng)質心的位置矢量；=+為系統(tǒng)相對于點的慣性張量;表示剛體所受外力矩;表示其它外力;表示貯箱內(nèi)流體對剛體的力矩;表示其它外力。

基于交錯法，建立了航天器剛-液耦合系統(tǒng)的數(shù)值模擬。該系統(tǒng)有液體模塊和剛體模塊兩個模塊。采用迭代法求解耦合模型，在時間步[,+1]中，對式(19)和(20)采用龍格-庫塔法進行數(shù)值求解。

3 算例分析與討論

3.1 充液圓柱貯箱耦合問題

首先驗證剛-液耦合方法建模的正確性。選取彈簧-阻尼作用下的充液圓柱貯箱剛-液耦合模型，進行數(shù)值仿真并與Peterson等的試驗結果進行對比。

如圖6所示，剛性圓柱貯箱與彈簧和阻尼器相連，貯箱僅能沿著軸方向做平動運動。貯箱質量=0.156 kg，貯箱半徑=0.0155 m，充液深度=2=0.031 m，液體密度=1000 kg/m，黏性系數(shù)=1×10Ns/m，表面張力系數(shù)=0.0357 N/m，接觸角=0°，重力環(huán)境=9.8 m/s，彈簧剛度系數(shù)=158.59 N/m，阻尼器阻尼系數(shù)=1.2 kg/s。初始時刻，貯箱具有位移，貯箱與液體均靜止，然后釋放貯箱，貯箱沿軸方向運動。

圖6 充液圓柱貯箱耦合系統(tǒng)Fig.6 Coupling system of liquid-filled cylindrical tanks

圖7所示為數(shù)值仿真與試驗下的貯箱位移對比結果。文獻[4]中，Peterson分別在設置初位移為3.9%、12.8%、17.0%三種工況下得到試驗結果。這里分別做出同等條件下的數(shù)值仿真，并將兩者的結果進行對比。結果表明，三種工況下的數(shù)值仿真結果均與試驗結果吻合良好，這充分表明了本文建立的剛-液耦合動力學模型具有很高的適用性。

圖7 貯箱位移量Fig.7 Tank displacement

為了研究液體晃動對耦合系統(tǒng)的影響，對考慮液體晃動和不考慮液體晃動的情況進行了數(shù)值仿真。不考慮液體晃動，即將液體凍結成剛體并與貯箱固連；考慮液體晃動即對液體作正常處理。圖8將兩種情況下貯箱位移的時間歷程進行了對比?？梢钥吹剑瑑煞N情況下貯箱位移的變化趨勢有著明顯的區(qū)別。不考慮液體晃動時，貯箱振幅持續(xù)衰減，考慮液體晃動時，第三個周期的振幅比前一周期有所增大，而且初始位移越大，振幅增大越明顯?？紤]液體晃動時，貯箱運動衰減得更慢，且振動周期更長，頻率更小。

圖8 貯箱位移量Fig.8 Tank displacement

圖9給出了數(shù)值仿真所得三種不同工況下貯箱位移的對比，并將位移做了歸一化處理。可以更明確地看到，當初始位移越大時，貯箱振動周期越大，頻率越小，表現(xiàn)出非線性軟彈簧的特性。這也和Peterson等在試驗中觀察到的現(xiàn)象一致。

圖9 三種不同x0工況下貯箱位移數(shù)值結果對比Fig.9 Comparison of numerical results of tank displacements under three different x0 conditions

圖10給出了初始位移12.8%工況下考慮晃動與不考慮晃動時晃動力的對比。不考慮晃動時，液體凍結為剛體并隨著貯箱一起運動，此時由于慣性效應，凍結的液體仍對貯箱產(chǎn)生作用力?？梢钥吹?，考慮液體晃動時，晃動力要大得多，其峰值大約是不考慮液體晃動時的兩倍，并且考慮液體晃動時晃動力存續(xù)時間更長。表明研究充液耦合系統(tǒng)問題時必須考慮液體晃動的影響。

圖10 x0=12.8%R時晃動力Fig.10 Slosh force (case: x0=12.8%R)

3.2 帶球形貯箱充液航天器的起旋機動

起旋機動是衛(wèi)星在軌運行時一種常見的機動形式。同時考慮球形貯箱和微重力環(huán)境。充液衛(wèi)星示意圖如圖11所示。隨體系原點位于貯箱球心。航天器的部分參數(shù)如下：

圖11 充液球形貯箱耦合系統(tǒng)Fig.11 Coupling system of liquid-filled spherical tanks

剛體部分:質量=24 kg，慣量主軸分別為R=58 kg·m、R=60 kg·m、R=5 kg·m，質心坐標=[0, 0, 1.5]m;

液體部分:貯箱半徑= 0.25 m,充液比40%,接觸角= 60°，黏性系數(shù)= 8.5×10Ns/m，表面張力系數(shù)= 0.0339 N/m，貯箱內(nèi)液體密度= 874.4 kg/m，容器的中心坐標為[0, 0, 0]m。

設定為零重力環(huán)境。數(shù)值仿真開始時，航天器處于穩(wěn)定自旋狀態(tài)，ω=[0, 0.1, 0]rad/s，且系統(tǒng)質心速度為零。在=2 s、=10 s、=30 s時，航天器分別受到持續(xù)時間為 1 s的階躍外力矩:

來增大航天器的角速度。在= 0 s時，指定慣性系與隨體系重合。

圖12所示為剛體角速度的時間歷程，并對比考慮液體晃動和不考慮液體晃動的結果。、、分別表示角速度矢量沿、、軸的分量。由于初始角速度和外力矩僅為沿軸的分量，因此和的值均為零。經(jīng)過了外力矩的三次加速，從0.1 rad/s增加到約0.14 rad/s。不考慮液體晃動時，呈現(xiàn)出平整的折線式變化，有外力矩時即加速，無外力矩時即恒定。考慮液體晃動時，在加速過程中會超過穩(wěn)態(tài)值，隨后回落，并經(jīng)歷振蕩之后，達到最終的穩(wěn)態(tài)值，在的變化過程中，可以清楚地看到液體晃動造成的影響。

圖12 剛體角速度Fig.12 Angular velocity of rigid body

圖13給出了貯箱壁面處的液體波高。圖13(a)此處給出的是沿軸方向的波高相對貯箱半徑的值?？梢钥吹?，前2 s內(nèi)液體沒有晃動，在施加外力矩之后，液體開始晃動。液體晃動時，由表面張力與離心力提供回復力。平面內(nèi)的液體波高較大，最大波高達到將近 60%。平面內(nèi)的液體波高較小，且基本上是同相位的，表明液體晃動基本上是關于平面對稱的，這與剛體角速度、外力矩的作用方式相一致。

圖13 貯箱壁面處液體波高Fig.13 Liquid wave-high at the tank wall

圖14給出了沿軸負方向觀察不同時刻的液體構形圖。初始時刻，在表面張力與離心力的共同作用下，自由液面既有下凹又有整體的弧形彎曲。液體晃動過程中，自由液面的形狀變化較為復雜，如=7.2 s時可以看到一階晃動模態(tài)(反對稱晃動，振型表現(xiàn)為1/2個周期正弦波形)與二階晃動模態(tài)(反對稱晃動，振型表現(xiàn)為3/2個周期正弦波形)的疊加，=48.5 s時可以看到三階晃動模態(tài)(對稱晃動，振型表現(xiàn)為1個周期正弦波形)的影響。=120 s時，液體穩(wěn)定到新的位置。

圖14 不同時刻的液體構形Fig.14 Liquid configuration at different times

圖15 晃動力和晃動力矩Fig.15 Slosh force and slosh torque

4 結 論

本文基于液體晃動的ALE有限元方法，對航天器剛-液耦合動力學進行了研究。推導了貯箱內(nèi)流體對剛體的作用力與力矩，以及受液體晃動影響時三維空間剛體動力學方程。為了避免液體質量較大時引起數(shù)值不穩(wěn)定，對剛-液耦合動力學方程加以變形，得到了數(shù)值穩(wěn)定的形式，并采用交錯方法對液體晃動模塊與剛體動力學模塊進行迭代求解，通過液體受慣性力、剛體受流體作用力與力矩來實現(xiàn)系統(tǒng)的耦合。最后進行了算例仿真與分析，研究了彈簧-阻尼作用下充液圓柱貯箱耦合問題，成功通過數(shù)值仿真觀察到耦合系統(tǒng)中貯箱橫向振動的非線性軟彈簧特性，數(shù)值結果與文獻中的試驗結果進行了對比，二者吻合良好。接下來研究了帶球形貯箱充液航天器的起旋機動，對比了考慮液體晃動與不考慮液體晃動兩種情況，明顯觀察到液體的高階模態(tài)以及液體晃動對耦合系統(tǒng)的影響。研究發(fā)現(xiàn)液體晃動對航天器的作用力與力矩存在復雜變化，其影響需要針對不同的情況作具體的定量分析。