唐浩等

摘 要：隨著航空航天技術的發(fā)展，研究液體晃動顯得越來越重要。液體晃動對飛行器燃油供給、飛行狀態(tài)以及機體結構變形產生較大影響。如何準確模擬液面變化、計算液體對壁面的沖擊力是目前存在的難點。針對這些問題，簡要介紹了SPH方法的基本原理，對導彈油箱內燃油的晃動進行了仿真計算，模擬出液面變化。計算了液體重心位置的變化曲線，同時得到油箱與燃油之間的相互作用力，并對結果做出了分析。對導彈油箱內燃油晃動研究提供了一定的理論和工程參考價值。

關鍵詞：光滑粒子動力學；燃油晃動；數值仿真；接觸力

貯箱內液體的晃動是指液體自由表面由于受到外加擾動或激勵而產生的運動。液體的晃動問題在航空航天、船舶、石油化工以及核動力等領域較為常見，從上世紀50年代起，各國工程師對這一問題已經給予了廣泛關注[1]。儲箱內帶自由液面的液體晃動不僅會對貯箱壁產生沖擊載荷，引起結構的疲勞破壞。對于載液比重大的載液系統(tǒng)，如航天器、大型民機以及油罐車等，在運動過程中，短時間內液體的大幅晃動引發(fā)的液體重心變化還會對這些載液系統(tǒng)的運動軌跡或飛行姿態(tài)產生不利的影響。依據相應的要求，油箱的幾何尺寸應避免在飛行中使燃油產生劇烈晃動，使燃油晃動造成的重心變化在可接受范圍內；由于飛機姿態(tài)變化所帶來的燃油重心移動也處于可以接受的范疇[2]。在航空領域，飛機在起飛、著陸及飛行過程中，因氣動紊流、發(fā)動機轉子的不平衡，武器的發(fā)射等原因，會引起油箱內燃油的晃動。研究晃動產生的沖擊力、晃動力矩所引發(fā)的油箱結構強度問題、液體晃動重心的變化對全機穩(wěn)定性的影響是十分必要的[3]。

液體晃動是一種特殊的流動現(xiàn)象，擁有自由表面流動的相關特點，同時還與約束液體運動的壁面有著較為復雜的相互作用，此類問題的難點就在于對液面大幅晃動的準確模擬以及對液體與壁面相互作用的處理如何準確模擬液面大幅晃動以及處理[4]。因此，對液體晃動問題的求解極具挑戰(zhàn)性。

光滑粒子力學（Smoothed Particle Hydrodynamic， SPH）方法廣泛用于求解流體動力學問題，可以很好的適應自由邊界的變化，并且對于自由表面大變形流動問題的處理效果也比較好。對于液體的邊界條件處理，本文考慮了自由液面以及液體與固體壁面兩種情況。自由邊界問題，文章采用Dirichlet方法處理，根據表面粒子密度來判斷特定粒子是否位于自由邊界上。而對于液體與固體壁面的相互作用的區(qū)域，則在固體壁面上布置一組虛粒子，用來對內部粒子產生排斥力，從而對粒子與固體壁面之間的相互作用進行模擬。

文章基于SPH方法數值模擬導彈油箱受到激勵時液體的大幅晃動，模擬了燃油液面變化，計算了燃油晃動過程中重心坐標的變化以及油箱和燃油相互作用力的變化，并對結果作了分析。

1 理論簡介

1.1 SPH方法基本原理

SPH的理論基礎是粒子方法[5-7]，使用一系列離散點或者粒子來描述連續(xù)的系統(tǒng)，系統(tǒng)中的離散的粒子將承載所有的物理信息。

對于粒子i的任意一個函數值都可以應用光滑函數對其支持域內所有粒子相應的函數值進行加權平均后獲得，支持域是一個半徑為kh的圓，k是確定粒子影響半徑的系數，如圖1所示，其中S是影響區(qū)域的邊界。

通過粒子近似式，可以將連續(xù)積分形式下的函數轉換成在支持域的粒子的離散化求和，通過這種方式實現(xiàn)無網格化。

1.2 核函數

本文研究選擇的核函數為Monghan和Lattanzio提出的三次樣條函數[8]，也稱B樣條函數。

為了滿足一致性要求，上式中的系數αd在一維、二維和三維空間中分別取為 和 。

1.3 控制方程

拉格朗日描述下的粘性流體控制方程就是Navier-Stokes方程。

用希臘字母上標α、β表示坐標系的方向，可以通過張量理論中的指標法用來表示方程的疊加，并且可以在變化的拉格朗日空間內引入對總時間的導數，將N-S方程作粒子近似，并進行空間離散，得到粒子近似的N-S方程。

連續(xù)方程：

動量守恒方程：

式中，σ為總應力張量，由各向同性壓力p和粘性應力 τ兩部分組成，v為速度，ρ為密度，x為距離，σ為總應力張量。

1.4 邊界條件

在液體和固體壁面的相互作用的區(qū)域內，液體不可以穿透貯箱壁，同時也不能脫離貯箱壁存在。因此在固壁邊界上布置一組虛粒子，對鄰近粒子施加排斥力，這樣可以防止內部粒子穿透固壁邊界，如圖2所示。排斥力的方向和兩個粒子的中心連線相同，排斥力的大小可以根據公式（5）計算得到[4]。

式中的n1和n2分別取值為12和4；r0是排斥力的作用范圍，其值為離散后粒子的初始間距；D為常數，大小與最大速度的平方相近。

處理自由表面邊界，則采用Dirichlet邊界條件[11]，根據粒子i的粒子密度來判斷是否位于自由邊界，若粒子 密度滿足式（6），那么認定其在自由表面上，并使該粒子密度等于流體表面的實際密度。

式中，ρ0是自由表面流體的實際密度，參數β的大小在0.85-0.98之間，本文將其取值為0.95。

2 模型的建立

2.1 模型說明

貯箱結構取自某型導彈的油箱[9-10]，幾何模型如圖3所示，左圖為貯箱結構圖，油箱的長約為830mm，寬約為650mm，高度約為560mm。油箱內部液體為煤油，液體體積約為0.148m3。

為方便將油箱內液體離散成均勻的SPH粒子，采用體積等效的方式，將油箱中的燃油用體積與其相等的形狀簡單的幾何體代替。計算得到油箱內燃油的體積約為0.148m3，為保證體積等效的幾何體完全包含于油箱內部，設計如圖4所示長方體形狀的水柱，燃油劃分為5460個SPH單元，油箱劃分為殼單元，不考慮油箱的變形，系統(tǒng)采用右手坐標系。

2.2 計算工況說明

油箱晃動的俯仰角度隨時間變化的曲線，如圖5所示?？紤]到燃油從圖4所示的液柱狀態(tài)運動到平穩(wěn)狀態(tài)需要一段時間才能穩(wěn)定，因此采用圖5的角度激勵曲線，初始的0-3.5s內俯仰角度為0，液柱在重力作用下運動。3.5s后燃油達到基本穩(wěn)定，俯仰角度從負方向的-45度線性變化到正方向的45度，變化時間為3s，6.5s后俯仰角度重新變?yōu)?。在計算過程中，豎直方向始終作用地面的重力加速度G，值取9.8m/s2。

2.3 仿真流程說明

文章的仿真計算基于PAM Crash軟件平臺。完成了模型相關參數的設置后，將其導入Crash Pam Solvers求解器中進行計算，時間步長取0.0001s，仿真計算10s后停止。

3 仿真計算結果

3.1 油箱內液面的變化

在不同時刻油箱內燃油晃動的液面變化情況如圖6所示，其中左側為文章采用的數值方法計算得到的液面位置。

前3.5s燃油在重力作用下趨于穩(wěn)定，液面平穩(wěn)。3.5s-7.5s時，箱體做俯仰運動，液面發(fā)生波動，SPH方法能夠比較真實地模擬出液面隨油箱運動的變化情況，晃動過程中產生的波浪也能比較準確的模擬出來。

3.2 燃油的重心位置變化

油箱中燃油在晃動過程中，粒子的運動是相當無序的。由于液體已經離散成一系列有質量的SPH單元，因此重心位置可以根據劃分得到的SPH單元加權求和得到。

燃油重心在X、Z三個方向的變化曲線，如圖7所示。由于模型是關于XZ平面對稱，同時貯箱系統(tǒng)也在XZ平面內運動，因此，重心在Y方向的變化不大，主要表現(xiàn)為X方向和Z方向的運動。

從以上重心位置的時間歷程曲線可知，0-3.5s內燃油從液柱狀態(tài)下落運動到穩(wěn)定狀態(tài)的過程中，X方向的坐標基本保持不變，Z方向坐標也基本保持不變。在3.5-6.5s內，油箱在外界激勵下運動，使得液體X方向的坐標有較大的變化，近似于一個與所給激勵類似的正弦變化，Z方向的坐標有一個比較小波動變化。6.5s之后，油箱停止運動，由于是繞Y軸的俯仰運動，箱體也對稱于XZ平面，因此停止運動后，重心在X和Z方向坐標基本無變化。

3.3 油箱與燃油相互作用力變化

油箱內燃油晃動過程中，燃油對油箱有不同程度的沖擊作用力，而且兩者之間的相互作用力是不斷變化的。由于對稱面XZ兩側燃油對油箱的沖擊作用力可以相互抵消，因此Y方向的作用力相對很小。X方向和Z方向相互作用力較大，變化曲線如圖8所示。

圖中可以看出，前3.5s內，燃油液柱運動到穩(wěn)定狀態(tài)的過程中，X方向的沖擊力除了一開始液體下落時一直在很小的量級波動，而Z方向的沖擊力逐漸增大?；蝿舆^程中，X方向作用力波動明顯，Y方向變化始終較小，Z方向也有較小波動。6.5s后X方向趨于0，Z方向逐漸減小到1130N左右，接近燃油重力。

4 結語

文章從貯箱內液體晃動問題研究的必要性出發(fā)，討論了研究的意義，簡要闡述了SPH方法的基本理論。隨后針對某型導彈油箱中的燃油晃動問題展開了仿真計算和分析。

通過文章的相關研究工作可以得出如下結論：

⑴油箱內燃油晃動的結果較好，并且晃動過程中的波浪也能比較準確地描述出來，對油箱內防晃板的布置優(yōu)化和燃油供給方式的選擇提供一定參考。

⑵計算得到了液體晃動過程中重心坐標的變化，有利于對飛行器飛行姿態(tài)進行準確的控制調整。

⑶由于油箱燃油體積為0.148m3，通過計算可知重力為1138N，與仿真得到的穩(wěn)定值1130N差別很小，驗證了仿真的可靠性。

文章的研究可以應用在航天器液體推進劑儲箱晃蕩或者大型LNG船儲箱液體晃蕩問題，對飛機燃油晃動等研究也可以提供一定的參考價值。

[參考文獻]

[1]劉富.貯箱內液體晃動動力學分析及結構防晃技術研究[D].南京：南京航空航天大學，2010.

[2]冷飛.基于SPH方法的飛機油箱燃油晃蕩研究[D].南京：南京航空航天大學，2009.

[3]夏恒新，王為，寶音賀西，等.帶球形底的圓柱容器中液體小幅晃動的實驗研究和有限元分析[J].動力學與控制學報，2009，7（1）：66-70.

[4]孫云舫.SPH方法的改進及其在水波數值模擬中的應用[D].天津：天津大學機械工程學院，2009.

[5]Jiannong Fang，Aurèle Parriaux，Martin Rentschler，Christophe Ancey.Improved SPH methods for simulating free surface flows[J].Applied Numerical Mathematics，2009（59）：251-271..

[6]Sudarshan Tiwari，J?rg Kuhnert.Modeling of two-phase flows with surface tension by finite pointset[J].Journal of Computational and Applied Mathematics，2007（203）：376–386.

[7]Youn Kyung Song，Kuang-An Chang，Yonguk Ryu，Sun Hong Kwon.Experimental study on flow kinematics and impact pressure in liquid sloshing[J].Exp Fluids，2013（54）：1592-1611.

[8]Monaghan J J，Lattanzio J C.A refined particle method for astrophysical problems[J].Astronomy and Astrophysics， 1985，149（1）：135-143.

[9]Takehiro Himeno，Toshinori Watanabe，Satoshi Nonaka，Yoshihiro Naruo，Yoshifumi Inatani，Hiroshi Aoki.Numerical and experimental investigation on sloshing in rocket tanks with damping devices[J].AIAA，2007-5557：1-17.

[10]Sudarshan Tiwari，Axel Klar，Steffen Hardt，Alexander Donkov.Coupled solution of the Boltzmann and Navier–Stokes equations in gas–liquid[J].Computers & Fluids， 2013，10（71）：283-296.

[11]Koshizuka S，Nobe A，Oka Y.Numerical analysis of breaking waves using the moving particle semi-implicit method[J]. International Journal for Numerical Methods in Fluid，1998，26（7）：751-769.