導(dǎo)彈油箱燃油晃動仿真分析

唐浩等

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

關(guān)鍵詞：光滑粒子動力學(xué)；燃油晃動；數(shù)值仿真；接觸力

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

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

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

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

1 理論簡介

1.1 SPH方法基本原理

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

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

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

1.2 核函數(shù)

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

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

1.3 控制方程

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

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

連續(xù)方程：

動量守恒方程：

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

1.4 邊界條件

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

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

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

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

2 模型的建立

2.1 模型說明

貯箱結(jié)構(gòu)取自某型導(dǎo)彈的油箱[9-10]，幾何模型如圖3所示，左圖為貯箱結(jié)構(gòu)圖，油箱的長約為830mm，寬約為650mm，高度約為560mm。油箱內(nèi)部液體為煤油，液體體積約為0.148m3。

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

2.2 計算工況說明

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

2.3 仿真流程說明

文章的仿真計算基于PAM Crash軟件平臺。完成了模型相關(guān)參數(shù)的設(shè)置后，將其導(dǎo)入Crash Pam Solvers求解器中進(jìn)行計算，時間步長取0.0001s，仿真計算10s后停止。

3 仿真計算結(jié)果

3.1 油箱內(nèi)液面的變化

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

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

3.2 燃油的重心位置變化

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

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

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

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

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

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

4 結(jié)語

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

通過文章的相關(guān)研究工作可以得出如下結(jié)論：

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

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

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

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

[參考文獻(xiàn)]

[1]劉富.貯箱內(nèi)液體晃動動力學(xué)分析及結(jié)構(gòu)防晃技術(shù)研究[D].南京：南京航空航天大學(xué)，2010.

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

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

[4]孫云舫.SPH方法的改進(jìn)及其在水波數(shù)值模擬中的應(yīng)用[D].天津：天津大學(xué)機(jī)械工程學(xué)院，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.