含隔板球形貯箱液體晃動激勵頻率影響分析

馬 亮， 劉 昊， 魏 承， 湯 亮， 趙 陽

(1.哈爾濱工業大學 航天工程系，哈爾濱 150001； 2.北京控制工程研究所，北京 100190)

光滑粒子流體動力學方法(Smoothed Particle Hydrodynamics,SPH)自1977年由Lucy等[1-2]首次提出，在經歷了40年的發展和改進過程后，現已在航空航天[3-5]、船舶海洋[6-8]、公路運輸[9-11]等多個領域取得了廣泛應用。

在航天領域，隨著現代航天器對高精度載荷及長壽命飛行的深入需求，航天器貯箱內液體晃動所產生的晃動力對系統動力學行為的影響已經成為總體設計重點考慮的問題。

基于SPH方法的流體動力學問題研究經歷了逐步完善的過程。Monaghan[12]于1994年首次將SPH方法引入流體動力學問題的求解，隨后多位學者基于此方法對鈍體繞流[13]、剪切流動[14]、潰壩[15]等問題進行了仿真分析。同時，SPH方法在計算精度[16]、計算穩定性[17]及邊界處理[18]等方面經過不斷地改進和修正，使其成為一種解決流體動力學問題的成熟建模方法。

現代航天器貯箱多為球形結構，而針對球形貯箱的液體晃動問題，現有研究多采用等效模型法或有限元法：黃華等[19]通過建立三維質心面等效模型，利用質心動力學方程及接觸碰撞模型對橢球形貯箱的晃動力進行求解；趙陽等[20]采用任意拉格朗日-歐拉法(Arbitrary Lagrange-Euler,ALE)，對常重力及微重力環境下球形貯箱內的液體晃動現象進行了仿真分析；王天舒等[21]提出適用于微重力環境的球形貯箱液體晃動彈簧-質量等效模型，同時利用虛功率原理對晃動力進行了求解；苗楠等[22]建立了靜平衡面垂直于等效重力的等效模型，從而將液體大幅晃動分解為等效重力的整體運動及在此基礎上的小幅晃動；岳寶增等[23]提出網格移動算法，在ALE方法基礎上對球形貯箱大幅晃動進行了數值模擬。

由于航天器在執行任務過程中，伴隨著頻繁的軌道及姿態機動，使得航天器貯箱存在大范圍運動或大載荷受力，從而導致貯箱內液體出現大幅度晃動。現有研究所采用的等效模型方法[24]或有限元方法均無法適用于大幅晃動下液體出現破碎的現象，而針對航天器球形貯箱的液體大幅晃動問題又亟待解決。

通過以上分析，本文開展了基于SPH方法的球形貯箱液體晃動問題研究：首先根據SPH方法的基本方程，給出N-S方程的求解形式，同時確定人工黏度項及壓力項方程，并采用動態邊界以提高邊界計算精度；然后基于SPH方法建立球形貯箱液體晃動模型，采用兩種激勵模式及三種激勵頻率進行仿真分析，分別記錄艙壁受力和觀測點壓強情況；最終設計一種帶隔板的球形貯箱結構，并分析隔板對液體晃動的抑制作用。

1 基于SPH方法的N-S方程求解

光滑粒子流體動力學方法的基本思想是通過兩步近似，將空間函數及空間函數導數表示為核函數支持域內各離散點處函數值與核函數及核函數導數乘積求和的形式。有利于簡化空間函數及空間函數導數求解的難度，使得該方法適用于求解N-S方程。

由Gingold等的研究可知，在經歷核近似、粒子近似兩步近似過程后，空間函數f(x)在粒子i處的函數值可以表示為

(1)

式中:W為核函數;h為光滑長度;N為核函數支持域內粒子總數。

空間函數f(x)的導數在粒子i處的函數值可以表示為

(2)

核函數的支持域由光滑長度h及光滑因子κ決定，支持域半徑為κh，如圖1所示。本文采用三次樣條核函數[25]，其具體表達式為

式中:R=|x-x′|/h;光滑因子取為κ=2。

圖1 核函數支持域及粒子近似Fig.1 Support domain of kernel function and the particle approximation

采用SPH方法的基本方程，可將N-S方程組中的連續性方程、動量方程表示為

(4)

(5)

式中:vij=vi-vj;Wij=W(xi-xj,h);fe為單位質量外力;α和β為坐標軸方向。

(6)

在弱可壓縮模型下[27]，壓強項可表示為

(7)

式中:ρ0為流體參考密度；壓強-密度剛性系數γ=7;聲速cs取決于流體的最大速度，取為cs=10Vmax。

SPH方法存在邊界問題：在距邊界一個κh長度的范圍內(包括邊界處)會出現核函數支持域與問題域交錯的現象。如圖2所示，方形粒子代表邊界粒子，三角形粒子代表發生截斷的流體粒子，圓形粒子代表未發生截斷的流體粒子。為使這一范圍內的粒子計算保持精度，需對邊界問題進行特殊處理。

為解決SPH方法的邊界問題，本文采用動態邊界法[28]，該方法的思想是直接構造與流體實粒子性質相同的邊界虛粒子，在計算過程中將邊界虛粒子代入計算，解決邊界附近(一個κh長度范圍內)流體粒子核函數支持域截斷問題。但邊界虛粒子自身參數不隨計算結果發生變化，而是保持固定或隨規劃運動軌跡進行運動，該方法下虛粒子厚度與核函數支持域半徑相關。

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

2 仿真工況設計

采用SPH方法對N-S方程進行離散化，并給出人工黏度項、壓力項的經驗方程，針對邊界問題提出動態邊界法，從而實現對N-S方程的完整求解。在此基礎上完成球形貯箱液體晃動仿真模型的建立，主要包括貯箱及液體建模、激勵形式及工況設計。

2.1 貯箱及液體

貯箱設計分為無隔板模型及有隔板模型，兩種結構的唯一區別在于隔板結構，其他參數均相同。球形貯箱半徑r=0.125 m，液體深度h=0.1 m，隔板高度hb=0.05 m，隔板位于xoz平面內。為測量液體晃動對艙壁產生的壓力，設置了4個測量點，測量點均位于yoz平面的艙壁面上，不同測量點的高度不同。其中，p2點高度與隔板高度相同，p3點高度與液面高度相同，具體位置如圖3所示，貯箱與隔板均設計為剛體。

圖3 貯箱及液體建模(m)Fig.3 Model of spherical tank and fluid(m)

2.2 激 勵

由于航天器貯箱在實際飛行過程中存在大范圍運動或大載荷受力，仿真設計兩種類型激勵模式：力模式、運動模式。力模式的加載對象為液體，控制量為液體所受加速度；運動模式的加載對象為貯箱，控制量為貯箱位移。兩種模式下均采用簡諧激勵形式：x=Asin(2πf·t+φ)，力模式下x表示液體y方向所受加速度，運動模式下x表示貯箱y方向位移，仿真的重力環境選為常重力。設計6組仿真工況，激勵模式及簡諧運動參數如表1所示。

表1 仿真工況設計Tab.1 Design of simulation and parameters

3 激勵頻率影響分析

圖4所示為兩種貯箱模型下，4個時刻點處晃動液體的側視圖。每個時刻下有隔板貯箱的液體晃動幅度均小于無隔板情況，為了定量分析隔板影響，分別從艙壁受力、測點壓強兩方面進行對比分析。

3.1 艙壁受力

力模式激勵下液體的加載方向為y方向，液體在yoz平面的晃動行為明顯，因此主要分析艙壁y,z兩方向的受力情況。

艙壁y方向初始受力為零，隨著液體受到簡諧激勵，開始出現周期性晃動力。不同的激勵頻率造成艙壁受力不同的曲線形式，當激勵頻率f=1.5 Hz時晃動力幅度達到最大，而三種激勵頻率下隔板的存在均使得晃動幅度減小，如圖5(a)和圖5(b)所示。

艙壁z方向初始受力為負值，這是由液體重力造成的，隨著液體受到簡諧激勵，同樣出現周期性晃動力。不同的激勵頻率對艙壁z方向晃動力影響巨大，且同樣在激勵頻率f=1.5 Hz時晃動力幅度達到最大，隔板的存在同樣使得晃動幅度減小，如圖5(c)和圖5(d)所示。可以看出，由于存在液體晃動現象，艙壁在y,z兩方向的受力均增大，且為周期性受力。

表2給出了各工況下，艙壁受力穩定后的平均晃動幅值，體現出的規律與曲線圖一致，同時可以看出隔板對激勵頻率大的工況晃動抑制效果更明顯。

圖4 晃動過程中各時刻下液面形態Fig.4 Shape of the fluid at different moment in sloshing

圖5 艙壁受力隨振動頻率變化情況Fig.5 Force applied on bulkhead varies with vibrational frequency

表2 艙壁受力平均晃動幅值Tab.2 Average force amplitude applied on bulkhead

3.2 測點壓強

運動模式激勵下貯箱的加載方向同樣為y方向，液體在yoz平面的晃動行為明顯，因此在yoz平面設置了4個壓力測量點：p1點高度低于隔板高度，p2點高度與隔板高度相同，p3點高度與液面高度相同，p4點高度高于液面。

對于p1點，在無隔板工況下，不同的激勵頻率使得壓力曲線呈現出截然不同的形式：當f=1 Hz時，晃動行為與激勵形式一致，表現為簡諧波動；當f=1.5 Hz時，出現“雙波峰”現象，前壓力峰值略高于后壓力峰值，這是由液體撞擊艙壁后沿球形內輪廓上升后回落所形成；當f=2 Hz時，也會出現“雙波峰”現象，但兩次壓力峰值相差較大，這是由于液體撞擊艙壁后的上升幅度較小，沒有形成翻卷、破碎，如圖6(a)所示。在有隔板工況下，三種激勵頻率下的測點壓力值均略有下降。當f=1 Hz時，壓力曲線形式同樣表現為簡諧波動；當f=1.5 Hz時，“雙波峰”的前后壓力峰值更加接近；當f=2 Hz時，后壓力峰值反高于前壓力峰值，如圖6(b)所示。這是由于p1點位于隔板下方，隔板的存在會對液體首次沖擊艙壁存在阻礙作用，而當液體回落時隔板同樣阻礙了液體的回流，正是這種作用使得前壓力峰值有所減小而后壓力峰值明顯增大。但隔板的存在使得壓力峰值達到“平均化”效果，這樣不會使艙壁局部在某時刻的壓力過大，而導致結構變形或損壞。

對于p2點，其高度高于p1點，但仍位于液面之下，因此在有無隔板及三種激勵頻率下的壓力曲線形式與p1點都相似，但局部現象存在不同之處。在無隔板及激勵頻率為1.5 Hz,2 Hz的工況下，沒有出現明顯的“雙波峰”，而變為前波峰加“后壓力平面”的組合,如圖6(c)所示。這是由于p2點高度增加，液體不存在回落或回落幅度極小；在有隔板及激勵頻率為1.5 Hz,2 Hz的工況下，前壓力峰值高于后壓力峰值，如圖6(d)所示。這是由于p2點高度與隔板高度相同，隔板對液體首次沖擊艙壁及液體回落回流的阻礙作用都減小，但隔板的存在同樣使得測點壓力值有所下降。

對于p3點，其高度位于液面處，壓力曲線形式與前兩個位于液面下的測點存在較大差異。在三種激勵頻率下壓力曲線均呈現“單波峰”加零壓力組合。當不存在隔板時，激勵頻率越大壓力峰值越大,如圖6(e)所示；當存在隔板時，隨著激勵頻率的增加隔板對測點處壓力的抑制效果越明顯，如圖6(f)所示。

對于p4點，其高度位于液面之上，壓力曲線形式與p3點相似。在激勵頻率為1 Hz工況下已檢測不到壓力值，而在激勵頻率為1.5 Hz,2 Hz的工況下，壓力曲線也變為“脈沖波”加零壓力組合。隔板的存在同樣對激勵頻率越大的工況測點處壓力抑制越明顯,如圖6(g)和圖6(h)所示。

圖6 測點壓強隨振動頻率變化情況Fig.6 Pressure of measuring point varies with vibrational frequency

4 結 論

通過建立有無隔板兩種球形貯箱模型，針對兩種激勵模式及三種激勵頻率進行仿真分析，最終得到以下結論：

(1) 力模式激勵下，無論貯箱有無隔板均在激勵頻率f=1.5 Hz時，得到艙壁受力的最大峰值，而隔板的存在會減弱艙壁所受晃動力。

(2) 力模式激勵下，由于液體晃動現象的產生，艙壁在y,z兩方向的受力均增大，且為周期性受力。

(3) 運動模式激勵下，位于液面下的艙壁測點p1,p2表現出一致性，液體沖擊艙壁及回落作用導致“雙波峰”現象的出現，測點高度的不同又會使得該現象在兩測點處存在差異。而隔板的存在對艙壁所受壓力具有“平均化”效果。

(4) 運動模式激勵下，位于液面處及液面以上的艙壁測點p3,p4表現出一致性，出現“單波峰”或“脈沖波”加零壓力組合形式，且隔板的存在對激勵頻率越大的工況抑制效果越明顯。