船舶有航速時三維時域格林函數計算方法研究

魯 江，張 楠，張新曙，顧 民

（1.中國船舶科學研究中心，江蘇無錫 214082；2.上海交通大學海洋工程國家重點實驗室，上海 200240）

0 引 言

船舶CAE 軟件是國產船舶行業工業軟件體系的重要組成部分，對基于三維時域面元法的勢流求解器開發提出了需求。在上個世紀80年代和90年代，基于邊界元法的三維時域面元法得到了快速發展，但局限于計算機計算速度，主要停留在算法驗證階段，工程應用上仍以頻域為主。2000年以來，隨著計算機軟硬件的快速發展，三維時域面元法的工程應用得到快速發展。2020年以來，船舶工業界把基于勢流求解器和粘流求解器的國產船舶流體CAE軟件開發提上了議程。三維時域面元法是直接在時域內建立初邊值問題進行求解，在處理瞬態問題（抨擊、沖擊作用）、非線性失穩大幅運動、波浪中機動操縱等問題時具有頻域方法不可替代的優勢，但三維時域邊界積分方程，不僅包含面積分和線積分，而且還包含整個歷程的時間卷積計算，計算量巨大。由于三維時域面元法的工程應用優點以及近代計算機的快速發展，基于三維時域面元法的勢流求解器開發重新引起船舶界的重視。

Finkelstain（1957）[1]提出了線性自由面下的三維時域格林函數；Wehausen 和Laitone（1960）[5]給出了有航速深水三維時域格林函數積分形式；Cummins（1962）[2]和Ogilvie（1964）[3]首先討論了非定常運動問題的時域直接求解方法；Wehausen（1967）[4]詳細討論了零航速時非定常運動問題，給出了零航速三維時域問題的Haskind關系。目前，針對該格林函數主要有三種求解方法，介紹如下。

美國密歇根大學Beck 團隊的Liapis（1986）[6]、King（1987）[7]、Magee（1991）[8]研究了船舶有航速時的三維時域格林函數法，提出一種三維時域格林函數計算方法，根據μ和β值的范圍將時域格林函數分成（μ＜0.7,β＜10）、（μ＜0.7,β≥10）、（0.7≤μ＜0.98）、（0.98≤μ≤1,β＜10）、（0.98≤μ≤1,β≥10）5個區域，每個區域用不同的級數展開和漸進展開進行計算。

Newman（1985，1990）[9-10]提出滿足10E-6和10E-10計算精度的兩套三維時域格林函數計算方法，Newman 團隊的Bingham（1994）[11]采用三維時域面元法進一步驗證了該三維時域格林函數方法解決Newman-Kelvin 線性運動的有效性。Lin 和Yue（1990）[12]基于Newman（1985）[9]的三維時域格林函數計算方法，根據μ和β值的范圍將時域格林函數分成5個區域，每個區域用不同的級數展開、漸進展開和剩余區間Chebyshev 正交逼近進行計算，依據該格林函數的計算方法，Lin 和Yue（1999）等[13]進一步提出一種滿足線性自由面條件的數值方法，將流場分為內場和外場，以內場Rankine源方法為核心，外場匹配時域格林函數來滿足輻射條件，即混合源法。以此為基礎，美國海軍開發了船舶大幅運動評估軟件LAMP，分為LAMP1線性版本、LAMP2弱非線性版本和LAMP4全非線性版本（Shin，2003）[14]。

黃德波（1992）[15]導出了三維時域格林函數及其導數的簡化計算公式，對于β＞8+1.515μ采用漸進展開，對于β＜8+1.515μ使用級數展開，提出造表插值方案和相應的節點劃分辦法，三維時域格林函數法開始在國內得到學者的廣泛應用；周正全（1993）等[16]提出了反向輻射勢表示繞射力的關系式，可認為是Wehausen（1967）[4]的無航速時域Haskind關系的發展，稱之為有航速時的Haskind關系，其三維時域格林函數的計算采用黃德波（1992）[15]方法；王大云（1996）[17]提出了一種三維水彈性時域分析方法，其三維時域格林函數的計算采用黃德波（1992）[15]方法；卜淑霞等（2018）[18]和儲紀龍（2019）等[19]采用三維時域混合源法研究了參數橫搖現象，其三維時域格林函數的計算采用Newman（1985）[9]和Lin和Yue（1990）[12]的計算方法；陳紀康、段文洋等（2021）[20]采用三維時域混合源法構建邊界積分方程預報船舶運動，其外域三維時域格林函數及其導數采用黃德波（1992）[15]的計算方法，并利用泰勒展開邊界元法求解該邊界積分方程，解決了有航速定解問題中涉及的非光滑邊界處切向誘導速度精度低的問題。

Clement（1998）[21]提出了三維時域格林函數及其導數的四階微分方程方法；Duan 和Dai（2001）[22]、申亮和朱仁傳等（2007）[23]分別對三維時域格林函數及其導數的四階微分方程方法做了推導和部分改進；Li、Ren 等（2015）[24]提出了三維時域格林函數及其導數的四階微分方程精細積分方法，對Clement（1998）[21]的四階微分方程法做出了顯著改進；楊鵬（2016）[25]給出了一種三維時域水彈性運動方程的數值求解方法，其三維時域格林函數的計算采用Clement（1998）[21]的四階微分方程法；周文俊等（2021）[26]提出了基于垂向積分形式時域格林函數的多域高階面元法，其三維時域格林函數的計算采用Clement（1998）[21]的四階微分方程法。

Newman（1990）[10]指出Beck團隊對三維時域格林函數計算方法做出了實質性改進，但國內公開應用的很少，本文依據美國密歇根大學Beck 團隊的三維時域自由面格林函數理論方法，推導出可直接用于程序開發的三維時域自由面格林函數記憶項完整數學展開表達式，推導給出可直接用于程序開發的三維時域自由面格林函數記憶項導數的完整數學展開表達式，給出可編程應用的參數取值和求解說明，給出了本文計算結果和公開的試驗結果、計算結果以及作者使用的加強切片法計算結果的對比驗證，使得三維時域勢流求解器的工程開發簡單化，可促進波浪中操縱性、快速性和穩性等學科的發展。

1 數學模型

空間固定坐標系O0-X0Y0Z0原點位于靜水面，X0軸為波浪傳播方向，Y0軸指向左為正，Z0軸垂直水面向上為正。參考坐標系o-xyz隨著船體以恒定速度U0沿著x正方向前進，初始時刻和空間坐標系重合，不隨船體轉動而轉動。隨船坐標系o-x'y'z'固定于船體上，原點和參考坐標相同，船體正浮平衡時，與參考坐標系重合，隨船體的轉動而轉動。坐標系如圖1 所示。設入射波浪向角為β，船舶航向角為χ，則β=-χ。

圖1 坐標系Fig.1 Coordinate system

1.1 輻射力時域直接積分法

在勢流理論里總的速度勢ΦT可寫成如下表達式：

式中，p(x,y,z,t)為t時刻坐標(x,y,z)處的壓力，SB( )t為t時刻浮體濕表面，ρ為液體密度，g為重力加速度。

根據壓力和力/力矩的關系，忽略二階速度勢小量，由船體i模態運動引起的作用在j自由度方向上輻射力/力矩可以表示為

式中，j=1,2,…,6分別對應縱蕩、橫蕩、垂蕩、橫搖、縱搖和首搖6個方向。

為了消除船體物面速度勢空間導數的計算，許多學者采用Tuck[27]理論假定。公式（3）的第二項是關于速度勢的前進速度影響項，適用于分部積分法，該積分項在船首到船尾整船積分為零，公式（3）可以寫成下面表達式：

式中，nj為浮體濕表面上在j方向矢量，n1是小量，n2、n3在x方向變化很慢，

式（4）離散成可編程的數學表達式：

式中，M是浮體濕表面上面元的個數，Ap是第p個面元的面積，[nj]p是第p個面元上的nj，[mj]p是第p個面元上的mj。

在勢流理論里，輻射速度勢?i(x,y,z,t)物面邊界條件滿足如下表達式，本文計算采用船體靜水中濕表面。

式中，xi(t),x?i(t)分別是t時刻浮體i方向模態運動的位移和速度。

時域輻射力求解的關鍵是求出輻射速度勢?i，下述章節將論述輻射速度勢的求解方法。

1.2 三維時域自由面格林函數法

在勢流理論里，時域求解非定常速度勢的方法主要有時域自由面格林函數法和時域Rankine 源法。自由面格林函數法滿足線性自由面條件以及除物面條件外的其他所有邊界條件，僅需在浮體濕表面上分布奇點，數值離散化的方程只含有浮體濕表面及水線上的未知量。

Wehausen 和Laitone（1960）[5]給出了深水有航速三維時域格林函數積分形式，Liapis（1986）[6]和King（1987）[7]采用脈沖函數法。三維時域格林函數積分形式可寫成如下表達式：

1.3 三維時域自由面格林函數及其導數計算方法

三維時域自由面格林函數的瞬時項采用Hess&Smith 方法，參見文獻[28]，不再論述。下面主要描述本文中三維時域自由面格林函數的記憶項具體采用的計算方法及推導給出的編程數學表達式。

對三維時域自由面格林函數記憶項表達式作如下無因次處理（King，1987）[7]：

三維時域自由面格林函數記憶項對空間導數的無因次表達式如下，和（King，1987）[7]一致：

參考文獻[7]，本文時域格林函數記憶項求解分為五個區域，第一、二區域的分界線β為8.8，第三、四區域的分界線為0.97，并推導了可直接用于程序開發的三維時域格林函數記憶項及其導數的數學展開表達式。

第一區域，0≤μ＜0.7，0≤β＜8.8，時域格林函數記憶項采用Legendre 多項式的級數展開，文獻[6]給出了公式的前4項，從軟件開發角度，推導給出可編程的數學表達式如下：

n的取值以計算結果收斂到計算精度為準，文獻[7]沒有給出具體值，n大于80時已滿足計算精度10E-10的要求。

第二區域，0≤μ＜0.7，β≥8.8，采用漸進展開，格林函數公式分成兩部分來論述。

文獻[7]在第一項中給出了前3個表達式，從軟件開發角度，推導給出該區域第一項可編程的數學表達式為

n的取值以計算結果收斂為準，文獻[7]沒有給出具體值，n大于100 時已滿足計算精度10E-6 的要求。

該區域格林函數第二項表達式（King，1987）[7]如下：

式中，θ=sin-1(μ)。

根據公式（16），從軟件開發角度，推導給出該區域對應的三維時域格林函數記憶項第一項導數的數學表達式為

根據公式（17），從軟件開發角度，推導給出該區域格林函數記憶項公式第二項第一部分導數的可編程數學表達式如下，其他部分導數推導原理類同。

文獻[7]定義(γh)min為小值，但沒有給出具體值。本文(γh)min的值設為0.35，Δk=h=0.05，n的值為kmax/(2h)取整數。

根據公式（21），推導給出該區域格林函數記憶項的第一項導數的可編程數學表達式：

第四區域，0.97≤μ＜1，β＜10，采用貝塞爾函數[7]展開。

公式（27）中n最大值為5，公式（28）中對應的m最大值為22。

根據公式（27）和公式（28），推導給出該區域格林函數記憶項的導數的可編程數學表達式：

第五區域，0.97≤μ＜1，β≥10，采用誤差函數的漸進展開（King，1987）[7]，推導給出可編程的數學表達式：n的取值以計算結果收斂為準，文獻[7]沒有給出具體值，n大于30時已滿足計算精度10E-6的要求。

當μ≈1時，格林函數記憶項及其導數表達式如下：

1.4 間接法求解源密度速度勢

同時在物面上布置源和偶極子的方法稱為直接法，僅布置源的方法稱為間接法。間接法可以通過邊界積分方程得到船體濕表面的切向速度，更適合求解考慮瞬時濕表面的非線性大幅運動問題。文獻[8]的公式（2.11）、（2.12）基于勢流理論給出了如下求解分布源密度σ(P,t)的邊界積分方程和速度勢?(P,t)表達式：

式中：SB(t)為t時刻浮體濕表面；SB(τ)為τ時刻浮體濕表面；Γ(τ)為τ時刻水面與船體濕表面的交線；n?為浮體濕表面面元的單位法向矢量，方向為指向浮體濕表面、離開流體域；N?是自由面上單位法向量，方向為指向Γ(τ)、離開流體域；VN是Γ(τ)在N?方向的法向速度，VN和Vn的關系為VN=Vn/N?·n?，VN可以近似取線元臨近面元的法向速度在該線元切向速度的投影。

公式（35）～（36）離散寫成可編程的數學表達式如下：

式中：M為浮體濕表面上面元個數，q為浮體濕表面上第q個面元，p為浮體濕表面上第p個面元；L為水線上的線元個數，l為水線上第l個線元；tN表示當前的時間步N對應的時間，tn表示開始計及記憶效應的第n步對應的時間，Δt為時間步長。

2 數值計算方法驗證

2.1 三維時域自由面格林函數記憶項及其導數計算及驗證

圖2 給出了無量綱化格林函數記憶項及其導數值。圖3 中“Present”表示本文計算結果，“Yang，2016”表示文獻Yang（2016）[25]的結果，后續圖中表示方法類似。本文計算結果和Yang（2016）[25]采用Clement（1998）的微分方程法計算的結果吻合。

圖2 無量綱化格林函數及其導數Fig.2 Nondimensional Green function and its derivative

圖3 無量綱化格林函數Fig.3 Nondimensional Green function

從圖中可以看出，在μ接近0 時（場點P和源點Q靠近水線），隨時間β的增大，格林函數及其導數的振蕩特性非常明顯，能夠影響時域波浪力求解的穩定性，是三維時域格林函數計算方法需要不斷改進的原因之一。隨著μ變大，格林函數及其導數值隨時間β的增大而減小，最后趨于零，μ越大，其趨于零的速度越快。

2.2 Wigley I船型輻射力計算及驗證

Wigley Ι 船型型值參見文獻[29]，本文船長設為9.81 m，船寬設為0.981 m，吃水設為0.613 125 m。ω為搖蕩頻率，L為船長，Δ 為排水體積。圖4 和圖5 中“Present”表示本文計算結果，“Exp_Journee_1992”表示文獻[29]的試驗結果，“Bingham_1994”表示文獻[11]的計算結果，“Estrip”表示本文采用加強切片法計算的結果。

圖4 Wigley I船在Fn=0.3時垂蕩附加質量、阻尼系數以及垂蕩縱搖耦合附加質量、阻尼系數Fig.4 Heave-heave added mass and damping coefficient,heave-pitch added mass and damping coefficient of Wigley I at Fn=0.3

圖5 Wigley I船在Fn=0.3時縱搖附加質量、阻尼系數以及縱搖垂蕩耦合附加質量、阻尼系數Fig.5 Pitch-pitch added mass and damping coefficient,pitch-heave added mass and damping coefficient of Wigley I at Fn=0.3

通過三維時域格林函數計算Wigley I 船型的垂蕩、縱搖時域輻射力，根據時域輻射力曲線的振幅及其和波浪的相對相位，整理對應附加質量和阻尼系數（Present），并和公開的試驗結果（Journée，1992）[29]、計算結果（Bingham，1994）[11]以及作者采用加強切片法（Kashiwagi et al，2010）[30]計算的結果（EStrip）進行對比。本文網格數設為240，時間步長設為0.05 s，垂蕩方向時域輻射力計算采用公式（3），縱搖方向時域輻射力計算采用公式（4）。圖4給出了強制垂蕩運動引起的垂向附加質量、阻尼系數以及在縱搖方向的附加質量和阻尼系數，圖5給出了強制縱搖運動引起的縱搖方向的附加質量、阻尼系數以及垂蕩方向附加質量和阻尼系數。本文計算結果和（Journée，1992）[29]試驗結果、（Bingham，1994）[11]計算結果相比，整體吻合較好，但偏差還是存在。和（Journée，1992）[29]試驗結果出現偏差的原因主要在于真實流體是黏性的，強制搖蕩試驗時還有自由面影響，而計算是基于理想的理論假定，且只考慮平均濕表面。和（Bingham，1994）[11]計算結果出現偏差的原因主要是計算資源的限制，網格、時間步長取值不一致，水線的計算處理不完全一致。作者同時采用加強切片法（EStrip）計算了水動力系數，和切片法計算結果相比，三維時域方法整體上優于切片理論的計算結果。后續采用三維時域混合源法進一步驗證三維時域理論的有效性。

2.3 Wigley I船型強制運動產生的時域輻射力

通過三維時域格林函數進一步計算Wigley I船型在不同航速時的時域輻射力，圖6給出了強制垂蕩/縱搖振幅為0.05 m/0.05 rad，頻率為3 rad/s，航速Fn=0.0、0.3、0.6 時在垂蕩/縱搖方向產生的時域力/力矩。盡管計算是基于靜水濕表面，但力的求解公式（4）和物面條件公式（6）都含有和航速有關的m項，同時時域格林函數中場點和源點的距離和航速有關。航速增大后，該方法仍可以計算出穩定的垂蕩和縱搖方向的時域輻射力/力矩，只是力/力矩振幅增大，相位變化微小。根據時域輻射力/力矩曲線的振幅及其和波浪的相位差，可以整理出對應的附加質量和阻尼系數。

圖6 Wigley I船不同航速時強制垂蕩/縱搖運動在垂蕩/縱搖方向產生的力/力矩Fig.6 Force/moment in the heave and pitch directions generated by the forced heave and pitch motions of Wigley I at different Fn

3 結 論

本文依據美國密歇根大學Beck 團隊的三維時域自由面格林函數理論公式，推導出可直接用于程序開發的三維時域自由面格林函數記憶項及其導數的數學展開表達式，通過Wigley I船型驗證了本文方法的可靠性，得出如下結論：

（1）本文采用的方法能夠準確得出三維時域自由面格林函數記憶項及其導數值，重現格林函數及其導數的振蕩特性。

（2）本文采用的三維時域自由面格林函數法能夠可靠地求解Wigley I船型的垂蕩、縱搖及其耦合方向的輻射力。

（3）本文推導給出了三維時域自由面格林函數記憶項及其導數的編程數學展開表達式，有助于從事此方面研究的人員進行編程應用，實用性強。

本文方法和代碼可用于船舶工業CAE 軟件勢流求解器，實船工程應用上需要在規避三維時域格林函數的振蕩特性以及外飄船型、尾部淺吃水的速度勢發散方面做深入研究。