棱形液艙內(nèi)液體晃蕩問(wèn)題的SPH數(shù)值模擬

劉謀斌，周 冉，邵家儒

(中國(guó)科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100190)

在未充滿液體的容器內(nèi)，液體在外界激勵(lì)下的運(yùn)動(dòng)稱為液體晃蕩(或晃動(dòng)，sloshing)。實(shí)際工程中存在很多晃蕩現(xiàn)象，如載液貨船(航行中的大型原油貨輪(VLCC)、液化石油氣船(LPG)、液化天然氣船(LNG))、飛行中的火箭液體燃料艙、地震時(shí)核反應(yīng)爐和水庫(kù)等都涉及晃蕩問(wèn)題［1］?；问幜黧w由具有輕重差異、互不相溶的2種流體組成(一般是氣液兩相流)，因此液體晃蕩的一個(gè)顯著特點(diǎn)是存在自由液面或運(yùn)動(dòng)界面［2］。當(dāng)外界激振力振動(dòng)幅度較大，或者激振力頻率與晃蕩系統(tǒng)固有頻率接近甚至一致時(shí)，容器內(nèi)的液體可能產(chǎn)生劇烈的振蕩并對(duì)承載容器產(chǎn)生巨大的沖擊壓力。液體晃蕩作為一類相當(dāng)復(fù)雜的運(yùn)動(dòng)現(xiàn)象，已成為工程流體力學(xué)領(lǐng)域的一個(gè)基礎(chǔ)性與前沿性研究課題，具有極其重要的科學(xué)意義和工程應(yīng)用價(jià)值。

晃蕩研究最初是在航空航天與核工業(yè)領(lǐng)域里進(jìn)行，經(jīng)過(guò)幾十年的發(fā)展，已經(jīng)形成了三大類研究方法——理論研究、實(shí)驗(yàn)研究和數(shù)值計(jì)算。近年來(lái)隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，數(shù)值模擬方法在晃蕩問(wèn)題中得到了廣泛應(yīng)用，如有限元法、有限差分法、邊界元法及粒子法［3-4］等。Nakayama和Washizu［5］用有限元法計(jì)算了勢(shì)流的晃蕩問(wèn)題。朱仁慶等［6］采用VOF(volume of fluid)法模擬了盛液容器內(nèi)液體的二維晃蕩。Faltisen［7］用邊界元法計(jì)算了矩形容器的液體晃蕩，考慮了非線性自由表面條件并引入人工黏性克服瞬態(tài)項(xiàng)的影響。

光滑粒子動(dòng)力學(xué)SPH方法(以下簡(jiǎn)稱SPH方法)最早是由Gingold等［8］以及Lucy［9］各自獨(dú)立地提出，最初用于模擬分析天體物理問(wèn)題，是一種拉格朗日無(wú)網(wǎng)格粒子法［10］。它利用核函數(shù)對(duì)物理問(wèn)題進(jìn)行近似處理，把連續(xù)的物質(zhì)空間離散到一系列無(wú)序分布的可運(yùn)動(dòng)粒子上，用離散的粒子代替連續(xù)分布的流體，在模擬具有自由表面的液體晃蕩問(wèn)題中，粒子是運(yùn)動(dòng)的，可以有效捕捉液面的位置和運(yùn)動(dòng)，處理大變形問(wèn)題［11］。Iglesias等［12］用SPH方法對(duì)液艙晃蕩所產(chǎn)生的力矩幅值進(jìn)行了計(jì)算。崔巖等［13］運(yùn)用SPH方法對(duì)二維矩形水槽的縱蕩過(guò)程和縱搖過(guò)程進(jìn)行數(shù)值模擬，分析激勵(lì)頻率接近容器一階固有頻率時(shí)的晃蕩現(xiàn)象，討論縱蕩和縱搖綜合作用對(duì)流體運(yùn)動(dòng)的影響。Shao等［14］應(yīng)用改進(jìn)的SPH方法研究了二維矩形水箱內(nèi)的液體晃蕩，結(jié)果顯示改進(jìn)的SPH方法能夠獲取光滑、準(zhǔn)確的壓力場(chǎng)，能有效模擬大幅度液體晃蕩。

筆者應(yīng)用SPH方法對(duì)棱形液艙內(nèi)的液體晃蕩問(wèn)題進(jìn)行數(shù)值模擬，將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，并進(jìn)一步研究不同充液比、激勵(lì)頻率對(duì)晃蕩行為的影響。

1 SPH方法基本理論

在SPH方法中，計(jì)算域用一系列的粒子來(lái)描述，粒子是整個(gè)場(chǎng)變量近似的計(jì)算框架，每個(gè)粒子包含獨(dú)自的材料性質(zhì)如質(zhì)量、密度等，在內(nèi)外力的作用下按照守恒控制方程的規(guī)律運(yùn)動(dòng)。傳統(tǒng)SPH方法對(duì)偏微分方程的近似由核近似和粒子近似2部分組成。核近似是對(duì)函數(shù)及其導(dǎo)數(shù)函數(shù)進(jìn)行連續(xù)形式的近似，通過(guò)對(duì)函數(shù)本身與某個(gè)核函數(shù)之積的積分得到。粒子近似則是對(duì)粒子支持域內(nèi)的所有粒子加權(quán)求和。

SPH形式的控制方程為

式中:ρ——流體密度;t——時(shí)間;N——支持域內(nèi)的粒子數(shù);m——粒子質(zhì)量;v——粒子的速度;W——光滑函數(shù);p——壓強(qiáng);μ——?jiǎng)恿︷ば韵禂?shù);h——光滑長(zhǎng)度;x——粒子的位置;g——重力加速度;rij——粒子i和粒子j之間的距離。

傳統(tǒng)的SPH方法精度較低，尤其是在粒子分布不均勻的區(qū)域(如邊界處粒子支持域被截?cái)?，不能準(zhǔn)確重構(gòu)二次和線性函數(shù)，甚至無(wú)法精確重構(gòu)一個(gè)常數(shù)，導(dǎo)致傳統(tǒng)SPH方法求解精度和穩(wěn)定性降低。筆者采用密度修正和核梯度修正改進(jìn)傳統(tǒng)SPH方法［14］。通過(guò)恢復(fù)核函數(shù)在粒子分布不均勻區(qū)域的歸一化性質(zhì)對(duì)密度進(jìn)行修正，使密度趨近初始密度，滿足不可壓流體的特性。核梯度修正在二維情況下可以表示為

式中:Vj——粒子j的體積。

基于核梯度修正的SPH算法具有二階精度，并且可以在不改變傳統(tǒng)SPH程序整體框架的條件下方便實(shí)施。

在SPH數(shù)值模擬中，對(duì)于固壁邊界很難像網(wǎng)格方法那樣嚴(yán)格實(shí)施。邊界上或鄰近邊界處粒子的支持域往往被邊界截?cái)啵@些粒子通常只受到邊界內(nèi)粒子的影響，這種單邊影響作用會(huì)使誤差持續(xù)積累，影響最終的模擬結(jié)果。筆者應(yīng)用改進(jìn)的耦合動(dòng)力學(xué)邊界法［15］，用兩類粒子完成整個(gè)固壁邊界的構(gòu)建。第1層虛粒子布置在邊界處，對(duì)于靠近邊界的流體粒子將產(chǎn)生排斥力。另2層虛粒子布置在固壁邊界外部，初始時(shí)刻即生成，且不需與內(nèi)部流體粒子關(guān)于邊界鏡像對(duì)稱，位置不需要隨時(shí)間推進(jìn)而改變。

2 棱形液艙液體晃蕩

圖1 棱形液艙模型示意圖(單位:mm)Fig.1 Model of prismatic tank(units:mm)

隨著海洋運(yùn)輸?shù)陌l(fā)展和能源需求的提高，棱形液艙被廣泛應(yīng)用于大型載液貨船儲(chǔ)液系統(tǒng)設(shè)計(jì)中。Mikelis等［16］在實(shí)驗(yàn)中通過(guò)對(duì)棱形液艙施加橫搖及俯仰激勵(lì)，測(cè)量液體晃蕩對(duì)液艙側(cè)壁及頂部產(chǎn)生的沖擊壓力。液艙的晃蕩實(shí)驗(yàn)滿足模型相似律，即幾何相似、運(yùn)動(dòng)相似、動(dòng)力相似、邊界條件相似及初始條件相似。實(shí)驗(yàn)?zāi)Ｐ桶?∶40的比例制作，在液艙的側(cè)壁上布置壓力測(cè)試點(diǎn)。筆者采用二維模型，模擬和分析了有無(wú)內(nèi)部隔板時(shí)棱形液艙內(nèi)的液體晃蕩行為。圖1給出了液艙尺寸及文中計(jì)算工況所用到的壓力測(cè)試點(diǎn)布置。測(cè)點(diǎn)A、B、C分別距離液艙底部75 mm、113.5 mm、205 mm，壓力分別表示為P1、P2、P3。外加激勵(lì)為繞轉(zhuǎn)軸的橫搖運(yùn)動(dòng)，轉(zhuǎn)軸距離艙底194 mm，激勵(lì)運(yùn)動(dòng)方程為［17］

式中:θ0——最大轉(zhuǎn)角，取0.1 rad;ω——角速度。

2.1 不帶隔板棱形液艙液體晃蕩

筆者在2種工況下對(duì)不帶隔板的棱形液艙進(jìn)行數(shù)值模擬，追蹤壓力P1、P2、P3的時(shí)間歷程，并將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。

工況1中，外界激勵(lì)周期T=1.517s，充液比h/H=0.33(H為液艙的高度)，水體粒子數(shù)為9330，總粒子數(shù)為11439。圖2給出t=3.0～4.5s約一個(gè)周期的液面形狀變化。在不同的激勵(lì)周期和充液比等參數(shù)影響下，液體晃蕩可以表現(xiàn)為多種復(fù)雜的形式。常見(jiàn)的液體晃蕩現(xiàn)象有駐波、行進(jìn)波、水躍以及組合波，還有可能在沖擊作用下表現(xiàn)出液體的翻卷和飛濺等強(qiáng)非線性現(xiàn)象。圖2直觀地顯示了液艙在橫搖激勵(lì)作用下較明顯的行進(jìn)波及輕微的飛濺現(xiàn)象。圖3為壓力P1、P2、P3的歷程與實(shí)驗(yàn)結(jié)果對(duì)比，壓力的最大峰值及周期變化基本吻合，在時(shí)間上略有滯后可能是由于實(shí)驗(yàn)數(shù)據(jù)選取時(shí)刻上的誤差所致。由圖3可知，壓力呈現(xiàn)雙峰特性，這是晃蕩過(guò)程中液體遲滯特性的典型表現(xiàn)。前一個(gè)峰值是由液體對(duì)艙壁沖擊產(chǎn)生的高脈沖壓力，具有瞬時(shí)性，主要發(fā)生在水躍或行進(jìn)波中;后一個(gè)峰值也是由沖擊產(chǎn)生的脈沖壓力，但是時(shí)間相對(duì)較長(zhǎng)，峰值較小，主要是由于液體快速連續(xù)地作用在未被完全淹沒(méi)的艙壁，壓力具有連續(xù)性［17］。

圖2 工況1典型時(shí)刻的自由液面形狀變化Fig.2 Evolutions of free surface at typical time instants obtained from case 1

圖3 壓力P1、P2、P3的時(shí)間歷程Fig.3 Pressure time histories of three observation points，P1，P2，and P3

工況2中，外界激勵(lì)周期T=1.112s，充液比h/H=0.61，水體粒子數(shù)為19416，總粒子數(shù)為21525。圖4給出了t=3.0～4.1 s約一個(gè)周期的液面形狀變化。可以看到在大充液比下，由沖擊產(chǎn)生較劇烈的液體飛濺，并且模擬出液面的翻卷和碎波現(xiàn)象。圖5為P2的歷程與實(shí)驗(yàn)結(jié)果的對(duì)比，壓力的最大峰值及周期變化趨勢(shì)基本吻合，在時(shí)間上略有滯后可能是由于實(shí)驗(yàn)數(shù)據(jù)讀取時(shí)刻上的差異所致。通過(guò)與工況1中P2的對(duì)比，可以看到在大充液比情況下晃蕩劇烈，液體對(duì)艙壁產(chǎn)生較大的沖擊，壓力值明顯比工況1增大許多。

圖4 典型時(shí)刻工況2的自由液面形狀變化Fig.4 Evolutions of free surface at typical time instants obtained from case 2

2.2 帶隔板棱形液艙液體晃蕩

為進(jìn)一步研究液艙結(jié)構(gòu)對(duì)液體晃蕩特性的影響，筆者對(duì)帶隔板的棱形液艙液體晃蕩問(wèn)題進(jìn)行了數(shù)值模擬。隔板布置在液艙的中部，隔板高度為95 mm，見(jiàn)圖1。其他參數(shù)與不帶隔板的工況1情況相同。

圖6為t=3.8～4.5 s約一個(gè)周期的速度矢量圖。通過(guò)與不帶隔板的情況(圖2)進(jìn)行對(duì)比，可知加設(shè)隔板后液體晃蕩幅度大幅度減弱，晃蕩較平穩(wěn)，沒(méi)有出現(xiàn)強(qiáng)烈的沖擊壁面以及飛濺等現(xiàn)象。圖7為P2歷程與不加隔板情況的對(duì)比，明顯看到加設(shè)隔板后的壓力值大幅度減小，并且雙峰特性基本消失，只有一個(gè)比較光滑平穩(wěn)的壓力峰值。由此表明，隔板的阻撓減緩了液體向另一側(cè)的流動(dòng)，大幅度消除了液體對(duì)艙壁的沖擊作用，能夠有效地抑制液艙的晃蕩。在生產(chǎn)實(shí)際中，人們往往選用在艙內(nèi)加隔板，達(dá)到抑制晃蕩的目的。加設(shè)隔板不僅有防晃作用，而且不影響艙容，因此在生產(chǎn)實(shí)際中得到廣泛應(yīng)用。不過(guò)一些具體問(wèn)題還應(yīng)當(dāng)具體分析，不能一概而論［18］。

圖5 P2壓力歷程Fig.5 Pressure time history at point P2

圖7 P2壓力歷程對(duì)比Fig.7 Comparison of pressure time histories at point P2under different conditions

3 結(jié) 語(yǔ)

應(yīng)用改進(jìn)的SPH方法對(duì)棱形液艙中的液體晃蕩現(xiàn)象進(jìn)行了數(shù)值模擬與研究。首先對(duì)棱形液艙在外加正弦激勵(lì)作用下的自由液面晃蕩進(jìn)行了數(shù)值模擬，研究不同激勵(lì)周期及充液比的液體晃蕩問(wèn)題，模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，且能夠捕捉到壓力的雙峰特性及波浪翻卷、飛濺等復(fù)雜現(xiàn)象。另外，通過(guò)改變艙體結(jié)構(gòu)，對(duì)帶有隔板的棱形液艙進(jìn)行了研究，結(jié)果表明通過(guò)合理改變艙體結(jié)構(gòu)可以有效地抑制晃蕩現(xiàn)象。

［1］娜日薩.VLCC液艙晃蕩仿真及結(jié)構(gòu)強(qiáng)度評(píng)估方法研究［D］.哈爾濱:哈爾濱工程大學(xué)，2006.

［2］李曉明.船舶液艙晃蕩載荷及結(jié)構(gòu)強(qiáng)度的工程計(jì)算方法［D］.武漢:華中科技大學(xué)，2007.

［3］朱仁慶，吳有生，彭興寧，等.船舶液體晃蕩動(dòng)力學(xué)的研究方法及進(jìn)展［J］.華東船舶工業(yè)學(xué)院學(xué)報(bào)，1999，13(1):45-47. (ZHU Renqing，WU Yousheng PENG Xingning，et al.Study and advance of liquid sloshing dynamics of ship［J］.Journal of East China Shipbuilding Institute，1999，13(1):45-47.(in Chinese))

［4］陸志妹，范佘明.船舶液艙晃蕩研究進(jìn)展［J］.上海造船，2010(2):14-17.(LU Zhimei，F(xiàn)AN Sheming.Study advance of liquid sloshing of ship［J］.Shanghai Shipbuilding，2010(2):14-17.(in Chinese))

［5］NAKAYAMA T，WASHIZU K.Nonlinear analysis of liquid motion in a container subjected to forced pitching oscillation［J］. International Journal for Numerical Methods in Engineering，1980，15(8):1207-1220.

［6］朱仁慶，顏開(kāi)，吳有生.盛液容器內(nèi)液體二維晃蕩的數(shù)值模擬［J］.華東船舶工業(yè)學(xué)院學(xué)報(bào)，1998，12(2):14-21.(ZHU Renqing，YAN Kai，WU Yousheng.Two dimensions numerical simulation of liquid sloshing in tanks［J］.Journal of East China Shipbuilding Institute，1998，12(2):14-21.(in Chinese))

［7］FALTINSEN O M.A numerical nonlinear method of sloshing in tank with two dimensional flow［J］.Journal of Ship Research，1978，22(3):193-202.

［8］GINGOLD R A，MONAGHAN J J.Smoothed particle hydrodynamics-theory and application to non-spherical stars［J］.Monthly Notices Royal Astronomical Society，1977，181:375-389.

［9］LUCY L B.A numerical approach to the testing of the fission hypothesis［J］.Astronomy Journal，1977，82:1013-1024.

［10］LIU Guirong，LIU Moubin.Smoothed particle hydrodynamics:a meshfree particle method［M］.Singapore:World Scientific，2003.

［11］謬吉倫，陳景秋，張永祥.SPH方法在自由表面流體研究中的應(yīng)用［J］.水利水電科技進(jìn)展，2011，31(3):20-23.(MIU Jilun，CHEN Jingqiu，ZHANG Yongxiang.Application of SPH method to studies on free surface flows［J］.Advances in Science and Technology of Water Resources，2011，31(3):20-23.(in Chinese))

［12］IGLESIAS A S，DELORME L，ROJAS L P，et al.Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics［J］.Ocean Engineering，2006，33(11/12):1462-1484.

［13］崔巖，吳為，龔凱，等.二維矩形水槽晃蕩過(guò)程的SPH方法模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展，2008，23(6):618-624.(CUI Yan，WU Wei，GONG Kai，et al.Numerical simulation of sloshing in two dimensional rectangular tanks with SPH［J］.Chinese Journal of Hydrodynamics，2008，23(6):618-624.(in Chinese))

［14］SHAO Jiaru，LI Huiqi，LIU Guirong，et al.An improved SPH method for modeling liquid sloshing dynamics［J］.Computer and Structures，2012，100-101:18-26.

［15］LIU Moubin，SHAO Jiaru，CHANG Jianzhong.On the treatment of solid boundary in smoothed particle hydrodynamics［J］. Science China Technological Sciences，2012，55(1):244-254.

［16］MIKELIS N E，MILLER J K，TAYLOR K V.Sloshing in partially filled liquid tanks and its effect on ship motions:numerical simulations and experimental verification［J］.The Royal Institution of Naval Architects，1984，126:267-277.

［17］劉富，童明波，陳建平.基于SPH方法的三維液體晃動(dòng)數(shù)值模擬［J］.南京航空航天大學(xué)學(xué)報(bào)，2010，42(1):122-126. (LIU Fu，TONG Mingbo，CHEN Jianping.Numerical simulation of three-dimensional liquid sloshing based on SPH method［J］. Journal of Nanjing University of Aeronautics＆Astronautics，2010，42(1):122-126.(in Chinese)

［18］端木玉，朱仁慶，陳正云，等.不同液艙結(jié)構(gòu)形式對(duì)晃蕩的影響分析［J］.水動(dòng)力學(xué)研究與進(jìn)展，2006，21(6):760-769. (DUAN Muyu，ZHU Renqing，CHEN Zhengyun，et al.Analysis on sloshing of different ship tanks［J］.Chinese Journal of Hydrodynamics，2006，21(6):760-769.(in Chinese))