唐浩等

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

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

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

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

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

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

1 理論簡(jiǎn)介

1.1 SPH方法基本原理

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

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

通過粒子近似式，可以將連續(xù)積分形式下的函數(shù)轉(zhuǎn)換成在支持域的粒子的離散化求和，通過這種方式實(shí)現(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)引入對(duì)總時(shí)間的導(dǎo)數(shù)，將N-S方程作粒子近似，并進(jìn)行空間離散，得到粒子近似的N-S方程。

連續(xù)方程：

動(dòng)量守恒方程：

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

1.4 邊界條件

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

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

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

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

2 模型的建立

2.1 模型說明

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

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

2.2 計(jì)算工況說明

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

2.3 仿真流程說明

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

3 仿真計(jì)算結(jié)果

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

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

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

3.2 燃油的重心位置變化

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

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

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

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

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

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

4 結(jié)語

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

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

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

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

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

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

[參考文獻(xiàn)]

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

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

[3]夏恒新，王為，寶音賀西，等.帶球形底的圓柱容器中液體小幅晃動(dòng)的實(shí)驗(yàn)研究和有限元分析[J].動(dòng)力學(xué)與控制學(xué)報(bào)，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.