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

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

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

0 引 言

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

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

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

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

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

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

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

1 數(shù)學(xué)模型

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

圖1 坐標系Fig.1 Coordinate system

1.1 輻射力時域直接積分法

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

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

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

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

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

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

式（4）離散成可編程的數(shù)學(xué)表達式：

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

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

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

時域輻射力求解的關(guān)鍵是求出輻射速度勢?i，下述章節(jié)將論述輻射速度勢的求解方法。

1.2 三維時域自由面格林函數(shù)法

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

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

1.3 三維時域自由面格林函數(shù)及其導(dǎo)數(shù)計算方法

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

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

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

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

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

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

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

文獻[7]在第一項中給出了前3個表達式，從軟件開發(fā)角度，推導(dǎo)給出該區(qū)域第一項可編程的數(shù)學(xué)表達式為

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

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

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

根據(jù)公式（16），從軟件開發(fā)角度，推導(dǎo)給出該區(qū)域?qū)?yīng)的三維時域格林函數(shù)記憶項第一項導(dǎo)數(shù)的數(shù)學(xué)表達式為

根據(jù)公式（17），從軟件開發(fā)角度，推導(dǎo)給出該區(qū)域格林函數(shù)記憶項公式第二項第一部分導(dǎo)數(shù)的可編程數(shù)學(xué)表達式如下，其他部分導(dǎo)數(shù)推導(dǎo)原理類同。

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

根據(jù)公式（21），推導(dǎo)給出該區(qū)域格林函數(shù)記憶項的第一項導(dǎo)數(shù)的可編程數(shù)學(xué)表達式：

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

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

根據(jù)公式（27）和公式（28），推導(dǎo)給出該區(qū)域格林函數(shù)記憶項的導(dǎo)數(shù)的可編程數(shù)學(xué)表達式：

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

當μ≈1時，格林函數(shù)記憶項及其導(dǎo)數(shù)表達式如下：

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

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

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

公式（35）～（36）離散寫成可編程的數(shù)學(xué)表達式如下：

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

2 數(shù)值計算方法驗證

2.1 三維時域自由面格林函數(shù)記憶項及其導(dǎo)數(shù)計算及驗證

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

圖2 無量綱化格林函數(shù)及其導(dǎo)數(shù)Fig.2 Nondimensional Green function and its derivative

圖3 無量綱化格林函數(shù)Fig.3 Nondimensional Green function

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

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

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

圖4 Wigley I船在Fn=0.3時垂蕩附加質(zhì)量、阻尼系數(shù)以及垂蕩縱搖耦合附加質(zhì)量、阻尼系數(shù)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時縱搖附加質(zhì)量、阻尼系數(shù)以及縱搖垂蕩耦合附加質(zhì)量、阻尼系數(shù)Fig.5 Pitch-pitch added mass and damping coefficient,pitch-heave added mass and damping coefficient of Wigley I at Fn=0.3

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

2.3 Wigley I船型強制運動產(chǎn)生的時域輻射力

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

圖6 Wigley I船不同航速時強制垂蕩/縱搖運動在垂蕩/縱搖方向產(chǎn)生的力/力矩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 結(jié) 論

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

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

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

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

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