基于SPH方法的球形貯箱液體晃動動力學分析

馬駿驍1，馬亮，魏承，*，胡月，趙陽

1. 中國空間技術研究院 通信衛星事業部，北京 100094 2. 哈爾濱工業大學 航天工程與力學系，哈爾濱 150001

隨著航天器的持續飛行，貯箱內的燃料也在不斷消耗，理想情況為貯箱內燃料含量由充滿狀態至耗盡狀態。不同的燃料含量意味著液體質量的不同，從而導致液體由于晃動對貯箱產生的力和壓強的不同。貯箱內液體燃料的含量總是處于變化之中，不同含量的液體燃料產生的晃動行為會具有怎樣的差別，需要進行對比分析。光滑粒子流體動力學方法(SPH方法)自1977年由Lucy[1]和Gingold[2]首次提出。基于SPH方法的流體動力學問題研究經歷了逐步完善的過程：Monaghan[3]于1994年首次將SPH方法引入流體動力學問題的求解，隨后多位學者基于此方法對鈍體繞流[4]、剪切流動[5]、潰壩[6]等問題進行了仿真分析。同時，SPH方法在計算精度[7]、計算穩定性[8]及邊界處理[9]等方面經過不斷的改進和修正，使其成為一種解決流體動力學問題的成熟建模方法。

液體晃動問題主要具有三類研究方案：理論研究、數值研究、試驗研究。理論研究中，最為常用的方法為等效模型方法，但等效模型不足以準確地捕捉自由表面的形狀，使其無法描述晃動液體的實際幾何形狀。數值研究以計算流體動力學(CFD)為主，但無論哪種CFD方法，都無法對大幅液體晃動進行準確描述，尤其當液體發生破碎等強非線性現象。試驗研究多用于驗證建模方法的準確性或對比分析防晃結構的優劣，但試驗方法成本太高，且受多種試驗因素制約，對部分試驗條件難以模擬，因此無法大規模應用。由此對于液體晃動問題，無網格的SPH方法得到了快速的發展與應用。劉富[10]基于SPH方法對飛機油箱進行晃動分析，并提出晃動抑制方案。朱小松[11]則將并行計算引入SPH方法，試圖提升仿真效率。張凱凱[12]建立液體晃動試驗平臺以驗證SPH數值模型的正確性，在此基礎上對機翼油箱的防晃結構進行優化設計。章恒[13]針對流固耦合問題，采用SPH方法通過改進流動模擬算法及耦合處理算法，對常見場景進行模擬仿真，驗證了算法的正確性。杜林霏[14]通過改進傳統SPH方法中壓強、邊界等方面的處理方法，同時給出SPH模型的晃動力、力矩計算方法，對二維半圓形貯箱的晃動問題進行仿真分析。聶霄[15]通過并行計算理論分析了不可壓縮流(ISPH)的晃動現象。李清[16]則分別采用了理論方法、SPH方法、試驗方法對貯箱晃動穩定性進行了分析。上述分析均缺乏對三維球形貯箱的仿真研究，而這正是本文的重點研究內容。

本文針對航天器球形貯箱內的液體晃動問題，采用SPH方法進行建模仿真，建立不同仿真需求的航天器球形貯箱模型。通過記錄貯箱受力情況及設置測量點壓強情況，對不同仿真參數進行影響分析，并最終得到相應規律，從而為航天器球形貯箱的設計提供技術依據。

1 基于SPH方法的流體動力學建模

1.1 SPH方法近似原理

SPH方法的核心近似原理包括核近似、粒子近似，如圖1所示，圖中k和h表示光滑因子和光滑長度。核近似的目標是通過構造適合的光滑函數(核函數)，近似等效替代狄拉克函數，從而將任意函數表示為積分形式。粒子近似的過程就是將積分計算離散為求和計算的過程，光滑半徑控制求和范圍。通過上述兩步近似，任意函數及其導數即可在SPH理論框架內由結構更簡單、求導更方便的核函數進行計算求解。

圖1 SPH方法核心近似原理Fig.1 The core principle of approximate SPH method

(1)核近似

核近似是SPH方法近似過程的第一步，利用狄拉克函數的性質，將任意函數表示為積分形式。設任意函數f(x)，x為空間位置向量，該函數可表示為一種含有狄拉克函數的積分形式：

(1)

式中：x′為定義域內某一空間位置向量，Ω為體積分域，δ(x-x′)為狄拉克函數，δ函數定義為：

(2)

由于δ函數“雖連續但不可導”，構造一種光滑函數(核函數)W(x-x′,h)替代δ函數，如圖1(a)所示。將核函數代替式(1)中的δ函數，則任意函數f(x)表示為：

(3)

式(3)中方程為約等，這是由于核函數不能完全等價于狄拉克函數，可以通過泰勒展開證明采用核函數替代的任意函數積分表達式具有二階精度。

·[f(x′)W(x-x′,h)]dx′-

(4)

采用散度定理，對式(4)中第一項進行積分變換，將體積分變為面積分，可以得到：

(5)

式中:n是面域S的單位法向量。當核函數支持域處于問題域內部時，面積分項為零,可以得到：

(6)

當核函數支持域與問題域存在重疊時，面積分項不為零，由此引發SPH方法的邊界問題。

(2)粒子近似

粒子近似是SPH方法近似過程的第二步，粒子近似的核心就是離散化，將積分運算離散為求和運算。式(3)通過粒子近似可推導為：

(7)

式中：j為粒子i支持域內的任一粒子;mj為粒子j的質量;ρj為粒子j的密度;j=1,2,…,N,N為粒子i支持域內粒子總數，如圖1(b)所示。因此，粒子i處的函數值可最終表達為：

(8)

式中：Wij=W(xi-xj,h)。式(6)也可通過粒子近似推導為：

W(x-xj,h)

(9)

粒子i處的函數空間導數值可最終表達為：

iWij

(10)

SPH方法實際是將對任意函數的運算轉化為對核函數的運算，由于核函數結構簡單,易于求導，從而達到了簡化運算的效果。

1.2 N-S方程求解

流體屬于連續介質，因此需遵守連續介質方程組，即Navier-Stokes方程組。N-S方程組包括：

1)連續性方程,

(11)

2)動量方程,

(12)

式中：t為時間；v為速度向量；ρ為密度；σ為應力張量；α和β均代表三個坐標軸方向。上述方程未體現外力及能量交換的情況，但并不影響后續SPH方法的應用。對N-S方程組分別進行SPH離散近似，將方程等式右側的偏微分項轉換為核函數的偏微分，從而得到SPH形式的離散化N-S方程組。

(1)連續性方程

應用式(10)對方程(11)等式右側的速度偏微分項進行SPH離散近似，可以得到：

(13)

由于單位空間導數存在如下表達形式：

W(x-x′,h)dx′=

(14)

經整理計算，可以得到SPH離散近似后的連續性方程：

(15)

(16)

方程(16)為最常用的SPH離散化連續性方程。從方程的結構可以看出，流體密度的時間變化率與粒子質量、粒子間的相對速度有關，同時也與核函數的偏導數有關。

(2)動量方程

應用式(10)對方程(12)等式右側的應力張量偏微分項進行SPH離散近似，可以得到：

(17)

由式(14)可知存在以下恒等式：

(18)

將式(18)與式(17)相加，則得到經SPH離散近似后的動量方程：

(19)

經整理后最終可以得到：

(20)

方程(19)與方程(20)均為常用的SPH離散化動量方程。從方程的結構可以看出，流體速度的時間變化率與粒子質量、粒子密度、粒子應力張量均有關，同時也與核函數的偏導數有關。

1.3 SPH方程相關計算

采用SPH方法通過離散近似得到了連續性方程、動量方程、能量方程的離散形式，為得到最終的計算結果，相關的計算問題需要得到進一步解決：核函數的構造與選擇、應力張量的分解與計算、邊界問題的處理等問題。

(1)核函數

核函數是SPH方法的核心概念，準確地構造滿足條件的核函數是SPH方法不斷發展的核心驅動力。本文采用B樣條核函數[17]。核函數的支持域由光滑長度h及光滑因子κ決定，支持域半徑為κh。B樣條核函數表達式：

W(R,h)=αd×

(21)

式中：αd=1/4πh3。

(2)應力張量

對于離散化連續性方程，核函數表達式確定后即可進行求解計算。但對于離散化動量方程，還需解決應力張量的計算問題。應力張量一般代表流體介質的本構方程，由壓強項及粘度項組成，目前主流的計算方案是通過狀態方程及人工粘度對其進行計算求解。以動量方程(19)為例，可將應力張量分解，得到如下形式的動量方程：

(22)

式中：pi,pj分別表示粒子i,j處的壓強；Πij代表粒子i,j間的人工粘度。

(3)邊界處理

為解決邊界問題，SPH方法引入虛粒子，這種類型的粒子一般位于問題域邊界之外，主要用于解決距問題域邊界一個支持域半徑長度內發生截斷的流體粒子的邊界問題。

圖2描述了SPH方法中的幾種粒子類型及其相對位置關系。其中，粒子3，8代表邊界粒子；粒子4，7代表發生截斷的流體粒子；粒子5，6代表未發生截斷的流體粒子；空心粒子(粒子1，2，9，10)代表虛粒子。粒子5，6距問題域邊界不小于kh，粒子4，7距問題域邊界小于kh，虛粒子位于問題域邊界以外。虛粒子一般不具有物理參數且不帶入SPH方程中進行計算，部分特殊虛粒子即使具有物理參數且代入SPH方程，下一仿真時刻虛粒子也不會按照計算結果移動，而是根據初始設置保持固定或按照一定規律運動。

不同邊界問題處理方案的主要差別在于設置具有不同作用的虛粒子。

圖2 支持域與問題域截斷導致邊界問題Fig.2 Boundary problem caused by the overlap of support domain and computainal domain

圖3所示為動態粒子邊界。在此邊界問題解決方案中，虛粒子被設置為多層動態粒子，虛粒子位于問題域邊界及邊界外的區域，虛粒子厚度與支持域半徑相關。動態粒子邊界的特殊性在于動態粒子與流體粒子共同代入SPH方程，用于場變量的離散近似計算，從而提高邊界處的SPH計算精度。但動態粒子在每時間步長下的計算結果不代入下一時間步長，而是按照仿真設置保持靜止或按照一定規律運動，這也使得動態粒子邊界適合用于含有復雜邊界或運動邊界的流體仿真問題。

圖3 動態粒子邊界Fig.3 Dynamic particle boundary

后續仿真任務涉及球形貯箱，并會在貯箱內設置不同類型的防晃擋板，因此本文選擇動態粒子邊界[18]作為邊界問題解決方案。

1.4 貯箱受力計算

液體晃動現象同樣對貯箱產生力的作用，進而影響攜帶貯箱的火箭、衛星等航天器的動力學行為。因此，需同樣監測貯箱的受力情況，設置一定的防晃措施，減小貯箱受力。建模仿真過程中，需要計算貯箱受力，貯箱結構由離散的邊界粒子等效，計算貯箱受力實際上是計算邊界粒子所受合力。邊界粒子作為虛粒子會被代入流體粒子的動量方程中進行計算，反過來同樣可以將流體粒子代入邊界粒子的動量方程進行計算，方程表達為：

(23)

式中：下標b和f分別代表邊界粒子量和流體粒子量。上式得到一個邊界粒子的數值加速度值，將所有邊界粒子的加速度與自身質量相乘再求和，即可得到貯箱的受力F：

(24)

2 仿真試驗結果與分析

2.1 模型建立

(1)液體體積

本文采用的球形貯箱半徑設置為r=0.125 m，貯箱設置為剛體。由于液體體積的變化會改變液體對貯箱的作用力，針對仿真需求設計如圖4所示的5種含有不同體積液體的球形貯箱模型。5種模型的貯箱運動模式設置為正弦運動，速度曲線方程為vy=0.12πcos(3πt)，且保持不變。定義液面高度h為貯箱內液面至貯箱底部的距離，通過設置不同的液面高度值實現不同的液體體積。液面高度分別設置為h=0.05 m,0.075 m,0.1 m,0.125 m,0.15 m，如圖4(a)～(e)所示，圖中數值單位為m。

(2)擋板設置

貯箱內的液體晃動現象需要受到抑制，是否設置擋板及擋板的安裝位置變化均會改變液體對貯箱的作用力，針對仿真需求設計如圖5所示的一種不含擋板及4種含有不同安裝位置擋板的球形貯箱模型。擋板同樣設置為剛體，且位于xoz平面。以上5種模型的貯箱運動模式設置為正弦運動，速度曲線方程為vy=0.12πcos(3πt)。液面高度設置為h=0.1 m，且均保持不變。定義擋板高度hb為擋板上邊至下邊的距離，擋板位置hp為擋板下邊至貯箱底部的距離。擋板高度設置為hb=0.05 m，且保持不變。擋板位置分別設置為hp=0 m,0.025 m,0.05 m,0.075 m，如圖5(a)～(e)所示，圖中數值單位為m。

圖4 不同液體體積貯箱模型Fig.4 Tank models of different liquid volume

圖5 不同擋板設置貯箱模型Fig.5 Tank models of different baffle positions

2.2 仿真結果分析

本文主要分析球形貯箱的受力情況及局部壓強情況。由于貯箱運動僅發生y軸方向，液體在yoz平面的晃動行為明顯，因此主要分析貯箱y,z兩軸方向的受力情況。貯箱內設置4個測量點，均位于yoz平面。測量點在高度范圍需分布均勻且涵蓋貯箱重點受壓區域，因此測量點按照圖4(a)～(e)及圖5(a)～(e)中所示位置進行設置，各仿真工況的粒子數及仿真時間如表1所示。

表1 粒子數及仿真時間

(1)液體體積

測量點p1處壓強情況如圖6(a)所示。液體體積越大(液面高度越高)，壓強越大，當h=0.05 m時，該測量點處存在周期性不受壓時段。對于5種不同液面高度，壓強曲線基本表現為“雙波峰”形態。測量點p2處壓強情況如圖6(b)所示。5種不同液面高度工況在該測量點處壓強值排序與測量點p1處壓強值排序相同，當h=0.05 m,0.075 m時，該測量點處存在周期性不受壓時段。隨著液面高度的降低，壓強曲線逐漸由“前波峰”加“后平臺”形態過度為“單波峰”形態。測量點p3處壓強情況如圖6(c)所示。5種不同液面高度工況在該測量點處壓強值排序與測量點p1,p2處壓強值排序相同，當h=0.05 m,0.075 m,0.1 m,0.125 m時，該測量點處存在周期性不受壓時段。對于5種不同液面高度，壓強曲線基本表現為“單波峰”形態。測量點p4處壓強情況如圖6(d)所示。5種不同液面高度工況在該測量點處壓強值排序與測量點p1,p2,p3處壓強值排序相同，且全部存在周期性不受壓時段。對于5種不同液面高度，壓強曲線同樣表現為“單波峰”形態。

貯箱y軸方向受力情況如圖7(a)所示。5種不同液面高度工況均存在受力方向周期性改變現象(力值為正表示力方向沿y軸正向，力值為負表示力方向沿y軸負向)。隨著液面高度的增加，y軸方向受力增大。5種不同液面高度工況y軸正負方向受力峰值均基本相同，受力曲線均表現為周期性均勻變化，無力震蕩情況。

圖6 不同液體體積測量點壓強Fig.6 Pressure of measuring point of different liquid volume

續圖6Fig.6 Continued

貯箱z軸方向受力情況如圖7(b)所示。不存在受力方向改變情況，受力方向保持為z軸負向。隨著液面高度的增加，z軸方向受力增大。5種不同液面高度工況受力曲線均表現為周期性均勻變化，無力震蕩情況。

通過對不同液面高度工況下貯箱受力及壓強情況仿真結果的描述可以作出如下分析：5種不同液面高度工況的根本區別在于液體的質量不同，液面高度越高，液體質量越大。而正是這個原因，導致了貯箱受力及壓強表現出仿真結果展示的規律和現象。不同測點位置的壓強區別與貯箱在y,z兩軸受力方向上的不同則體現出與正弦、階躍運動模式下的仿真結果相同的規律。不同液體體積晃動情況如圖8所示。

(2)擋板設置

測量點p1處壓強情況如圖9(a)所示。無擋板工況壓強最大，hp=0.05 m工況壓強最小，其他3種工況壓強值由大到小依次為：hp=0 m，hp=0.025 m，hp=0.075 m。當hp=0.075 m時，壓強曲線表現為“前波峰”加“后平臺”形態；當hp=0.05 m時，壓強曲線表現為“單波峰”形態；另3種工況下壓強曲線表現為“雙波峰”形態。測量點p2處壓強情況如圖9(b)所示。hp=0 m工況壓強最大，hp=0.05 m工況壓強最小，其他3種工況壓強值由大到小依次為：hp=0.025 m，無擋板，hp=0.075 m。當hp=0 m及hp=0.025 m時，壓強曲線表現為“雙波峰”形態；當hp=0.05 m時，壓強曲線表現為“單波峰”形態；另2種工況下壓強曲線表現為“前波峰”加“后平臺”形態。測量點p3處壓強情況如圖9(c)所示。5種不同擋板設置在該測量點處壓強值排序與測量點p1處壓強值排序相同，壓強曲線全部表現為“單波峰”形態，且全部存在周期性不受壓時段。測量點p4處壓強情況如圖9(d)所示。5種不同擋板設置在該測量點處壓強值排序與測量點p1,p3處壓強值排序相同。hp=0.05 m及hp=0.075 m兩工況在該測量點處仿真全程不受壓，另三種工況壓強曲線均表現為“單波峰”形態。

圖7 不同液體體積貯箱受力Fig.7 Tank force of different liquid volume

圖8 不同液體體積晃動情況Fig.8 Sloshing state of different liquid volume

圖9 不同擋板設置測量點壓強Fig.9 Pressure of measuring point of different baffle position

貯箱y軸方向受力情況如圖10(a)所示。5種不同擋板設置均存在受力方向周期性改變現象(力值為正表示力方向沿y軸正向，力值為負表示力方向沿y軸負向)。無擋板工況受力最大，hp=0.05 m工況受力最小，其他3種工況受力值由大到小依次為：hp=0 m，hp=0.025 m，hp=0.075 m。5種不同擋板設置y軸正負方向受力峰值均基本相同，受力曲線均表現為周期性均勻變化，無力震蕩情況。貯箱z軸方向受力情況如圖10(b)所示。不存在受力方向改變情況，受力方向保持為z軸負向。5種不同擋板設置在z軸方向受力值排序與在y軸方向受力值排序相同。5種不同擋板設置受力曲線均表現為周期性均勻變化，無力震蕩情況。

通過對不同擋板設置工況下貯箱受力及壓強情況仿真結果的描述可以作出如下分析：設置擋板的目的在于抑制晃動，4種含有擋板的工況基本起到了抑制晃動的效果。在測量點p2處，hp=0 m及hp=0.025 m工況的壓強都略大于無擋板工況，與期望效果正相反。這種增強晃動的現象僅發生于該測量點處的壓強情況，說明是一種特例，而原因在于擋板的存在會改變液體的流動方式，致使貯箱部分區域受到的壓強增大。抑制效果最好的工況為hp=0.05 m，這與參數設置存在一定關系，液面高度設置為h=0.1 m，擋板高度設置為hb=0.05 m，所以當擋板位置設置為hp=0.05 m時，貯箱內液體被分割為兩個相對獨立的空間，液體在其間不能自由流動，從而獲得了最好的抑制效果，導致了貯箱受力及壓強表現出仿真結果展示的規律和現象。不同測點位置的壓強區別與貯箱在y,z兩軸受力方向上的不同則體現出與正弦、階躍運動模式下的仿真結果相同的規律。不同擋板晃動情況如圖11所示。

圖10 不同擋板設置貯箱受力Fig.10 Tank force of different baffle position

圖11 不同擋板晃動情況Fig.11 Sloshing state of different baffle position

3 結束語

由于晃動問題發生于航天器球形貯箱，貯箱的特殊形狀使得液體極易發生強非線性晃動，傳統網格類計算流體動力學方法及等效模型方法均無法描述這一現象,因此本文采用無網格類的SPH方法進行建模分析。首先掌握SPH方法的核近似、粒子近似原理，然后將該原理應用于流體介質控制方程(N-S方程)，得到離散近似的連續性方程、動量方程。在此基礎上，通過選取核函數、拆分應力張量、處理邊界問題、確定鄰域粒子搜索方案等過程，最終建立基于SPH方法的流體動力學模型。通過對航天器球形貯箱不斷變化的燃料含量及防晃擋板的設置方式等方面的分析與討論，分別進行不同燃料液體體積及不同擋板安裝位置仿真分析。記錄仿真數據，分析貯箱受力及局部壓強，仿真結果表明：運動加速度越大、變化越突然，貯箱受力及壓強均越大、越易發生震蕩。液體體積越大，貯箱受力及壓強均越大。針對本文所建立的球形貯箱模型，設置擋板均會抑制液體的晃動行為，但在hp=0.05 m工況下，抑制效果最好。