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

（中國(guó)船舶及海洋工程設(shè)計(jì)研究院 上海200011）

引 言

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

線性勢(shì)流理論水動(dòng)力計(jì)算無(wú)法計(jì)入粘性和非線性效應(yīng)，導(dǎo)致在計(jì)算部分問(wèn)題時(shí)會(huì)出現(xiàn)數(shù)值共振問(wèn)題從而引起計(jì)算結(jié)果失真。為處理這一問(wèn)題，可在勢(shì)流方法中引入人工數(shù)值耗散項(xiàng)［7］，即數(shù)值阻尼來(lái)解決。本文將自由面數(shù)值阻尼引入限制水域的浮體運(yùn)動(dòng)計(jì)算，并通過(guò)自由面積分的解析解較好地解決了這一問(wèn)題。本文的所有計(jì)算均通過(guò)Fortran編程計(jì)算求解。

1 控制方程和定解條件

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

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

式（5）等號(hào)右側(cè)為開(kāi)敞有限水深自由面格林函數(shù)的無(wú)窮級(jí)數(shù)（m為整數(shù)），其表達(dá)式為：

Gm的級(jí)數(shù)表達(dá)式為：

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

K0為0 階修正的貝塞爾函數(shù)；

H0為0 階漢克爾函數(shù)。

2 自由面池壁格林函數(shù)計(jì)算

由式（5）可知，池壁格林函數(shù)可通過(guò)有限水深自由面格林函數(shù)展開(kāi)的無(wú)窮級(jí)數(shù)得到。但是由于式（5）級(jí)數(shù)收斂非常慢，直接進(jìn)行求和計(jì)算的成本非常高，實(shí)際工程應(yīng)用中幾乎是不可行的，所以需要找到高效的數(shù)值方法進(jìn)行計(jì)算。

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

式中：包括近場(chǎng)部分GN，中場(chǎng)部分GM和遠(yuǎn)場(chǎng)部分GF。近場(chǎng)采用級(jí)數(shù)直接求和的方法：

中場(chǎng)部分利用式（7），并注意到修正的貝塞爾函數(shù)漸近表達(dá)式為

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

式中：

er為剩余偏差項(xiàng)；

B=k0b，表示無(wú)因次化的航道寬度，系數(shù)

3 數(shù)值方法與切比雪夫多項(xiàng)式逼近

Newman［8］給出了近場(chǎng)部分GN的詳細(xì)計(jì)算方法，其中為了提高計(jì)算效率，對(duì)近場(chǎng)部分進(jìn)行了切比雪夫多項(xiàng)式逼近，中場(chǎng)部分GM計(jì)算較簡(jiǎn)單，可通過(guò)級(jí)數(shù)求和以及有理分式逼近或直接采用現(xiàn)有的函數(shù)庫(kù)來(lái)求解，以上計(jì)算方法無(wú)需贅述。本節(jié)主要介紹遠(yuǎn)場(chǎng)GF的高效計(jì)算方法以及對(duì)其進(jìn)行切比雪夫多項(xiàng)式逼近方法。

其中，剩余項(xiàng)的漸近表達(dá)式LerExtract（z，s，a）的奇異性與Φ（z，s，a）一致，所以原函數(shù)減去奇異項(xiàng)之后的剩余項(xiàng)可通過(guò)多變量切比雪夫多項(xiàng)式來(lái)逼近：

式中：Cjk為逼近多項(xiàng)式的系數(shù)；Tj和Tk分別是關(guān)于z和a的切比雪夫多項(xiàng)式；截?cái)囗?xiàng)數(shù)J和K根據(jù)不同的區(qū)域來(lái)選取，通?？赏ㄟ^(guò)試算確定。

為更好說(shuō)明式（31）， 現(xiàn)取s=1/2，a=2/3，b∈（0，π/2）分別對(duì)Φ（z，s，a）、LerExtract（z，s，a）、Φ（z，s，a）-LerExtract（z，s，a）作圖，對(duì)應(yīng)下頁(yè)圖1-圖3?？梢钥闯觯害担▃，s，a）和LerExtract（z，s，a）在b= 0 處存在奇異性，剩余項(xiàng)不再存在奇異性，且為光滑函數(shù)，可用切比雪夫多項(xiàng)式高效地逼近。

圖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 自由面人工數(shù)值阻尼施加

由于勢(shì)流理論無(wú)法計(jì)入粘性效應(yīng)等因素，采用池壁格林函數(shù)求解浮體運(yùn)動(dòng)時(shí)在某些特定的頻率會(huì)產(chǎn)生數(shù)值共振，造成計(jì)算結(jié)果與實(shí)際情況有較大的差異。為了較好地處理數(shù)值計(jì)算所帶來(lái)的奇異性問(wèn)題，可通過(guò)引入人工數(shù)值耗散項(xiàng)來(lái)解決。通常的做法是在自由面條件中引入數(shù)值阻尼ε［7］。

當(dāng)場(chǎng)點(diǎn)在自由面上時(shí)，積分方程為：

式中：σ為源強(qiáng)；G和Gn分別用池壁格林函數(shù)及其法向?qū)?shù)代入。由于池壁格林函數(shù)在自由面上存在奇異性，為保證計(jì)算精度，需要小心處理。可通過(guò)去奇點(diǎn)法進(jìn)行數(shù)值積分計(jì)算，即積分核中先減去一奇異項(xiàng)，此奇異項(xiàng)的積分可給出解析表達(dá)式。阻尼系數(shù)ε取一小值，通常可通過(guò)試算和試驗(yàn)方法確定。

5 算例與比較分析

為驗(yàn)證池壁格林函數(shù)計(jì)算和本文方法的準(zhǔn)確性，通過(guò)算例進(jìn)行對(duì)比驗(yàn)證。算例為計(jì)算一直立圓柱（見(jiàn)圖4）在方形航道中的水動(dòng)力系數(shù)，航道寬100 m、水深50 m、圓柱半徑為15 m、圓柱吃水為40 m。將水動(dòng)力系數(shù)計(jì)算結(jié)果與BV 船級(jí)社的HydroStar 軟件的結(jié)果進(jìn)行比較，此項(xiàng)比較計(jì)算中不增加自由面數(shù)值阻尼（Hydrostar 的池壁格林函數(shù)計(jì)算無(wú)法考慮自由面數(shù)值阻尼）。

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

附加質(zhì)量和阻尼系數(shù)的計(jì)算結(jié)果如圖5-圖14所示（圖中標(biāo)記“hstar”為HydroStar 的計(jì)算結(jié)果；“my”為本文的計(jì)算結(jié)果，水動(dòng)力系數(shù)進(jìn)行了無(wú)因次化處理。ρ為流體密度，V為圓柱排水體積，d為圓柱直徑）。由于在勢(shì)流理論下圓柱艏搖的附加質(zhì)量和阻尼系數(shù)為0，所以結(jié)果中不包含艏搖模態(tài)。

圖5 縱蕩附加質(zhì)量

圖6 橫蕩附加質(zhì)量

圖7 垂蕩附加質(zhì)量

圖8 橫搖附加質(zhì)量

圖9 縱搖附加質(zhì)量

圖10 縱蕩阻尼系數(shù)

圖11 橫蕩阻尼系數(shù)

圖12 垂蕩阻尼系數(shù)

圖13 橫搖阻尼系數(shù)

圖14 縱搖阻尼系數(shù)

從結(jié)果可以看出，本文的計(jì)算結(jié)果和HydroStar的計(jì)算結(jié)果非常吻合，驗(yàn)證了本文計(jì)算方法的可靠性和計(jì)算精度，特別是池壁格林函數(shù)計(jì)算的準(zhǔn)確性。結(jié)果中的橫蕩和橫搖曲線中出現(xiàn)了明顯的數(shù)值共振的頻率點(diǎn)，此共振現(xiàn)象是由于基于勢(shì)流理論的池壁效應(yīng)造成的。

本文通過(guò)引入自由面數(shù)值阻尼來(lái)消除數(shù)值共振點(diǎn)，為驗(yàn)證自由面數(shù)值阻尼的有效性，本文根據(jù)算例來(lái)進(jìn)行驗(yàn)證。在圖4 所示的計(jì)算模型下（數(shù)值阻尼層的尺寸為96 m（橫向）×60 m（縱向），數(shù)值阻尼系數(shù)ε= 0.5），考慮自由面數(shù)值阻尼之后的水動(dòng)力系數(shù)計(jì)算對(duì)比結(jié)果如圖15～圖16 所示。圖中標(biāo)記為T(mén)GF 的曲線為沒(méi)有考慮自由面數(shù)值阻尼的池壁效應(yīng)計(jì)算結(jié)果；標(biāo)記為T(mén)GF_damp_rec 的曲線為考慮自由面數(shù)值阻尼的池壁效應(yīng)計(jì)算結(jié)果；標(biāo)記為FG_Opensea 的曲線為在普通開(kāi)敞水域自由面格林函數(shù)下的水動(dòng)力系數(shù)計(jì)算結(jié)果。文中沒(méi)有列出阻尼系數(shù)的比較結(jié)果，計(jì)算表明，在引入人工數(shù)值阻尼之后，阻尼系數(shù)對(duì)數(shù)值阻尼系數(shù)ε較為敏感，后續(xù)將繼續(xù)研究如何合理地確定數(shù)值阻尼系數(shù)的大小。

圖15 橫蕩附加質(zhì)量

圖16 橫搖附加質(zhì)量

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

6 結(jié) 語(yǔ)

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