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

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

（1.中國(guó)船舶科學(xué)研究中心，江蘇無(wú)錫 214082；2.上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240）

0 引 言

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

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

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

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

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

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

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

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

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

圖1 坐標(biāo)系Fig.1 Coordinate system

1.1 輻射力時(shí)域直接積分法

在勢(shì)流理論里總的速度勢(shì)ΦT可寫成如下表達(dá)式：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)三維時(shí)域自由面格林函數(shù)記憶項(xiàng)表達(dá)式作如下無(wú)因次處理（King，1987）[7]：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.4 間接法求解源密度速度勢(shì)

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

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

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

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

2 數(shù)值計(jì)算方法驗(yàn)證

2.1 三維時(shí)域自由面格林函數(shù)記憶項(xiàng)及其導(dǎo)數(shù)計(jì)算及驗(yàn)證

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

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

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

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

2.2 Wigley I船型輻射力計(jì)算及驗(yàn)證

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

圖4 Wigley I船在Fn=0.3時(shí)垂蕩附加質(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時(shí)縱搖附加質(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

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

2.3 Wigley I船型強(qiáng)制運(yùn)動(dòng)產(chǎn)生的時(shí)域輻射力

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

圖6 Wigley I船不同航速時(shí)強(qiáng)制垂蕩/縱搖運(yùn)動(dòng)在垂蕩/縱搖方向產(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ù)美國(guó)密歇根大學(xué)Beck 團(tuán)隊(duì)的三維時(shí)域自由面格林函數(shù)理論公式，推導(dǎo)出可直接用于程序開發(fā)的三維時(shí)域自由面格林函數(shù)記憶項(xiàng)及其導(dǎo)數(shù)的數(shù)學(xué)展開表達(dá)式，通過(guò)Wigley I船型驗(yàn)證了本文方法的可靠性，得出如下結(jié)論：

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

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

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

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