限制水域船體水動力計算分析與自由面數值阻尼的應用

王偉飛 彭亞康 濮 駿 趙曉斌 楊震峰

（中國船舶及海洋工程設計研究院 上海200011）

引 言

浮體在限制水域的運動計算具有很多的工程用途（如船在河流或航道中的運動等），目前商業軟件HydroStar 提供此項功能，但在此功能下無法引入自由面數值阻尼。池壁格林函數的計算一直是浮體在限制水域運動計算的一個重點問題，Newman［1］討論了浮體在航道和水池中以定常速度航行下的水動力計算問題，其采用了無界流假定，池壁格林函數由1/r關于池壁的無窮個鏡像疊加得到，此方法無法應用于自由面問題。C.M. Linton［2］導出了一種新的滿足自由面條件、基于有限水深池壁格林函數的表達形式，但其表達方式與目前常用的開敞水域格林函數的形式差別很大，需要重新編制算法，且其收斂性依賴部分系數，實施起來并不方便。目前較為流行的做法為通過將無窮個開敞水域自由面格林函數關于池壁進行鏡像疊加［3-5］，該方法可充分利用現有的開敞水域格林函數計算法，但需要計算一個收斂較慢的無窮級數求和問題。對于無窮級數求和，可以采用一些加速的辦法［6］，但也無法滿足本文中計算的要求。Xia Jinzhu［3-4］和Chen-Xiaobo［5］采取了對池壁格林函數的分區域計算方法克服了收斂慢的問題，本文在此基礎上采用漸近法，進一步分析了其奇異性，從中分離出了奇異項，并對剩余項進行切比雪夫逼近，實現了池壁格林函數的高效和高精度計算。

線性勢流理論水動力計算無法計入粘性和非線性效應，導致在計算部分問題時會出現數值共振問題從而引起計算結果失真。為處理這一問題，可在勢流方法中引入人工數值耗散項［7］，即數值阻尼來解決。本文將自由面數值阻尼引入限制水域的浮體運動計算，并通過自由面積分的解析解較好地解決了這一問題。本文的所有計算均通過Fortran編程計算求解。

1 控制方程和定解條件

根據三維線性勢流理論，在一寬度為b，深為h 的航道中，浮體以頻率ω和時間因子e-iωt作周期性運動；記源點坐標為（ξ，η，ζ），場點坐標為（x，y，z），池壁格林函數G滿足的定解條件如下［3-4］：

式中：g 為重力加速度。滿足上述定解條件的解為：

式（5）等號右側為開敞有限水深自由面格林函數的無窮級數（m為整數），其表達式為：

Gm的級數表達式為：

（6）式和（7）式中：J0為0 階第一類貝塞爾函數；

K0為0 階修正的貝塞爾函數；

H0為0 階漢克爾函數。

2 自由面池壁格林函數計算

由式（5）可知，池壁格林函數可通過有限水深自由面格林函數展開的無窮級數得到。但是由于式（5）級數收斂非常慢，直接進行求和計算的成本非常高，實際工程應用中幾乎是不可行的，所以需要找到高效的數值方法進行計算。

為提高計算效率，并保證計算精度，將式（5）分為3 部分［3-5］，即，

式中：包括近場部分GN，中場部分GM和遠場部分GF。近場采用級數直接求和的方法：

中場部分利用式（7），并注意到修正的貝塞爾函數漸近表達式為

式（16）中項數M0和式（17）中項數M1的取法在文獻［3-4］中已有詳細論述，在此不再贅述。利用Graf’s 加法公式［6］和漢克爾函數的漸近展開近似［6］，可得 ：

式中：

er為剩余偏差項；

B=k0b，表示無因次化的航道寬度，系數

3 數值方法與切比雪夫多項式逼近

Newman［8］給出了近場部分GN的詳細計算方法，其中為了提高計算效率，對近場部分進行了切比雪夫多項式逼近，中場部分GM計算較簡單，可通過級數求和以及有理分式逼近或直接采用現有的函數庫來求解，以上計算方法無需贅述。本節主要介紹遠場GF的高效計算方法以及對其進行切比雪夫多項式逼近方法。

其中，剩余項的漸近表達式LerExtract（z，s，a）的奇異性與Φ（z，s，a）一致，所以原函數減去奇異項之后的剩余項可通過多變量切比雪夫多項式來逼近：

式中：Cjk為逼近多項式的系數；Tj和Tk分別是關于z和a的切比雪夫多項式；截斷項數J和K根據不同的區域來選取，通常可通過試算確定。

為更好說明式（31）， 現取s=1/2，a=2/3，b∈（0，π/2）分別對Φ（z，s，a）、LerExtract（z，s，a）、Φ（z，s，a）-LerExtract（z，s，a）作圖，對應下頁圖1-圖3。可以看出：Φ（z，s，a）和LerExtract（z，s，a）在b= 0 處存在奇異性，剩余項不再存在奇異性，且為光滑函數，可用切比雪夫多項式高效地逼近。

圖1 Φ(2bi，1/2，3/2)，b∈(0，π/2)曲線

圖2 LerExtract(2bi，1/2，3/2)，b∈(0，π/2)曲線

圖3 Φ(2bi，1/2，3/2)-LerExtract(2bi，1/2，3/2)，b∈(0，π/2)曲線

4 自由面人工數值阻尼施加

由于勢流理論無法計入粘性效應等因素，采用池壁格林函數求解浮體運動時在某些特定的頻率會產生數值共振，造成計算結果與實際情況有較大的差異。為了較好地處理數值計算所帶來的奇異性問題，可通過引入人工數值耗散項來解決。通常的做法是在自由面條件中引入數值阻尼ε［7］。

當場點在自由面上時，積分方程為：

式中：σ為源強；G和Gn分別用池壁格林函數及其法向導數代入。由于池壁格林函數在自由面上存在奇異性，為保證計算精度，需要小心處理。可通過去奇點法進行數值積分計算，即積分核中先減去一奇異項，此奇異項的積分可給出解析表達式。阻尼系數ε取一小值，通常可通過試算和試驗方法確定。

5 算例與比較分析

為驗證池壁格林函數計算和本文方法的準確性，通過算例進行對比驗證。算例為計算一直立圓柱（見圖4）在方形航道中的水動力系數，航道寬100 m、水深50 m、圓柱半徑為15 m、圓柱吃水為40 m。將水動力系數計算結果與BV 船級社的HydroStar 軟件的結果進行比較，此項比較計算中不增加自由面數值阻尼（Hydrostar 的池壁格林函數計算無法考慮自由面數值阻尼）。

圖4 算例圓柱體示意圖（含自由面數值阻尼面元）

附加質量和阻尼系數的計算結果如圖5-圖14所示（圖中標記“hstar”為HydroStar 的計算結果；“my”為本文的計算結果，水動力系數進行了無因次化處理。ρ為流體密度，V為圓柱排水體積，d為圓柱直徑）。由于在勢流理論下圓柱艏搖的附加質量和阻尼系數為0，所以結果中不包含艏搖模態。

圖5 縱蕩附加質量

圖6 橫蕩附加質量

圖7 垂蕩附加質量

圖8 橫搖附加質量

圖9 縱搖附加質量

圖10 縱蕩阻尼系數

圖11 橫蕩阻尼系數

圖12 垂蕩阻尼系數

圖13 橫搖阻尼系數

圖14 縱搖阻尼系數

從結果可以看出，本文的計算結果和HydroStar的計算結果非常吻合，驗證了本文計算方法的可靠性和計算精度，特別是池壁格林函數計算的準確性。結果中的橫蕩和橫搖曲線中出現了明顯的數值共振的頻率點，此共振現象是由于基于勢流理論的池壁效應造成的。

本文通過引入自由面數值阻尼來消除數值共振點，為驗證自由面數值阻尼的有效性，本文根據算例來進行驗證。在圖4 所示的計算模型下（數值阻尼層的尺寸為96 m（橫向）×60 m（縱向），數值阻尼系數ε= 0.5），考慮自由面數值阻尼之后的水動力系數計算對比結果如圖15～圖16 所示。圖中標記為TGF 的曲線為沒有考慮自由面數值阻尼的池壁效應計算結果；標記為TGF_damp_rec 的曲線為考慮自由面數值阻尼的池壁效應計算結果；標記為FG_Opensea 的曲線為在普通開敞水域自由面格林函數下的水動力系數計算結果。文中沒有列出阻尼系數的比較結果，計算表明，在引入人工數值阻尼之后，阻尼系數對數值阻尼系數ε較為敏感，后續將繼續研究如何合理地確定數值阻尼系數的大小。

圖15 橫蕩附加質量

圖16 橫搖附加質量

從結果中可以看出，施加了自由面數值阻尼之后，附加質量系數隨著頻率變化是一條光滑的曲線，嚴重的數值共振現象不再出現，奇異性問題得到了較好的解決。從而使計算結果更為合理，與實際情況也更為吻合。

6 結 語

本文利用自由面池壁格林函數來考慮池壁效應，對其進行了基于切比雪夫多項式逼近的高效和高精度計算，并通過自由面數值阻尼來解決計算中出現的數值共振問題，目前公開商業軟件尚不具備此項功能。本文方法為浮體在限制水域的運動等響應計算提供了一個更為合理的方法，具備較大的工程應用價值。下一步的工作將著重解決如何確定數值阻尼的大小，使計算更便捷、結果更合理。