鄭亞雄，王明振，史圣哲

(中國(guó)特種飛行器研究所水動(dòng)力研究中心，湖北荊門(mén)448035)

基于FPM的液艙晃蕩研究

鄭亞雄，王明振，史圣哲

(中國(guó)特種飛行器研究所水動(dòng)力研究中心，湖北荊門(mén)448035)

摘要:隨著液艙晃蕩運(yùn)動(dòng)機(jī)理的研究深入，艙內(nèi)液體運(yùn)動(dòng)響應(yīng)和自由液面波形精細(xì)模擬成為研究的重點(diǎn)。本文應(yīng)用無(wú)網(wǎng)格方法中的FPM(Finite Point Method)來(lái)對(duì)液艙晃蕩中的動(dòng)邊界及自由液面大變形等問(wèn)題進(jìn)行數(shù)值求解和圖像模擬。以二維液艙晃蕩為例，校核運(yùn)動(dòng)邊界與內(nèi)部流體粒子點(diǎn)間的運(yùn)動(dòng)響應(yīng)，得到FPM對(duì)于邊界所受砰擊力和內(nèi)部興波波高的計(jì)算結(jié)果以及自由液面模擬圖像。將計(jì)算結(jié)果進(jìn)行無(wú)因次化處理，通過(guò)與試驗(yàn)結(jié)果進(jìn)行比對(duì)，驗(yàn)證了FPM在液艙晃蕩問(wèn)題中應(yīng)用的可行性和準(zhǔn)確性，也為FPM解決三維復(fù)雜液艙晃蕩問(wèn)題中提供鋪墊。

關(guān)鍵詞:FPM;液艙晃蕩;自由液面;動(dòng)邊界

Research on sloshing of liquid tank based on FPM

ZHENG Ya-xiong，WANG Ming-zhen，SHI Sheng-zhe
(Hydrodynamic Research Center，AVIC Special Vehicle Research Institute，Jingmen 448035，China)

Abstract:With the development of research on sloshing of liquid tank，the motion response and the free surface wave simulation have been focused．In this paper，F(xiàn)PM (Finite Point Method) has been applied into numerical calculation of sloshing in the moving boundary and simulation of free surface deformation．Taking the two－dimensional sloshing tank as an example，interactional response between moving boundary and internal fluid particle points has been checked．Slamming force on boundary，internal wave height and image of free surface have been attained．Taking the computational results into dimensionless，the comparison with the experimental results verify the feasibility and accuracy of application of FPM in the liquid sloshing，and pave the way for FPM to solve the complex 3D liquid sloshing．

Key words:FPM; sloshing of liquid tank; free surface;moving boundary

0　引言

隨著海上運(yùn)輸?shù)陌l(fā)展，液貨船得到廣泛的應(yīng)用，液艙內(nèi)液體晃蕩所引起的砰擊壓力對(duì)船舶性能有很大的影響?；问幰后w的運(yùn)動(dòng)具有很強(qiáng)的非線(xiàn)性特性，當(dāng)艙壁做規(guī)律性很強(qiáng)的簡(jiǎn)諧運(yùn)動(dòng)時(shí)，液艙內(nèi)液體受到簡(jiǎn)諧激勵(lì)，由此反饋激發(fā)的砰擊作用的持續(xù)時(shí)間和所產(chǎn)生的興波幅值跟著時(shí)間來(lái)變化。晃蕩所激發(fā)的砰擊力的幅值大小很難去做實(shí)時(shí)準(zhǔn)確的預(yù)報(bào)，只能從統(tǒng)計(jì)角度上進(jìn)行大小的估計(jì)。在容器晃蕩時(shí)，當(dāng)外界激勵(lì)頻率接近內(nèi)部帶自由液面液體的固有頻率時(shí)，液體的運(yùn)動(dòng)則會(huì)十分激烈，尤其是在接近最低固有頻率時(shí)，這會(huì)對(duì)容器壁面和容器結(jié)構(gòu)帶來(lái)嚴(yán)重的不利影響。對(duì)于液艙晃蕩的研究，主要集中在液艙所受的最大晃蕩沖擊載荷以及液艙晃蕩受力與時(shí)間之間的關(guān)系及其統(tǒng)計(jì)特征。

在早期受計(jì)算手段的限制，研究液艙晃蕩問(wèn)題都是用試驗(yàn)的方法來(lái)模擬和記錄的，但從70年代開(kāi)始，CFD技術(shù)的發(fā)展和興起，很多有效的數(shù)值計(jì)算方法得以實(shí)現(xiàn)，開(kāi)始用數(shù)值手段來(lái)進(jìn)行晃蕩問(wèn)題的模擬求解。有限元法、有限差分方法等被應(yīng)用于晃蕩問(wèn)題的求解中。通過(guò)很多研究者的計(jì)算模擬，發(fā)

現(xiàn)對(duì)晃蕩問(wèn)題模擬的最關(guān)鍵的問(wèn)題在對(duì)自由液面變化的精確捕捉。在整個(gè)模擬的時(shí)間過(guò)程中，自由液面的形狀和運(yùn)動(dòng)實(shí)時(shí)發(fā)生變化，同時(shí)由于壁面邊界也在做運(yùn)動(dòng)，相應(yīng)的液體在發(fā)生遷移，這造成了自由液面運(yùn)動(dòng)的強(qiáng)非線(xiàn)性特征，大大增加模擬仿真的難度。近年來(lái)，無(wú)網(wǎng)格方法逐漸興起，表現(xiàn)出有限元法和有限差分法等所不具備的優(yōu)勢(shì)［1］。FPM是無(wú)網(wǎng)格方法中很重要的一種，在計(jì)算流體中有著很好的適應(yīng)性，對(duì)于解決自由液面大變形問(wèn)題有其獨(dú)特之處。

在前人研究液艙晃蕩理論的成果中，給出不同浸深比下晃蕩的共振頻率與液艙形狀和大小之間的關(guān)系。對(duì)于二維簡(jiǎn)化的矩形液艙，其共振頻率公式可表述如下:

式中: fn為頻率; a為液艙內(nèi)的水深; b為液艙內(nèi)水的寬度; n為擬合的階數(shù)，其中一般情況都只按一階情況來(lái)處理即可; g為重力加速度。通過(guò)計(jì)算得到的二維液艙共振頻率，來(lái)制定進(jìn)行數(shù)值計(jì)算模擬時(shí)所采用的艙壁搖蕩頻率。

1　FPM計(jì)算概述

FPM (Finite Point Method )是由Onate和Idelsohn等［2－3］于1996年提出的。FPM采用了移動(dòng)最小二乘法(MLS)生成形函數(shù)，是不用預(yù)先生成網(wǎng)格作為背景的方法，利用配點(diǎn)法生成離散粒子點(diǎn)，將代數(shù)方程進(jìn)行離散求解，在計(jì)算流體力學(xué)領(lǐng)域中的應(yīng)用有獨(dú)特的優(yōu)勢(shì)。

FPM離散和近似的目的是為了求解流體力學(xué)控制方程［4］，控制方程包括質(zhì)量守恒方程、動(dòng)量守恒方程和能量守恒方程。其中離散和近似的方法簡(jiǎn)要介紹如下。

1.1配點(diǎn)法

FPM使用配點(diǎn)法來(lái)對(duì)空間場(chǎng)量進(jìn)行離散，根據(jù)幾何模型特征來(lái)確定所計(jì)算的流體域，在流體域內(nèi)產(chǎn)生計(jì)算節(jié)點(diǎn)進(jìn)行離散。由于計(jì)算流體域可能隨時(shí)間而變化，因此定義某時(shí)刻的計(jì)算流體域?yàn)棣? t)，邊界為Ω( t)。在計(jì)算域生成一系列的計(jì)算粒子點(diǎn)集，其中將計(jì)算粒子點(diǎn)進(jìn)一步劃分為計(jì)算域內(nèi)部的粒子點(diǎn)集，計(jì)算域邊界粒子點(diǎn)集。

離散得到的每個(gè)粒子點(diǎn)都攜帶相應(yīng)位置處的場(chǎng)量信息，粒子點(diǎn)的位置信息隨著時(shí)間實(shí)時(shí)發(fā)生變化。

1.2移動(dòng)最小二乘法(MLS)

分析流體域Ω中場(chǎng)變量u(x)，在x處u(x)的MLS近似表達(dá)式為

式中p(x)為對(duì)于2維空間坐標(biāo)xT=［x，y］的基函數(shù)。為保證基函數(shù)p(x)的最小完備性，采用Pascal三角形求解。式(2)中的a(x)為一系數(shù)向量，可表示為

注意，式(3)中的系數(shù)向量a(x)是x函數(shù)。系數(shù)a(x)可進(jìn)行加權(quán)離散取極小得到

式中: n為包含在權(quán)函數(shù)珚W(x－xi)≠0的x支持域中的計(jì)算節(jié)點(diǎn)數(shù); ui為u在x = xi處的計(jì)算點(diǎn)物理量。式(4)是一種泛函形式，主要是針對(duì)場(chǎng)函數(shù)的加權(quán)殘量，用粒子點(diǎn)坐標(biāo)參數(shù)和近似函數(shù)來(lái)得到場(chǎng)函數(shù)的近似形式。對(duì)泛函進(jìn)行求導(dǎo)如下:

導(dǎo)數(shù)為0的點(diǎn)即為駐點(diǎn)，可得到對(duì)空間場(chǎng)量分布的最佳近似。

1.3自由液面處理

在本文計(jì)算中，自由液面是FPM所需解決的主要問(wèn)題。在FPM中，自由液面通過(guò)自由邊界粒子來(lái)描述，在每個(gè)時(shí)間步長(zhǎng)開(kāi)始時(shí)，都需要進(jìn)行遍歷檢測(cè)自由面粒子點(diǎn)。在自由邊界上，表面張力占很大的比重，在壓力張量和表面張力張量之間存在著力的平衡。

τ·n= kκn。(6)

式中:τ為壓力張量; k為表面張力;κ為自由面曲率。壓力張量可以分解為壓力部分和衰減的壓力張量部分，后者可以進(jìn)一步表示為對(duì)稱(chēng)的速度變形張量，具體表示形式如下:

將式(7)乘以nT，得到壓力和速度的耦合方程為:

自由面粒子的速度－時(shí)間方程為

對(duì)式(9)需要補(bǔ)充不可壓縮條件·u = 0。通過(guò)上述自由液面粒子點(diǎn)運(yùn)動(dòng)方程即可對(duì)自由液面粒子點(diǎn)的運(yùn)動(dòng)進(jìn)行求解，使用粒子點(diǎn)對(duì)自由液面波形進(jìn)行精細(xì)化模擬。

2　二維液艙晃蕩計(jì)算模型

二維矩形艙模型如圖1所示。計(jì)算結(jié)果與Yang C［5］所得結(jié)果做比較。對(duì)液艙艙壁采用水平激勵(lì)，約束艙的運(yùn)動(dòng)，實(shí)現(xiàn)水平振蕩運(yùn)動(dòng)，其運(yùn)動(dòng)方程可表述為

x= Asin(ωt)，

式中: x為水平位移; A為振幅;ω為圓頻率，ω= 2π/T，T為周期［6］。

圖1　液艙晃蕩計(jì)算模型Fig.1　Computational model of sloshing tank

用式(10)來(lái)計(jì)算該液艙一階近似的共振周期

在計(jì)算中，取振幅A為0.05 m。在共振周期周?chē)?3個(gè)周期進(jìn)行計(jì)算，具體周期見(jiàn)表1，圖2中的計(jì)算圖像都是計(jì)算周期內(nèi)容。在觀測(cè)整個(gè)晃蕩過(guò)程中，壁面所受的橫向力Fx時(shí)間歷程，和觀測(cè)點(diǎn)處波高時(shí)間歷程，計(jì)算真實(shí)時(shí)間為20個(gè)周期，計(jì)算圖像如圖2所示。將12個(gè)工況下的結(jié)果進(jìn)行統(tǒng)計(jì)分析，得到每個(gè)工況下的最大橫向力和最大波高。此外，還取周期為1.265 s的工況，以驗(yàn)證共振周期(見(jiàn)圖3)。

3　計(jì)算圖像分析

從計(jì)算所得到的圖像可以看出，隨著周期的改變，其液面形狀也發(fā)生很大的改變，其中周期越大，所具有的完整波形就越少。由于圖片截取時(shí)刻的差異，在圖中波形存在不同。從總體上看，隨著周期的增長(zhǎng)，波長(zhǎng)幅值也在增大。

圖2　不同周期的計(jì)算圖像Fig.2　Image in different period

圖3　共振周期t = 1.265 sFig.3　Image in resonant period t = 1.265 s

此外，通過(guò)圖3的的測(cè)量證實(shí)，本二維矩形艙的共振周期為1.265 s，在此周期下，液艙內(nèi)共有4個(gè)完整波形。

4　計(jì)算結(jié)果分析

4.1激勵(lì)周期校對(duì)

由于在計(jì)算液艙晃蕩時(shí)，壁面運(yùn)動(dòng)與興波傳遞存在著一定的時(shí)間差。為了確定這個(gè)時(shí)間差，特選取如圖1中所示近壁面的A處作為測(cè)量點(diǎn)，測(cè)量不同時(shí)刻的波高和橫向力，通過(guò)波高的變化周期來(lái)得到內(nèi)部流場(chǎng)的周期。A是相對(duì)測(cè)量點(diǎn)，與兩側(cè)壁面保持固定位置，隨著液艙的晃動(dòng)而運(yùn)動(dòng)，不對(duì)內(nèi)部流場(chǎng)產(chǎn)生影響。

表1　給定周期與測(cè)量周期對(duì)比Tab.1 The comparison of given and measured period

從表1可看出，在計(jì)算過(guò)程中，預(yù)設(shè)周期與計(jì)算中實(shí)測(cè)周期相差很微小，保證了液面波形周期的可知性。從而也說(shuō)明，F(xiàn)PM計(jì)算方法對(duì)邊界的運(yùn)動(dòng)和內(nèi)流場(chǎng)激勵(lì)的實(shí)時(shí)性，數(shù)值延遲和失真現(xiàn)象幾乎可以忽略不計(jì)。

4.2最大橫向力與周期

為了更清晰地表征艙壁所受最大的橫向力與晃蕩激勵(lì)周期之間的變化關(guān)系，先將兩者分別進(jìn)行無(wú)因次化處理，如圖4所示。其中Tn為共振周期;ρ= 1 000 kg/m3; g= 9.8 m/s2; l= 1 m; b= 0.35 m。

圖4　最大橫向力與周期的曲線(xiàn)對(duì)比圖Fig.4 The comparison on the curves of maximum force with period

從圖4可看出，F(xiàn)PM方法計(jì)算結(jié)果與試驗(yàn)結(jié)果之間差異。大體上看，兩者走勢(shì)相同，其中當(dāng)最大橫向力超過(guò)峰值后，隨著周期的增大，F(xiàn)PM理論計(jì)算值與試驗(yàn)值之間的差異變得越來(lái)越小，誤差控制在10%以?xún)?nèi)。這表明FPM對(duì)于邊界所受力的計(jì)算合理可信，并未受自由液面變形的影響使計(jì)算誤差增大。

4.3最大波高與周期

同樣將最大波高與周期進(jìn)行無(wú)因次化處理，L為液艙的寬度，Tn為共振周期。

圖5　最大波高與周期的曲線(xiàn)對(duì)比圖Fig.5 The comparison on the curves of maximum wave height with period

從圖5可看出，對(duì)于最大波高的測(cè)量，F(xiàn)PM數(shù)值計(jì)算值與試驗(yàn)值的變化趨勢(shì)相同，在數(shù)值大小上存在一定差異，在周期較大的情況，誤差相對(duì)較小。兩者之間最大的不同在于峰值的不同，一方面是因?yàn)楸緢D采用折線(xiàn)圖，使得最大值的捕捉存在誤差，另一方面是在FPM計(jì)算中，是采用粒子點(diǎn)的位置來(lái)判斷波高位置，由于粒子點(diǎn)數(shù)目的限制，對(duì)于最高值的定位也就產(chǎn)生了不同。

5　結(jié)語(yǔ)

本文利用FPM在液艙晃蕩模擬過(guò)程中，其主要難點(diǎn)在于對(duì)劇烈變化的自由液面和動(dòng)邊界進(jìn)行模擬。通過(guò)上述計(jì)算分析，可得出如下結(jié)論:

1)通過(guò)監(jiān)視A處的變化周期，證實(shí)了在FPM中邊界運(yùn)動(dòng)與內(nèi)部流體粒子點(diǎn)運(yùn)動(dòng)相互作用的實(shí)時(shí)性;同時(shí)從計(jì)算圖像上可以看出，邊界例子點(diǎn)與內(nèi)部流體粒子點(diǎn)沒(méi)有發(fā)生干涉、分離等現(xiàn)象，自由液面形狀捕捉良好;

2)在液艙艙壁運(yùn)動(dòng)的激勵(lì)下，液艙內(nèi)的液體也隨之發(fā)生遷移運(yùn)動(dòng)，液面也在不斷的變化波形，并且波形的變化周期并不是與晃蕩周期相同，呈現(xiàn)處強(qiáng)烈的非線(xiàn)性特征;從計(jì)算數(shù)據(jù)上看，F(xiàn)PM對(duì)于動(dòng)邊界所受砰擊力及液體波高的計(jì)算誤差較小，適用于液艙晃蕩關(guān)鍵問(wèn)題的計(jì)算和模擬。

參考文獻(xiàn):

［1］張雄，劉巖．無(wú)網(wǎng)格法［M］．北京:清華大學(xué)出版社，2004: 137－210．

［2］ONATE E，IDELSOHN S，ZIENKIEWICZ O C，et al．A finite point method in computational mechanics．Applications to convective Transport and fluid flow［J］．Int．J．Numer．Methods Engrg，1996，39:3839－3866．

［3］ONATE E，IDELSOHN S，ZIENKIEWICZ O C，et al．A stbilized finite point method for analysis of fluid mechanics problems［J］．Comp．Meth．Appl．Mech．Eng．，1996，139: 315－346．

［4］王福軍．計(jì)算流體動(dòng)力學(xué)分析［M］．北京:清華大學(xué)出版社，2011:4－12．

［5］YANG C，LOHNER R．Computation of 3D flows with violent free surface motion［C］/ /Proceeding of the 15th ISOPE，Seoul，Korea．June，2005:19－24．

［6］祁江濤，顧民，吳乘勝．液艙晃蕩的數(shù)值模擬［J］．船舶力學(xué)，2008，15(8) :574－581．

［7］TIWARI S，KUHNERT J．Finite pointset method based on the projection method for simulations of the incompressible Navier-Stockes equations［J］．Springer LNCSE: Meshfree methods for Partial Differential Equations，Vol．26，Ed by M．Griebel，M．A．Schweitzer，2002．

［8］KUHNERT J．An upwind finite pointset method for compressible Euler and Navier-Stokes equations［J］．Springer Lecture Notes in Computational Science and Engineering: Meshfree Methods for Partial Differential Equations I，2002．

［9］LUO Yun-hua，ULRICH H C．A generalized finite difference method based on minmimizing global residual［J］．Comput．Methods Appl．Mech．Engrg，2002，191:1421－1438．

［10］LIU G R，LIU M B．Smoothed particle hydro-dynamics: a meshfree particle methods［M］．World Scientific，2003．

作者簡(jiǎn)介:鄭亞雄(1990－)，男，碩士，助理工程師，研究方向?yàn)榇八畡?dòng)力學(xué)。

基金項(xiàng)目:航空科學(xué)基金資助項(xiàng)目(20122305003)

收稿日期:2014－06－27;修回日期: 2014－09－09

文章編號(hào):1672－7649(2015) 07－0034－04doi:10.3404/j.issn.1672－7649.2015.07.008

中圖分類(lèi)號(hào):U661．32

文獻(xiàn)標(biāo)識(shí)碼:A