王慶豐，徐 剛，王樹(shù)齊，田 超

（1.江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江212003；2.中國(guó)船舶科學(xué)研究中心，江蘇 無(wú)錫214082）

0 引言

頻域理論的應(yīng)用可以說(shuō)已經(jīng)相當(dāng)成熟，但是頻域法存在一個(gè)本質(zhì)上的限制，即它通常只適用于周期穩(wěn)態(tài)問(wèn)題，對(duì)于瞬變或者強(qiáng)非線性問(wèn)題往往無(wú)法處理。但是時(shí)域理論是頻域理論的超集，是一種時(shí)間項(xiàng)不可分離的方法。理論上可以解決完全非線性以及物體任意運(yùn)動(dòng)的問(wèn)題，具有相當(dāng)大的自由度。本文采用基于簡(jiǎn)單格林函數(shù)的時(shí)域邊界元方法，對(duì)于傳統(tǒng)邊界元方法，奇異點(diǎn)的處理一直是難點(diǎn)。本文使用的去奇異邊界元法（Desingularized Boundary Integral Equation Method）將傳統(tǒng)方法中直接布置于流體計(jì)算域表面節(jié)點(diǎn)上的所有奇點(diǎn)移至計(jì)算域外部，很好地解決了奇異積分問(wèn)題。為了能夠更好地模擬液艙內(nèi)的運(yùn)動(dòng)，本文利用FORTRAN 研發(fā)了一套可以模擬任意尺度液艙晃蕩的程序。

目前，基于去奇異邊界元法對(duì)勢(shì)流問(wèn)題的研究已越來(lái)越多。Jensen 等[1]利用去奇異方法來(lái)分析在穩(wěn)定入射波情況下求解阻力的問(wèn)題；Cao，Schultz 和Beck[2]和Lee[3]基于時(shí)間步進(jìn)技術(shù)和去奇點(diǎn)邊界積分方程，結(jié)合二者來(lái)模擬簡(jiǎn)單潛體及浮體運(yùn)動(dòng)的完全非線性效應(yīng)問(wèn)題；李宗翰等[4]使用去奇異積分方程方法仿真了任意三維浮體的完全非線性流場(chǎng)；Scorpio 和Beck[5]，王大國(guó)等[6]對(duì)完全非線性數(shù)值水池進(jìn)行了探索；Kara 等[7]利用去奇異邊界元方法研究了完全非線性波物相互作用三維時(shí)域問(wèn)題；Zhang 等[8]基于去奇異方法研究了雙體船的耐波性；徐剛等[9]進(jìn)行了隨機(jī)波浪下三維液艙的全非線性數(shù)值分析。

1 去奇異邊界元理論

1.1 控制方程及初始條件

根據(jù)勢(shì)流理論，基于流體運(yùn)動(dòng)的連續(xù)性條件和無(wú)旋條件，速度勢(shì)φ 在流體域D 中滿(mǎn)足Laplace 方程，定解條件如下：

給定初始條件如下：

1.2 自由面條件差分格式

自由表面上的速度勢(shì)φ 時(shí)間步進(jìn)求解一般采用有限差分法。這種處理方法的缺點(diǎn)是相對(duì)比較復(fù)雜，而且在此過(guò)程中容易出現(xiàn)數(shù)值發(fā)散。本文采用一種新的積分格式的自由面條件（Integral form of free surface boundary condition[10]，IFBC）。通過(guò)此積分格式，可以由（3）式和（4）式得到自由面上一階速度勢(shì)的表達(dá)式，見(jiàn)（7）式。這種積分格式的自由面條件穩(wěn)定性高，誤差小，不會(huì)出現(xiàn)結(jié)果發(fā)散的現(xiàn)象。

對(duì)線性自由面條件采用先積分后離散的處理方法，在等式兩邊對(duì)t 在（0，τ1）重積分，然后再對(duì)τ1在（0，τ）中積分，可得到：

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

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

1.3 去奇異邊界元積分方法

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

式中：p（x，y，z）為場(chǎng)點(diǎn)，即流體域內(nèi)的動(dòng)點(diǎn)；q（ξ，η，ζ）為源點(diǎn)，即點(diǎn)源所在固定點(diǎn)；rpq為p 和q 兩點(diǎn)間距離，即；S1為流體域上積分表面。

因此上節(jié)提到的邊界條件可改寫(xiě)成如下形式：

將上述二式代入（11）式，那么關(guān)于未知量σ0（q，t）的積分方程可表示為：

對(duì)于（14）式和（15）式的積分方程，傳統(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í)際流體域外距離流體域很近的積分表面，則：

將（14）-（15）式代入上式，則：

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

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

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

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

1.4 基于DBIEM 晃蕩程序的流程

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

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

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

2.1 線性晃蕩解析解

2.1.1 水平方向激勵(lì)液艙晃蕩

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

2.1.2 橫蕩縱蕩耦合線性激勵(lì)液艙晃蕩

對(duì)于周期性振動(dòng)ue= Acosωt，我們將振動(dòng)速度ue投影到x 與y 方向并將此類(lèi)問(wèn)題當(dāng)作橫蕩與縱蕩耦合問(wèn)題考慮。將一維線性激勵(lì)下x 與y 方向的Faitinsen 解析解進(jìn)行結(jié)合，即得到橫蕩與縱蕩耦合線性激勵(lì)下三維液艙液面升高解析解（液艙水深h、長(zhǎng)度2a、寬度2w，且原點(diǎn)位置在液面中心處）。

2.2 數(shù)值解與解析解對(duì)比

在數(shù)值模擬中，取坐標(biāo)系O-XYZ，OXY 平面與靜水面重合，Z 軸垂直向上。液艙的長(zhǎng)和寬分別為L(zhǎng)和B，h 表示液艙內(nèi)液體高度。表1所示為液艙的主要參數(shù)。

表1 液艙主要模擬參數(shù)Tab.1 Main simulation parameters of tank

為驗(yàn)證程序是否準(zhǔn)確，本文在網(wǎng)格收斂性分析的基礎(chǔ)上，選取合適的網(wǎng)格大小、時(shí)間步長(zhǎng)以及去奇異距離（如表2所示），模擬不同激勵(lì)頻率作用下液艙晃蕩問(wèn)題。

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

表2 中Nx 為液艙長(zhǎng)度方向的劃分網(wǎng)格個(gè)數(shù)；Ny 為液艙寬度方向的劃分網(wǎng)格個(gè)數(shù)；Nz 為液艙高度方向的劃分網(wǎng)格個(gè)數(shù)；dt 為時(shí)間步長(zhǎng)；T 為激勵(lì)周期；ld為去奇異距離，相當(dāng)于一參數(shù)，取值范圍為0～1；A 為激勵(lì)幅值，單位為m。

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

考慮到微幅運(yùn)動(dòng)，單向激勵(lì)幅值取為0.01h，h 為液艙內(nèi)液面高度，入射頻率分別為0.5ω0（ω0為液艙在激勵(lì)方向上的固有頻率，此處即為液艙的縱向固有頻率），0.999 9ω0和1.1ω0。計(jì)算結(jié)果如圖2所示，圖中點(diǎn)劃線為數(shù)值解，實(shí)線為解析解。

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

從圖2 可以看出，各液艙激勵(lì)頻率下，數(shù)值解都能夠很好地與解析解吻合，誤差僅在1%以?xún)?nèi)，驗(yàn)證了去奇異邊界元法求解液艙晃蕩問(wèn)題的有效性。由以上各圖還可以看出計(jì)算點(diǎn)處波高時(shí)歷曲線隨著激勵(lì)頻率的不同其變化規(guī)律也出現(xiàn)變化，激勵(lì)頻率遠(yuǎn)離液艙固有頻率時(shí)如圖2（a）所示，曲線變化較快，有一定規(guī)律性，波高最大值在0.006 m 左右；當(dāng)激勵(lì)頻率接近液艙固有頻率時(shí)，計(jì)算點(diǎn)處波峰首先隨時(shí)間不斷增加達(dá)到最大值后又開(kāi)始下降，并以此規(guī)律不斷循環(huán)，有很強(qiáng)的規(guī)律性，圖2（c）顯示了這一特征。不難看出越接近液艙固有頻率波高峰值就越大，如圖2（b）所示。

2.2.2 橫蕩縱蕩耦合激勵(lì)作用下液艙晃蕩數(shù)值分析

針對(duì)單向微幅激勵(lì)下液艙內(nèi)流體運(yùn)動(dòng)模擬，程序已經(jīng)可以達(dá)到很高的精度。當(dāng)液艙受到雙向耦合激勵(lì)時(shí)流體運(yùn)動(dòng)模式會(huì)有很大的不同，程序中也要加入相應(yīng)的耦合模塊。因此，對(duì)于耦合激勵(lì)時(shí)液艙內(nèi)流體運(yùn)動(dòng)也要進(jìn)行相應(yīng)的對(duì)比分析。圖3-5給出了不同耦合激勵(lì)頻率組合下，液艙四個(gè)角點(diǎn)處波高時(shí)歷曲線。

圖3 ωpx=0.5ω0x，ωpy=0.5ω0y 時(shí)液面四角點(diǎn)的波高時(shí)歷曲線Fig.3 Wave elevation history at four corners when ωpx=0.5ω0x,ωpy=0.5ω0y

圖4 ωpx=0.9ω0x，ωpy=0.5ω0y 時(shí)液面四角點(diǎn)的波高時(shí)歷曲線Fig.4 Wave elevation history at four corners when ωpx=0.9ω0x,ωpy=0.5ω0y

圖5 ωpx=0.9ω0x，ωpy=0.9ω0y 時(shí)液面四角點(diǎn)的波高時(shí)歷曲線Fig.5 Wave elevation history at four corners when ωpx=0.9ω0x,ωpy=0.9ω0y

圖中ωpx，ωpy分別為沿x 方向與y 方向的激勵(lì)頻率；ω0x，ω0y分別為液艙沿x 方向與y 方向的固有頻率。與單一激勵(lì)相似，激勵(lì)頻率越接近固有頻率，液艙內(nèi)流體晃蕩越劇烈，ωpx=0.5ω0x，ωpy=0.5ω0y時(shí)晃蕩最大值約為0.013 m；ωpx=0.9ω0x，ωpy=0.9ω0y已經(jīng)達(dá)到0.1 m。ωpx=0.9ω0x，ωpy=0.5ω0y時(shí)時(shí)歷曲線表現(xiàn)出了很強(qiáng)的規(guī)律性，這與單一激勵(lì)結(jié)果相似；但ωpx=0.9ω0x，ωpy=0.9ω0y時(shí)卻比較復(fù)雜，這主要是由于當(dāng)耦合激勵(lì)頻率接近固有頻率時(shí)會(huì)產(chǎn)生共振并有強(qiáng)烈的相互干擾。

3 結(jié) 語(yǔ)

針對(duì)邊界元方法求解波浪中結(jié)構(gòu)物運(yùn)動(dòng)響應(yīng)，尤其是在處理非線性問(wèn)題時(shí)會(huì)遇到奇異性問(wèn)題，本文采用了去奇異邊界元方法，同時(shí)基于這一理論研發(fā)了一套可以針對(duì)不同模型不同尺寸液艙進(jìn)行晃蕩模擬的程序，程序中加入了六自由度運(yùn)動(dòng)耦合的模塊，可模擬液艙在六個(gè)自由度耦合激勵(lì)作用下內(nèi)部液體晃蕩各要素隨時(shí)間變化情況，文中僅針對(duì)單一激勵(lì)和橫蕩縱蕩耦合激勵(lì)下液艙晃蕩問(wèn)題的數(shù)值模擬，并將結(jié)果與解析解進(jìn)行對(duì)比。結(jié)果表明，數(shù)值解和解析解吻合很好，驗(yàn)證了文中方法的有效性及對(duì)此問(wèn)題的精度，可以進(jìn)一步拓展到二階或者全線性問(wèn)題。