楊鵬，寇冠元，朱學康，劉成義，高長華

1海南大學機電學院，海南?？?70228

2武漢第二船舶設計研究所，湖北武漢430205

0 引 言

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

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

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

1.1 入射勢求解

設一階波面升高表達式為

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

那么，入射勢ΦI的一階表達式為

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

1.2 輻射勢求解

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

在流域內，

在z=0上，

在物面不可穿透條件下，

在水底不可穿透條件下，

遠方輻射條件下，

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

1）積分形式［4］：

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

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

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

那么流域內各點的速度勢為

決定源強的積分表達式為

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

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

1.3 繞射勢求解

繞射勢?D的求解方程為

在流域內，

在z=0上，

在物面上，

在水底，

遠方輻射條件下，

得到繞射勢的計算公式為

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

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

在淺水輻射勢和繞射勢的求解過程中，首先需要對格林函數(shù)及其導數(shù)進行求解。

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

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

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

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

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

為書寫方便，引入

結合式（18）和式（19），可以得到

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

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

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

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

針對式（23）中的高階級數(shù)項，有

式（24）表明，式（23）的收斂速度取決于R H，當R/H=0時，式（23）將不收斂：

1）當R/H≥0.50時，如果式（23）的n=10，

2）當R/H≥0.25時，如果式（23）的n=20，

3）當R/H≥0.10時，如果式（23）的n=50，

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

2.3 具有對稱面的計算處理

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

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

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

對于這類具有2個對稱面（yz和xz平面）的物體，存在如下關系式：

式中：i和j分別為場點和源點的面元序號；角標?為速度勢。

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

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

定義4N階矩陣P4N和P2N為

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

P2N可將式（36）對角化為

源強的線性求解方程組變換為

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

3 波浪載荷計算方法

3.1 浮體表面的波動壓力

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

式中，η3，η4和η5分別為浮體的垂蕩、橫搖和縱搖運動。

3.2 浮體橫剖面的波浪載荷

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

3.2.1 部分長度浮體的剛體慣性力載荷

以浮體x軸與橫剖面的交點為坐標原點O，選取單位長度的浮體，設其質量為μ(x)，質心坐標為(x，y，z) ，其慣性質量矩陣為

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

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

式中，{}為浮體六自由度運動加速度。

3.2.2 浮體橫剖面的力和力矩

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

式中：Sx，M，zb和zg分別為橫剖面x處至船艏xf之間部分長度浮體的濕表面、質量、浮心和重心的z坐標；p為流體壓力；η為浮體六自由度運動幅值。

4 淺水數(shù)值計算結果對比分析

4.1 箱型浮體計算結果

針對ISSC（2006）［10］標準浮箱進行了運動傳遞函數(shù)的計算，其主尺度和計算結果如表1和圖2所示。

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

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

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

4.2 大型散貨船計算結果

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

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

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

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

5 結 語

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

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