浮體淺水波浪載荷數(shù)值計(jì)算方法研究

楊鵬，寇冠元，朱學(xué)康，劉成義，高長(zhǎng)華

1海南大學(xué)機(jī)電學(xué)院，海南海口570228

2武漢第二船舶設(shè)計(jì)研究所，湖北武漢430205

0 引 言

隨著國(guó)家對(duì)南海開(kāi)發(fā)和維護(hù)國(guó)家主權(quán)等需求的日益增強(qiáng)，需要在南海島礁附近布置大型浮式平臺(tái)，以供島礁開(kāi)發(fā)、綜合執(zhí)法補(bǔ)給和旅游等。島礁附近的水深一般較淺，位于淺水中浮體的水動(dòng)力運(yùn)動(dòng)和載荷響應(yīng)與深水中的存在較大差別，因此，需要研究淺水中浮體在波浪中的運(yùn)動(dòng)和載荷響應(yīng)。淺水中浮體與深水中浮體運(yùn)動(dòng)的最大差別和難點(diǎn)在于有限水深格林函數(shù)的準(zhǔn)確求解。只有準(zhǔn)確求解有限水深格林函數(shù)及其偏導(dǎo)數(shù)，才能得到浮體在波浪中的響應(yīng)。Li［1］、謝永和等［2］和劉日明等［3］提出了數(shù)值積分方法，用來(lái)計(jì)算有限水深的格林函數(shù)及其導(dǎo)數(shù)。有限水深格林函數(shù)求解方法分為積分形式［4］和級(jí)數(shù)形式［5］2 種。其中積分形式的計(jì)算精度高，遠(yuǎn)場(chǎng)和近場(chǎng)均適用，但計(jì)算效率低；級(jí)數(shù)形式計(jì)算效率高，但在近場(chǎng)附近很難收斂，同時(shí)在遠(yuǎn)方輻射半徑R=0處存在奇點(diǎn)。所以在計(jì)算有限水深格林函數(shù)及其導(dǎo)數(shù)時(shí)一般在近場(chǎng)采用積分形式，遠(yuǎn)場(chǎng)采用級(jí)數(shù)形式。

積分形式的解是主值積分且存在奇異性，這是有限水深格林函數(shù)及其導(dǎo)數(shù)求解的難點(diǎn)所在。傳統(tǒng)的有限水深格林函數(shù)積分計(jì)算公式采用多點(diǎn)Gauss-Laguerre公式直接進(jìn)行積分或多項(xiàng)式逼近。直接積分法一般需要取64個(gè)高斯積分點(diǎn)［1］方能滿足精度要求，但該方法耗時(shí)長(zhǎng)且誤差較大，尤其是在高頻率處計(jì)算失真。多項(xiàng)式逼近法（例如法國(guó)船級(jí)社的Hydrostar）需要計(jì)算大量的數(shù)據(jù)并選擇適當(dāng)?shù)谋平鼌^(qū)間，難以實(shí)施。本文擬通過(guò)推導(dǎo)和數(shù)值實(shí)驗(yàn)提出精確計(jì)算有限水深格林函數(shù)及其偏導(dǎo)數(shù)的方法，以解決高頻失真問(wèn)題。Newman［6］認(rèn)為在半徑和水深比R/H＞0.5時(shí)采用適當(dāng)項(xiàng)數(shù)的級(jí)數(shù)解可以得到理想的精度，但其并沒(méi)有給出小R/H的簡(jiǎn)化解。為了提高計(jì)算效率，本文擬給出對(duì)稱性的處理方法和簡(jiǎn)化的級(jí)數(shù)求解公式的實(shí)施方案，并編制計(jì)算程序?qū)⒈疚姆椒ㄅc商業(yè)軟件計(jì)算結(jié)果進(jìn)行比較分析，以驗(yàn)證本文方法的正確性和可行性。

1 淺水復(fù)雜格林函數(shù)法

1.1 入射勢(shì)求解

設(shè)一階波面升高表達(dá)式為

式中：A，K，ω和β分別為波幅、波數(shù)、波浪自然頻率和浪向，例如，頂浪為180°。坐標(biāo)系定義：x指向船艏，y指向船左舷，z垂直于靜水面向上，xyz符合右手法則。

那么，入射勢(shì)ΦI的一階表達(dá)式為

式中：?I為入射速度勢(shì)的幅值；H為水深；g為重力加速度。

1.2 輻射勢(shì)求解

無(wú)航速條件下，有限水深（均勻海底）的輻射勢(shì)求解條件如下。

在流域內(nèi)，

在z=0上，

在物面不可穿透條件下，

在水底不可穿透條件下，

遠(yuǎn)方輻射條件下，

為了滿足自由面邊界條件、無(wú)窮遠(yuǎn)輻射條件和水底不可穿透條件，有限水深格林函數(shù)G一般采用如下2種形式。

1）積分形式［4］：

式中：P.V.為取主值積分；P(x，y，z)和Q(ξ，η，ζ)分別為場(chǎng)點(diǎn)和源點(diǎn)；k為積分變量；關(guān)于池底的鏡像；J0()為第1類(lèi)零階Bessel函數(shù)。

2）級(jí)數(shù)形式［5］：

式中：kn為方程kntanknH+ν=0的正實(shí)根；為第1類(lèi)零階Hankel函數(shù)；K0()為第2類(lèi)零階修正的Bessel函數(shù)。

那么流域內(nèi)各點(diǎn)的速度勢(shì)為

決定源強(qiáng)的積分表達(dá)式為

那么，附加質(zhì)量μjk和附加阻尼λjk由下式求得：

式中：ρ為流體密度；?k為輻射勢(shì)。

1.3 繞射勢(shì)求解

繞射勢(shì)?D的求解方程為

在流域內(nèi)，

在z=0上，

在物面上，

在水底，

遠(yuǎn)方輻射條件下，

得到繞射勢(shì)的計(jì)算公式為

因此，可以不必求解繞射勢(shì)，直接由波浪入射勢(shì)?I和輻射勢(shì)?j求得波浪力。輻射勢(shì)在前面已求得，第j階波浪力（入射勢(shì)+繞射勢(shì)）的表達(dá)式為

2 淺水波浪載荷數(shù)值計(jì)算方法

在淺水輻射勢(shì)和繞射勢(shì)的求解過(guò)程中，首先需要對(duì)格林函數(shù)及其導(dǎo)數(shù)進(jìn)行求解。

2.1 格林函數(shù)的直接積分法

為了書(shū)寫(xiě)方便，式（4）和式（5）中的格林函數(shù)G及其偏導(dǎo)數(shù)的主值積分部分G0如下：

式中，J1()為第1類(lèi)一階Bessel函數(shù)。

式（14）～式（16）中不帶奇異性的積分通過(guò)借鑒下式得到：

式中，a和b為廣義變量。

為書(shū)寫(xiě)方便，引入

結(jié)合式（18）和式（19），可以得到

其中第1項(xiàng)積分無(wú)奇異，可以采用Gauss-Laguerre積分方法直接計(jì)算，第2項(xiàng)積分可以通過(guò)無(wú)限水深的格林函數(shù)計(jì)算公式得到，那么

式中：wj為Gauss-Laguerre積分的權(quán)重系數(shù)；，為無(wú)限水深格林函數(shù)的積分公式。

2.2 格林函數(shù)的級(jí)數(shù)解

根據(jù)前文所述，淺水格林函數(shù)的求解公式可以采用如下級(jí)數(shù)解形式：

針對(duì)式（23）中的高階級(jí)數(shù)項(xiàng)，有

式（24）表明，式（23）的收斂速度取決于R H，當(dāng)R/H=0時(shí)，式（23）將不收斂：

1）當(dāng)R/H≥0.50時(shí)，如果式（23）的n=10，

2）當(dāng)R/H≥0.25時(shí)，如果式（23）的n=20，

3）當(dāng)R/H≥0.10時(shí)，如果式（23）的n=50，

因此，當(dāng)R/H≥0.25時(shí)，可以選取適當(dāng)項(xiàng)數(shù)的級(jí)數(shù)解來(lái)得到較為精確的值。Newman［6］證明了當(dāng)R/H≥0.5時(shí)，采用級(jí)數(shù)解可以得到精度較高的近似值，級(jí)數(shù)項(xiàng)數(shù)一般取6H/R的整數(shù)。另外，當(dāng)當(dāng)knR≥8.0時(shí)，經(jīng)過(guò)上述簡(jiǎn)化，相比于直接積分，可以較快速地計(jì)算出格林函數(shù)及其導(dǎo)數(shù)的值。

2.3 具有對(duì)稱面的計(jì)算處理

對(duì)于具有對(duì)稱面的物體，在計(jì)算時(shí)可以節(jié)省大量的內(nèi)存和時(shí)間。在求解源強(qiáng)時(shí)，其系數(shù)矩陣為滿秩矩陣，只能使用計(jì)算效率較低的高斯消元法和LU分解法等，其時(shí)間度為O(n3)，所以利用矩陣的對(duì)稱性求解線性方程組能夠極大地提高求解速度。當(dāng)物體具有1個(gè)對(duì)稱面時(shí)，方程可以轉(zhuǎn)化為2個(gè)n/2階的線性方程組問(wèn)題，線性方程組的求解時(shí)間降為原來(lái)的1/4；當(dāng)物體具有2個(gè)對(duì)稱面時(shí)，方程可以轉(zhuǎn)化為4個(gè)n/4階的線性方程組問(wèn)題，線性方程組的求解時(shí)間降為原來(lái)的1/16。下面，以具有2個(gè)對(duì)稱面的浮體計(jì)算源強(qiáng)時(shí)的推導(dǎo)過(guò)程為例，具有1個(gè)對(duì)稱面的情況與此類(lèi)似。

為便于討論，將物面上的點(diǎn)分成4個(gè)區(qū)域（圖1）。

圖1 浮體示意圖Fig.1 Schematic diagram of floating body

對(duì)于這類(lèi)具有2個(gè)對(duì)稱面（yz和xz平面）的物體，存在如下關(guān)系式：

式中：i和j分別為場(chǎng)點(diǎn)和源點(diǎn)的面元序號(hào)；角標(biāo)?為速度勢(shì)。

設(shè)在區(qū)域I上有N個(gè)單元，那么整個(gè)物面的平均濕表面上共有4N個(gè)單元。從式（29）～式（33）可以看出，并不需要計(jì)算整個(gè)物面的速度勢(shì)矩陣G0，G0x，G0y，G0z和G0n，只需對(duì)每個(gè)矩陣計(jì)算4個(gè)N×N矩陣。對(duì)于具有2個(gè)對(duì)稱面（yz和xz平面）的濕表面來(lái)說(shuō)，可采用式（33）求解源強(qiáng)的系數(shù)矩陣。該矩陣為循環(huán)矩陣，可將原方程組的求解化為4個(gè)N階矩陣的方程組來(lái)求解［8］，下面介紹具體的計(jì)算方法。

式中，{ni}為第i個(gè)面元的法向。

定義4N階矩陣P4N和P2N為

式中，E2N和EN分別為2N階和N階單位矩陣。P4N可將G0n對(duì)角化，即

P2N可將式（36）對(duì)角化為

源強(qiáng)的線性求解方程組變換為

將式（37）代入式（38），可得

3 波浪載荷計(jì)算方法

3.1 浮體表面的波動(dòng)壓力

應(yīng)用面元法分別求解得到輻射勢(shì)?j（j=1，2，3，4，5，6）之后，便可通過(guò)規(guī)則波中的浮體運(yùn)動(dòng)方程獲得穩(wěn)態(tài)的運(yùn)動(dòng)響應(yīng)ηj（j=1，2，3，4，5，6），那么非定常輻射勢(shì)?R和繞射勢(shì)?D便隨之確定。浮體表面的波動(dòng)壓力除了由波浪入射勢(shì)、輻射勢(shì)和繞射勢(shì)引起的脈動(dòng)壓力，還將包含由浮體運(yùn)動(dòng)引起的浮體表面靜水壓力變化ps項(xiàng)，總的脈動(dòng)壓力為［9］

式中，η3，η4和η5分別為浮體的垂蕩、橫搖和縱搖運(yùn)動(dòng)。

3.2 浮體橫剖面的波浪載荷

確定了規(guī)則波中的浮體運(yùn)動(dòng)響應(yīng)之后，便可應(yīng)用達(dá)朗貝爾原理計(jì)算浮體橫剖面的力和力矩，包括波浪誘導(dǎo)的垂向彎矩與水平的剪力和彎矩以及軸力和扭矩。

3.2.1 部分長(zhǎng)度浮體的剛體慣性力載荷

以浮體x軸與橫剖面的交點(diǎn)為坐標(biāo)原點(diǎn)O，選取單位長(zhǎng)度的浮體，設(shè)其質(zhì)量為μ(x)，質(zhì)心坐標(biāo)為(x，y，z) ，其慣性質(zhì)量矩陣為

式中：Ixx，Iyy，Izz，Ixy，Ixz和Iyz為關(guān)于坐標(biāo)x，y，z的慣性矩（積）；Sx≡μx，Sy≡μy，Sz≡μz，均為剛體質(zhì)量關(guān)于坐標(biāo)平面的靜矩。

任意橫剖面x處至船艏xf之間的質(zhì)量對(duì)橫剖面的剛體慣性力載荷為

式中，{}為浮體六自由度運(yùn)動(dòng)加速度。

3.2.2 浮體橫剖面的力和力矩

作用于任意橫剖面x處至船艏xf之間部分長(zhǎng)度浮體上的真實(shí)流體載荷與橫剖面上的力相互平衡，那么橫剖面上的載荷為

式中：Sx，M，zb和zg分別為橫剖面x處至船艏xf之間部分長(zhǎng)度浮體的濕表面、質(zhì)量、浮心和重心的z坐標(biāo)；p為流體壓力；η為浮體六自由度運(yùn)動(dòng)幅值。

4 淺水?dāng)?shù)值計(jì)算結(jié)果對(duì)比分析

4.1 箱型浮體計(jì)算結(jié)果

針對(duì)ISSC（2006）［10］標(biāo)準(zhǔn)浮箱進(jìn)行了運(yùn)動(dòng)傳遞函數(shù)的計(jì)算，其主尺度和計(jì)算結(jié)果如表1和圖2所示。

本次計(jì)算的1/4模型網(wǎng)格數(shù)為155個(gè)。從圖2可以看出，針對(duì)無(wú)限水深和20 m水深，本文程序計(jì)算結(jié)果與采用商業(yè)軟件AQWA獲得的計(jì)算結(jié)果十分吻合，驗(yàn)證了本文計(jì)算方法的正確性。

表1 浮箱主尺度Table 1 Main dimensions of box

圖2 箱型浮體運(yùn)動(dòng)和載荷傳遞函數(shù)（0°浪向）Fig.2 Transfer function of motion and loads of box（0°wave direction）

4.2 大型散貨船計(jì)算結(jié)果

針對(duì)一艘20.5×104t大型散貨船，計(jì)算了零航速（壓載）工況下縱搖和垂蕩傳遞函數(shù)，其主尺度和計(jì)算結(jié)果如表2和圖3所示。

本次計(jì)算的半模型網(wǎng)格數(shù)為440個(gè)。從圖3可以看出，針對(duì)無(wú)限水深和40 m水深，本文程序計(jì)算結(jié)果與采用商業(yè)軟件Wadam獲得的計(jì)算結(jié)果十分吻合，驗(yàn)證了本文計(jì)算程序和方法的正確性。

表2 散貨船主尺度Table 2 Main dimensions of bulk carrier

圖3 散貨船運(yùn)動(dòng)傳遞函數(shù)Fig.3 Transfer function of motion of bulk carrier

5 結(jié) 語(yǔ)

基于勢(shì)流理論和邊界元方法，給出了速度勢(shì)求解方程，同時(shí)利用循環(huán)矩陣原理對(duì)復(fù)雜格林函數(shù)進(jìn)行了嚴(yán)謹(jǐn)?shù)那蠼馔茖?dǎo)和對(duì)稱性處理研究，從而建立了浮體在淺水中運(yùn)動(dòng)和波浪載荷的完整計(jì)算方法。本文給出的改進(jìn)方法和思路可以精確計(jì)算格林函數(shù)及其導(dǎo)數(shù)，同時(shí)通過(guò)研究有限水深格林函數(shù)的級(jí)數(shù)解形式，給出了簡(jiǎn)化的快速計(jì)算方法。針對(duì)具有對(duì)稱面的浮體，引入循環(huán)矩陣進(jìn)行嚴(yán)格數(shù)學(xué)推導(dǎo)，發(fā)現(xiàn)利用對(duì)稱性得到簡(jiǎn)化后的輻射勢(shì)求解方法可節(jié)省大量計(jì)算時(shí)間。

通過(guò)自主開(kāi)發(fā)程序進(jìn)行數(shù)值計(jì)算，并與成熟商業(yè)軟件結(jié)果進(jìn)行對(duì)比分析，驗(yàn)證了本文給出的有限水深格林函數(shù)計(jì)算方法和相應(yīng)計(jì)算程序的正確性。本文建立的計(jì)算方法對(duì)淺水中的浮體運(yùn)動(dòng)和波浪載荷評(píng)估具有重要意義。