去奇異邊界元方法在液艙晃蕩模擬中的應(yīng)用

王慶豐，徐 剛，王樹(shù)齊，朱仁慶

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

在船舶領(lǐng)域的應(yīng)用中，頻域理論已經(jīng)被運(yùn)用的相當(dāng)成熟，但是頻域法通常只適用于周期穩(wěn)態(tài)問(wèn)題，即它存在一個(gè)本質(zhì)上的限制就是對(duì)于瞬變或者強(qiáng)非線性問(wèn)題往往無(wú)法處理。而頻域理論的超集實(shí)際上是時(shí)域理論，它是一種不可分離時(shí)間項(xiàng)的方法，理論上可以解決物體任意運(yùn)動(dòng)的問(wèn)題，還能解決完全非線性問(wèn)題，在自由度方面存在明顯的優(yōu)勢(shì)。本文采用時(shí)域邊界元方法進(jìn)行數(shù)值模擬，但是對(duì)于邊界元方法，奇點(diǎn)的處理一直是關(guān)注的重點(diǎn)，傳統(tǒng)方法在交界面處的速度精度很難保障，為了消除積分的奇異性，學(xué)者們已經(jīng)研究出諸多成熟可用的方法，其一徑向積分邊界元法[1-2]，它基于純數(shù)學(xué)處理技術(shù)，不借助任何特解和微分算子可將任何問(wèn)題的域積分轉(zhuǎn)化為邊界積分，是一種無(wú)內(nèi)部網(wǎng)格的純邊界元算法，目前已應(yīng)用于斷裂力學(xué)、熱輻射等方面；其二就是無(wú)奇異邊界元方法，即將傳統(tǒng)方法中原本直接布置于流體計(jì)算域表面網(wǎng)格中心點(diǎn)上的所有奇點(diǎn)，移至計(jì)算域外部，可以減少存儲(chǔ)需求并提高計(jì)算效率，目前該方法主要應(yīng)用于船舶與海洋工程水動(dòng)力學(xué)分析方面。為了能夠更好地模擬液艙內(nèi)液體的運(yùn)動(dòng)，本文采用基于去奇異邊界元方法(Desingularized Boundary Integral Equation Method,DBIEM)研究了液艙晃蕩的基本數(shù)值解法。

傳統(tǒng)的邊界元積分方法，其奇點(diǎn)是分布在計(jì)算域的表面，在進(jìn)行相關(guān)數(shù)值計(jì)算的過(guò)程中需要對(duì)奇點(diǎn)積分進(jìn)行特殊的處理，表現(xiàn)在時(shí)域分析非線性問(wèn)題的時(shí)候這過(guò)程中會(huì)出現(xiàn)大量繁瑣的計(jì)算分析過(guò)程。當(dāng)我們采用無(wú)奇異邊界元方法的時(shí)候會(huì)得到兩個(gè)優(yōu)點(diǎn)[3]：① 可以得到控制面上高精度的計(jì)算結(jié)果；② 獲得離散邊界積分問(wèn)題代數(shù)系統(tǒng)的計(jì)算時(shí)間變少了。同時(shí)，該方法簡(jiǎn)單易實(shí)施。Jensen 等利用去奇異方法來(lái)分析在穩(wěn)定入射波情況下求解阻力的問(wèn)題[4]。Lee[5]將去奇點(diǎn)邊界積分方程和時(shí)間步進(jìn)技術(shù)進(jìn)行結(jié)合，模擬了簡(jiǎn)單潛體及浮體運(yùn)動(dòng)的完全非線性效應(yīng)問(wèn)題；李宗翰等[6]將去奇異積分方程方法運(yùn)用到任意三維浮體的完全非線性流場(chǎng)的模擬和仿真中；Scorpio等[7-8]對(duì)完全非線性數(shù)值水池進(jìn)行了探索；Kara等[9]利用去奇異邊界元方法在三維時(shí)域內(nèi)研究了完全非線性波物相互作用的問(wèn)題；Zhang等[10]基于去奇異方法研究了雙體船的耐波性；Xu等[11]進(jìn)行了隨機(jī)波浪下三維液艙的全非線性數(shù)值分析。

1 去奇異邊界元理論

1.1 控制方程及初始條件

根據(jù)勢(shì)流基本理論和圖1液艙晃蕩系統(tǒng)，速度勢(shì)φ在流體域D中滿足Laplace方程，定解條件如下

(1)

(2)

(3)

(4)

式中:η代表波高;x，y，z代表坐標(biāo)軸三個(gè)方向上的物理量。初始條件條件如下(φ0是初始速度勢(shì)、ξ是初始波高)

φ(x,y,z=ξ,t=0)=φ0(x,y,z)

(5)

η(x,y,t=0)=ξ(x,y)

(6)

圖1 液艙晃蕩坐標(biāo)系統(tǒng)Fig.1 The coordinate system of sloshing tank

1.2 自由面條件積分格式

有限差分法一般被運(yùn)用在自由表面上的速度勢(shì)φ的求解，有限差分法的不足是求解多維問(wèn)題時(shí)不易適應(yīng)不規(guī)則幾何邊界，并且在求解過(guò)程中會(huì)常常出現(xiàn)數(shù)值誤差的累積，從而導(dǎo)致計(jì)算結(jié)果發(fā)散，影響數(shù)值的穩(wěn)定性。本文采用一種新的積分格式的自由面條件(Integral Form of Free Surface Boundary Condition[12]，IFBC )，運(yùn)用此積分格式的自由面條件，能夠在線性自由面條件的基礎(chǔ)上得到自由面上一階速度勢(shì)的表達(dá)式，而且這種積分格式的自由面條件有相對(duì)較好的穩(wěn)定性，可以在很大程度上降低總體誤差，能防止出現(xiàn)結(jié)果發(fā)散的現(xiàn)象。參照文獻(xiàn)[13]中方法，在式(3)和(4)的基礎(chǔ)上，對(duì)線性自由面條件采用先積分后離散的處理方法，首先兩邊對(duì)t在(0,τ1)積分，然后再對(duì)τ1在(0,τ)中積分，得到

(7)

再由初始條件，上式可以寫(xiě)為

(8)

然后采用如下的離散模型將式(8)進(jìn)行離散：在空間上把自由面劃分成若干小平面單元，并假定每個(gè)單元中心點(diǎn)上的物理量為常數(shù)分布；在時(shí)間上以Δt為步距，等間隔前進(jìn)。采用梯形公式離散式(8)，可以得到

(p∈SF)

(9)

1.3 去奇異邊界元積分方法

采用邊界元方法的關(guān)鍵一步是積分方程的離散，本文采用分布源法，則流場(chǎng)中的速度勢(shì)可以表示為

(10)

因此式(2)和(9)中提到的邊界條件可改寫(xiě)成如下形式

φ(p,t)=φ0(p,t)

(11)

(12)

將式(11)和(12)代入式(10)，那么關(guān)于未知量σ0(q,t)的積分方程可表示為

(13)

(14)

對(duì)于式(13)、式(14)的積分方程，傳統(tǒng)的邊界元方法都是將源強(qiáng)直接分布在物體及自由表面上，也就是說(shuō)積分表面S1即為流域的邊界。這樣往往會(huì)出現(xiàn)由于格林函數(shù)1/rpq的分母趨于零所帶來(lái)的奇異積分的問(wèn)題。本文采用的方法為去奇異邊界元方法，即是將傳統(tǒng)作法中原本直接布置于流體計(jì)算域表面節(jié)點(diǎn)上的所有奇點(diǎn)，移至計(jì)算域外面，也就是將傳統(tǒng)作法中疊在一起的節(jié)點(diǎn)和奇點(diǎn)分開(kāi)了。在積分邊界SF和SW每個(gè)單元上分別布置強(qiáng)度為σF和σW的點(diǎn)源，SF和SW是實(shí)際流體域外距離流體域很近的積分表面，則：

(15)

將式(13)式(14)代入式(15)，則：

(p∈ΓF)

(16)

(17)

假設(shè)研究的問(wèn)題中總共有N個(gè)孤立的點(diǎn)源構(gòu)成，則式(16)和(17)可簡(jiǎn)化成如下的N個(gè)孤立點(diǎn)源的代數(shù)和的形式

(18)

(19)

當(dāng)采用上述離散方法對(duì)總共N個(gè)配置點(diǎn)列寫(xiě)方程后，可最終得到如下N×N的線性代數(shù)方程組

Aijσj=bi

(20)

式中：Aij為影響系數(shù)矩陣；bi為邊界條件所形成的確定矩陣；σj為待求的未知源強(qiáng)。

線性代數(shù)方程組(20)求解后即得出未知源強(qiáng)，由式(9)求得速度勢(shì)，可以通過(guò)下式求得水動(dòng)壓力p、波高η和水動(dòng)力F

(21)

1.4 基于DBIEM晃蕩程序的流程

基于DBIEM理論，以Fortran語(yǔ)言為工具，本文編寫(xiě)了一套能夠模擬任意尺寸任意形狀液艙晃蕩的程序，圖2是晃蕩程序的基本流程圖。

2 數(shù)值方法有效性驗(yàn)證

2.1 水平激勵(lì)下液艙晃蕩解析解

對(duì)于周期性激勵(lì)ue=Acosωt，其中ue是液艙激勵(lì)速度，A=bω速度幅值，b為位移幅值，ω為激勵(lì)頻率，F(xiàn)altinsen給出了求解速度勢(shì)φ的線性解析方法，通過(guò)下式可以很容易的將求得的φ值轉(zhuǎn)化為自由液面的位置分布[14]。(液艙水深h、長(zhǎng)度2a、寬度2w，且坐標(biāo)原點(diǎn)位置在自由面上液面的中心處。)

圖2 晃蕩程序基本流程Fig.2 The basic process of sloshing program

(22)

2.2 單向水平激勵(lì)下晃蕩數(shù)值解法驗(yàn)證

液艙的長(zhǎng)(L) 和寬(B)分別為1 m和0.5 m；液艙內(nèi)液體高度(h)為0.5 m。固有頻率ω0x和ω0y分別為5.316 rad/s和7.835 rad/s。為驗(yàn)證程序是否準(zhǔn)確，選取合適的網(wǎng)格大小、時(shí)間步長(zhǎng)以及去奇異距離(如表1所示)，模擬不同激勵(lì)頻率作用下液艙晃蕩問(wèn)題。

表1 主要計(jì)算參數(shù)Tab.1 The main calculation parameters

表1中，Nx為液艙長(zhǎng)度方向的劃分網(wǎng)格個(gè)數(shù)；Ny為液艙寬度方向的劃分網(wǎng)格個(gè)數(shù)；Nz為液艙高度方向的劃分網(wǎng)格個(gè)數(shù)；dt為時(shí)間步長(zhǎng)；T為激勵(lì)周期；ld為去奇異距離；A為激勵(lì)幅值。

考慮到液面條件為線性的，單向激勵(lì)幅值取為0.01h，h為液艙內(nèi)液面高度，入射頻率分別為0.5ω0(ω0為液艙在激勵(lì)方向上的固有頻率，此處即為液艙的縱向固有頻率)、0.9ω0、0.95ω0、0.999 9ω0、1.05ω0、1.1ω0。計(jì)算結(jié)果如圖3所示，圖中實(shí)線為數(shù)值解，虛線為解析解。

(a) ωP=0.5ω0(b) ωP=0.9ω0

(c) ωP=0.95ω0(d) ωP=0.999 9ω0

(e) ωP=1.05ω0(f) ωP=1.1ω0

圖3 不同激勵(lì)頻率下計(jì)算點(diǎn)波高時(shí)歷曲線
Fig.3 Time history of free surface elevation for different frequencies

可以看出，各激勵(lì)頻率下，數(shù)值解都能夠很好的與解析解吻合，誤差僅在1%以內(nèi)，這有效說(shuō)明了去奇異邊界元法在求解液艙晃蕩問(wèn)題的有效性。由以上各圖還可以看出計(jì)算點(diǎn)處波高時(shí)歷曲線隨著激勵(lì)頻率的不同其變化規(guī)律也發(fā)生了變化，外部激勵(lì)頻率遠(yuǎn)離液艙固有頻率時(shí)如圖3(a)所示，曲線變化較快，有一定規(guī)律性，波高最大值在0.006 m左右；當(dāng)入射頻率接近液艙固有頻率時(shí)，計(jì)算點(diǎn)處波峰首先隨時(shí)間不斷增加達(dá)到最大值后又開(kāi)始下降，并以此規(guī)律不斷循環(huán)，有很強(qiáng)的規(guī)律性，圖3(b)、(c)、(e)、(f)都顯示了這一特征，且波高最大值分別為0.06 m、0.11 m、0.15 m、0.1 m。不難看出越是接近液艙固有頻率波高最大值就越大，這是與共振密切相關(guān)的，理論上當(dāng)激勵(lì)頻率無(wú)限接近共振頻率時(shí)，如圖3(d)所示頻率為0.999 9ω0，隨著計(jì)算時(shí)間的推移該點(diǎn)處的波高將趨于無(wú)窮大，但由于程序算法和網(wǎng)格的限制，本文僅模擬了35 s，此時(shí)的最大幅值已經(jīng)約是圖3(a)的80倍。可見(jiàn)，為了保證所設(shè)計(jì)液艙的安全性，應(yīng)盡量避免出現(xiàn)共振現(xiàn)象，基于去奇異邊界元方法模擬液艙晃蕩問(wèn)題對(duì)實(shí)際工程問(wèn)題具有一定的指導(dǎo)作用。

3 不規(guī)則水平激勵(lì)作用下液艙晃蕩分析

在單向水平激勵(lì)下液艙晃蕩問(wèn)題研究的基礎(chǔ)上，本文也研究了不規(guī)則激勵(lì)作用下的液艙晃蕩，如圖4所示。

圖4 液艙晃蕩模型圖Fig.4 Model of sloshing tank

3.1 單向4成分疊加不規(guī)則激勵(lì)

首先采用四個(gè)不同于液艙固有頻率的單向激勵(lì)進(jìn)行研究，激勵(lì)參數(shù)選取如表2所示。

表2 四成分不規(guī)則激勵(lì)參數(shù)Tab.2 Parameters for irregular excitation with four

圖5為單向不規(guī)則激勵(lì)下點(diǎn)(-L/2,0)處波高時(shí)歷和相對(duì)體積誤差曲線，此處不規(guī)則激勵(lì)由四種不同單色激勵(lì)疊加而成。與單色激勵(lì)相比，不規(guī)則激勵(lì)下波高時(shí)歷曲線比較復(fù)雜，波高變化范圍為-0.004～0.004 m。相對(duì)體積誤差幾乎為0，符合質(zhì)量守恒物理規(guī)律。

(a) 波高時(shí)歷曲線(b) 相對(duì)體積誤差

圖5 不規(guī)則激勵(lì)下波高時(shí)歷和相對(duì)體積誤差曲線
Fig.5 Wave elevation history and relative volume error

for irregular excitation

圖6液艙波面圖，研究時(shí)間區(qū)間為13.252 0～16.345 9 s，約為一個(gè)波面運(yùn)動(dòng)周期，每隔Δt≈0.52的6個(gè)不同時(shí)刻波面圖，從圖中可以看到一個(gè)周期內(nèi)波面的基本變化情況。

3.2 單向512成分疊加不規(guī)則激勵(lì)

通?；诓ɡ俗V，一般200個(gè)左右規(guī)則波浪成分就已經(jīng)足夠模擬不規(guī)則波。為了更好的驗(yàn)證本文的數(shù)值方法，文中將采用512個(gè)成分的不規(guī)則波。本文將選用ITTC波浪譜研究隨機(jī)激勵(lì)作用下的液艙晃蕩情況，波浪譜如下

圖6 波面圖Fig.6 Free surface profiles

(23)

式中：Hs為有義波高；ωp為譜峰頻率。與隨機(jī)波浪組成相似，液艙運(yùn)動(dòng)速度和幅值寫(xiě)成

(24)

(25)

式中：Nω為激勵(lì)成分的個(gè)數(shù)，本文取512；φi為0～2π間隨機(jī)相位角，在程序中設(shè)定；ωi為激勵(lì)頻率。

選取有義波高Hs=0.005 h，譜峰頻率為ωp=0.9ω0，截?cái)鄥^(qū)域上限選取ωi≤5ω0波浪譜部分。

圖7為單向512成分不規(guī)則激勵(lì)下點(diǎn)(-L/2，0)處波高時(shí)歷和相對(duì)體積誤差曲線。其中，每個(gè)成分的波幅、頻率都不一樣，初始相位采用隨機(jī)數(shù)。模擬結(jié)果顯示，波高變化范圍為-0.015～0.015 m。由于模擬中使用了隨機(jī)相位，因此無(wú)法得到準(zhǔn)確的解析解，但根據(jù)相對(duì)體積誤差圖可以看出512成分不規(guī)則波激勵(lì)下模擬準(zhǔn)確度仍然很高，體積誤差接近于0。

(a) 波高時(shí)歷曲線(b) 相對(duì)體積誤差

圖7 不規(guī)則激勵(lì)下波高時(shí)歷和相對(duì)體積誤差曲線
Fig.7 Wave elevation history and relative volume error for irregular excitation

圖8是時(shí)間區(qū)間為13.252 0～16.345 9 s，每隔Δt≈0.52的6個(gè)不同時(shí)刻波面圖。

圖8 波面圖Fig.8 Free surface profiles

4 結(jié) 論

針對(duì)傳統(tǒng)邊界元方法求解波浪中結(jié)構(gòu)物運(yùn)動(dòng)響應(yīng)時(shí)會(huì)遇到奇點(diǎn)及交界面處速度分布精度較低的問(wèn)題，本文引入去奇異邊界元理論，利用Fortran語(yǔ)言開(kāi)發(fā)了一套可以針對(duì)不同形狀液艙晃蕩模擬的程序，程序中加入了各自由度運(yùn)動(dòng)耦合的模塊，可模擬液艙在六個(gè)自由度激勵(lì)聯(lián)合作用下內(nèi)部液體晃蕩各要素隨時(shí)間變化情況。本文主要對(duì)單自由度4成分不規(guī)則激勵(lì)、單自由度512成分不規(guī)則激勵(lì)下液艙內(nèi)部流體運(yùn)動(dòng)進(jìn)行了模擬，研究結(jié)果表明，文中方法可有效模擬隨機(jī)激勵(lì)下的液艙晃蕩問(wèn)題，且具有較高的精度。