姚朝幫，董文才

(海軍工程大學(xué)艦船工程系，湖北武漢430033)

開(kāi)爾文源格林函數(shù)近年來(lái)得到了較為廣泛的應(yīng)用［1-4］，它由Havelock在1928年導(dǎo)出，現(xiàn)在被稱為Kelvin(或Havelock)源，其物理意義是以恒定速度直航的單位點(diǎn)源在場(chǎng)點(diǎn)所產(chǎn)生的速度勢(shì)，表示為雙積分形式，被積函數(shù)具有奇異性、振蕩性，直接采用數(shù)值積分計(jì)算很難保證計(jì)算精度與效率，也難以應(yīng)用到船舶水動(dòng)力問(wèn)題的計(jì)算。為此，很多學(xué)者將其變換為易于計(jì)算的形式，常見(jiàn)的有Peters型，Michell型和Havelock型，Noblesse歸納了這3種形式的推導(dǎo)過(guò)程［5］。不同表達(dá)形式有所不同，本文選用的是Havelock型的改進(jìn)形，即將Kelvin源格林函數(shù)表達(dá)成Rankine源及其關(guān)于靜水面的映像源項(xiàng)、近場(chǎng)擾動(dòng)項(xiàng)N和遠(yuǎn)場(chǎng)波動(dòng)項(xiàng)W。Rankine源及其關(guān)于靜水面的映像源項(xiàng)可采用Hess-Smith提出的方法計(jì)算。因此，本文重點(diǎn)討論擾動(dòng)項(xiàng)N和遠(yuǎn)場(chǎng)波動(dòng)項(xiàng)W的計(jì)算方法對(duì)比分析。

近場(chǎng)擾動(dòng)項(xiàng)N的計(jì)算方法主要有4種:1)制表法［6-8］，通過(guò)直接數(shù)值積分，形成關(guān)于擾動(dòng)項(xiàng)N及其偏導(dǎo)數(shù)的大規(guī)模數(shù)據(jù)。該方法的優(yōu)點(diǎn)在于制表與船型本身無(wú)關(guān)，制表一旦形成，可以用于不同船型興波流場(chǎng)的計(jì)算，但為了提高計(jì)算的精度，通常需要規(guī)模較大的制表，通?？蛇_(dá)到expt;2)切比雪夫多項(xiàng)式擬合法［9］，采用切比雪夫多項(xiàng)式對(duì)N進(jìn)行擬合計(jì)算。采用該方法時(shí)，N的計(jì)算精度能達(dá)到5D-6D，但偏導(dǎo)數(shù)的精度會(huì)有所降低;3)有理函數(shù)法，將格林函數(shù)雙積分形式變換成只包含有理函數(shù)的單積分形式［10-11］，在此基礎(chǔ)上采用梯形法或辛普森方法計(jì)算積分。但實(shí)際計(jì)算時(shí)需要權(quán)衡積分精度與效率;4)使用多項(xiàng)式對(duì)N進(jìn)行擬合［12］，從而簡(jiǎn)化計(jì)算，但該方法的精度及有效性還有待進(jìn)一步驗(yàn)證。

遠(yuǎn)場(chǎng)波動(dòng)項(xiàng)W的計(jì)算方法主要有3種:1)Bessel函數(shù)法，通過(guò)引入Dawson型積分將W轉(zhuǎn)化為Bessel函數(shù)的表達(dá)形式［13-15］，利用Bessel函數(shù)的相關(guān)性質(zhì)進(jìn)行計(jì)算。采用該方法時(shí)需要分區(qū)進(jìn)行求解，在區(qū)域過(guò)渡的地方，計(jì)算結(jié)果有時(shí)會(huì)出現(xiàn)跳躍;2)有理函數(shù)法，與求解N項(xiàng)時(shí)的方法相同;3)變量代換法，對(duì)W的單積分形式進(jìn)行變形，利用變量代換得到易于積分的形式。

上述N與W的計(jì)算方法各有特點(diǎn)，也被廣泛采用，但公開(kāi)發(fā)表的文獻(xiàn)很少涉及它們計(jì)算精度、效率的對(duì)比及優(yōu)選，基于此，本文對(duì)上述常用的方法進(jìn)行總結(jié)分析，并評(píng)估它們的精度、效率，從而找到該格林函數(shù)的快速數(shù)值計(jì)算方法，為將該格林函數(shù)應(yīng)用于船舶水動(dòng)力的計(jì)算奠定基礎(chǔ)。

1 Kelvin源格林函數(shù)

設(shè)移動(dòng)源的前進(jìn)速度為U，ξ=(ξ，η，ζ)為源點(diǎn)，X=(X，Y，Z)為場(chǎng)點(diǎn)，L為特征長(zhǎng)度，g為重力加速度，采用L、g及U對(duì)相關(guān)參量進(jìn)行無(wú)因次化，x=X/L，F(xiàn)=，則Kelvin源格林函數(shù)的表達(dá)式為［8］

G*對(duì)場(chǎng)點(diǎn)的偏導(dǎo)數(shù)的表達(dá)式可參見(jiàn)文獻(xiàn)［8］，其和G*的計(jì)算相比并未引入計(jì)算難度。

2 近場(chǎng)擾動(dòng)項(xiàng)N(a)的計(jì)算方法

2.1 制表法

制表法是通過(guò)直接積分計(jì)算N(a)，進(jìn)而得到大規(guī)模的N隨(a1，a2，a3)的變化規(guī)律，并以此作為插值數(shù)據(jù)庫(kù)。

數(shù)值積分計(jì)算時(shí)，對(duì)式(2)中N(a)的表達(dá)式進(jìn)一步變形，得到在t∈［-1，1］時(shí)比積函數(shù)更為光順且積分端點(diǎn)處奇異性消失的N(a)及?aN(a)表達(dá)式［8］:

制表時(shí)需將a1、a2、a3、d投射到［0，1］區(qū)間內(nèi)，因此定義如下變換:

為了方便插值計(jì)算，采用式(9)對(duì)N(a)及?aN(a)進(jìn)行變換:

本文選取的制表規(guī)模為193×45×28。式(3)～(6)的積分計(jì)算采用自適應(yīng)辛普森法。

2.2 切比雪夫多項(xiàng)式擬合方法

Newman從開(kāi)爾文源格林函數(shù)雙積分形式出發(fā)，得到了Kelvin源格林函數(shù)近場(chǎng)項(xiàng)的切比雪夫多項(xiàng)式擬合形式，并計(jì)算得到了擬合系數(shù)［9］。該方法N的計(jì)算式為

式中:f(d)==arcsin(a1/d)，

φ=arctan(a2/a3);Ti，Tj，T2k是切比雪夫多項(xiàng)式，文獻(xiàn)［15］給出了具體的迭代求解方法。S由式 (11)計(jì)算得到:

Vm、Um可分別通過(guò)式(12)、(13)求解:

式中:初值U0=

式中，s=

采用該方法計(jì)算Kelvin源格林函數(shù)時(shí)，N、W與其他計(jì)算方法的對(duì)應(yīng)關(guān)系為

采用自適應(yīng)辛普森法計(jì)算式(14)中的積分，可同時(shí)得到N及W。

2.3 有理函數(shù)法

Wang［10］等在Bessho得到的開(kāi)爾文源格林函數(shù)表達(dá)的基礎(chǔ)上，進(jìn)一步整理分析，得到了開(kāi)爾文源格林函數(shù)的有理函數(shù)表達(dá)形式為

2.4 改進(jìn)的梯形法

圖1 改進(jìn)梯形法流程Fig.1 Solution strategy for improved trapezoidal rule

分析式(4)～(7)可以看出N(a)及其偏導(dǎo)數(shù)只是在積分端點(diǎn)附近變化比較劇烈，其他區(qū)間的曲線較為平滑，有利于直接積分。借鑒自適應(yīng)梯形法和自適應(yīng)辛普森法的思想，根據(jù)被積函數(shù)的特性，將被積函數(shù)的積分區(qū)間分為［-1，-1+Δh］，［-1+Δh，1-Δh］和［1-Δh，1］，Δh為積分端點(diǎn)處的微元段。在區(qū)間［-1，-1+Δh］及［1-Δh，1］采用局部加密的方法進(jìn)行離散，并利用梯形法求得積分值，文中將這種方法稱之為“改進(jìn)的梯形法”，其基本流程見(jiàn)圖1。

3 遠(yuǎn)場(chǎng)波動(dòng)項(xiàng)W(a)的計(jì)算方法

3.1 Bessel函數(shù)法

Bessho通過(guò)引入Dawson型積分，將開(kāi)爾文源格林函數(shù)的波動(dòng)項(xiàng)轉(zhuǎn)換成了Bessel函數(shù)的級(jí)數(shù)形式［13］為

式中:Jn、Kn、Yn、Ⅰn為Bessel函數(shù)，J'2n、Y'2n為關(guān)于a1的偏導(dǎo)數(shù)。

對(duì)式(16)關(guān)于a1、a2、a3求偏導(dǎo)便可得到對(duì)應(yīng)的導(dǎo)數(shù)項(xiàng);求和項(xiàng)數(shù)由指定的計(jì)算精度確定。

3.2 變量代換法

由式(2)可知W(a)及其偏導(dǎo)數(shù)的表達(dá)式為

W(a)是振蕩的，振蕩的頻率與a1、a2有關(guān)，振蕩的幅值與a3有關(guān)，對(duì)于這類積分，若按常規(guī)的方法進(jìn)行積分，積分效率較低且很難得到正確的結(jié)果［16］，本文提出變量代換的方法來(lái)處理，以期能同時(shí)兼顧積分效率與精度。

設(shè)W(a)的實(shí)際積分區(qū)間為［tg，th］，令 τ=，則W(a)的被積函數(shù)在［tg，th］的積分可化為

式(19)是一個(gè)復(fù)函數(shù)和指數(shù)函數(shù)的積分，設(shè)該復(fù)函數(shù)為f，則f=(l-im)/(l2+m2)。

設(shè)［tg，th］離散后的序列為t1、t2、t3、…、tn，對(duì)應(yīng)的A序列為A(t1)、A(t2)、A(t3)、…、A(tn)，可將f在相鄰兩點(diǎn)A(ti)和A(ti+1)(1≤i≤n-1)間的函數(shù)用關(guān)于A(ti)的一次函數(shù)近似，那么被積函數(shù)在相鄰兩點(diǎn)間的積分Ⅰi可直接積出，整個(gè)積分區(qū)間的積分值為各Ⅰi的和:

經(jīng)過(guò)上述變量代換后，將振蕩函數(shù)轉(zhuǎn)換成了易于積分的表達(dá)形式，積分精度取決于f本身的函數(shù)特性及其所對(duì)應(yīng)的離散方法。為此，開(kāi)展其函數(shù)的性狀分析，并重點(diǎn)分析m=0的鄰域內(nèi)的變化情況。圖2給出了f在積分區(qū)間內(nèi)的2種典型函數(shù)曲線(圖中框點(diǎn)為m=0對(duì)應(yīng)的點(diǎn))。從圖中可以看出，在m=0對(duì)應(yīng)的t附近，f函數(shù)值變化劇烈，出現(xiàn)了一些看似奇異的現(xiàn)象，但函數(shù)本身是連續(xù)的，這只是函數(shù)f局部區(qū)域變化劇烈的“偽奇異”，稱之為“偽奇異現(xiàn)象”，與之對(duì)應(yīng)的點(diǎn)稱之為“偽奇異點(diǎn)”。

圖2 f的典型特性曲線Fig.2 Typical characteristics of f

為了保證積分的精度和效率，在對(duì)積分區(qū)間進(jìn)行離散時(shí)必須將這種“偽奇異性”考慮進(jìn)來(lái)。區(qū)間離散時(shí)需要在這些“偽奇異點(diǎn)”處局部加密處理。進(jìn)一步分析可知“偽奇異點(diǎn)”出現(xiàn)的位置可由式(21)所示的解析形式得到:

需要指出的是，式(18)及(19)中N(a)的積分區(qū)間是(-∞，+∞)，不利于數(shù)值積分的實(shí)施，實(shí)際數(shù)值計(jì)算過(guò)程中需要對(duì)積分區(qū)間進(jìn)行適當(dāng)截?cái)?，具體截?cái)嗟姆椒蓞⒁?jiàn)文獻(xiàn)［17］。

4 數(shù)值計(jì)算結(jié)果及分析

為了比較各種計(jì)算方法的精度，本文選用如圖3所示的算例，算例中點(diǎn)源坐標(biāo)為(0.0，0.0，-1.0)，并以速度2 m/s沿ox軸正方向運(yùn)動(dòng)，場(chǎng)點(diǎn)位于水平面上與ox軸成角度ф=arctan(0.5)的直線上。

評(píng)估各種計(jì)算方法的效率時(shí)仍選用上述算例，此時(shí)改變?cè)袋c(diǎn)的垂向坐標(biāo)，使其逐漸接近自由表面(Z=0)，計(jì)算時(shí)源點(diǎn)的垂向坐標(biāo)分別為Z=-1，-0.1，-0.01，-0.001，統(tǒng)計(jì)分析計(jì)算一次格林函數(shù)及其偏導(dǎo)數(shù)需要的平均時(shí)間。

圖3 移動(dòng)源點(diǎn)與場(chǎng)點(diǎn)的位置Fig.3 The location of translating source and field points

4.1 N(a)計(jì)算結(jié)果分析及效率評(píng)估

圖4中給出了源點(diǎn)位于(0，0，-1)時(shí)不同方法計(jì)算的N(a)及偏導(dǎo)數(shù)的對(duì)比曲線。表1中給出了計(jì)算一次N(a)及其偏導(dǎo)數(shù)的總耗時(shí)，其中M1、M2、M3、M4分別代表切比雪夫多項(xiàng)式擬合方法、制表法、有理函數(shù)法及改進(jìn)的梯形法。文中耗時(shí)統(tǒng)計(jì)均在一CPU為Intel Core i3(頻率為2.27 GHz)的PC機(jī)上進(jìn)行。

圖4 N(a)及其偏導(dǎo)數(shù)的函數(shù)圖像Fig.4 Images of N(a)and its derivatives

表1 N(a)計(jì)算耗時(shí)對(duì)比Table 1 Computation time comparison of N(a)

從圖4中可以看出，不同方法計(jì)算的N(a)及其偏導(dǎo)數(shù)值基本相同，這說(shuō)明幾種方法的計(jì)算精度一致，從而互相驗(yàn)證了方法的有效性。從表1耗時(shí)對(duì)比可以看出，采用Newman導(dǎo)出的切比雪夫多項(xiàng)式擬合法及制表法效率較高，改進(jìn)的梯形法次之，采用有理函數(shù)法計(jì)算耗時(shí)較長(zhǎng)，主要原因在于3個(gè)偏導(dǎo)數(shù)計(jì)算公式繁瑣，求解較費(fèi)時(shí)。

4.2 W(a)計(jì)算結(jié)果分析及效率評(píng)估

仍采用上述算例，采用Bessel函數(shù)法、變量代換法及有理函數(shù)法分別計(jì)算W(a)，統(tǒng)計(jì)各種計(jì)算方法的耗時(shí)與計(jì)算精度。不同方法計(jì)算得到W(a)及其偏導(dǎo)數(shù)的對(duì)比曲線如圖5所示。表2中給出了各種計(jì)算方法計(jì)算一次W(a)及其偏導(dǎo)數(shù)的平均耗時(shí)，表中M4、M5、M6分別代表Bessel函數(shù)法，變量代換法及有理函數(shù)法。由圖5及表2可見(jiàn):3種方法的計(jì)算精度基本相同，變量代換方法的效率明顯優(yōu)于其他2種方法。

表2 W(a)計(jì)算耗時(shí)對(duì)比Table 2 Computation time comparison of W(a)

圖5 W及其偏導(dǎo)數(shù)的函數(shù)圖像Fig.5 Images of W and its derivatives

4.3 G*的計(jì)算結(jié)果及驗(yàn)證

圖6 G*及其偏導(dǎo)數(shù)的函數(shù)圖像Fig.6 Images of G*and its derivatives

從4.1 節(jié)、4.2 節(jié)的對(duì)比可以看出:N(a)的計(jì)算采用切比雪夫多項(xiàng)式擬合法和制表法效率較高;W(a)的計(jì)算采用“變量代換法”效率較高，因此本節(jié)將綜合采用切比雪夫多項(xiàng)式擬合法與“變量代換法”計(jì)算G*。仍采用上述算例，計(jì)算得到的G*及其偏導(dǎo)數(shù)與文獻(xiàn)［18］的對(duì)比如圖6所示。

從圖6中可以看出本節(jié)選用的近場(chǎng)項(xiàng)及波動(dòng)項(xiàng)的計(jì)算方法計(jì)算G*及偏導(dǎo)數(shù)的精度與文獻(xiàn)［18］基本相同，計(jì)算一次的平均耗時(shí)約為6～9 ms，能滿足船舶興波流場(chǎng)計(jì)算的工程需要。

5 結(jié)論

本文在分析開(kāi)爾文源格林函數(shù)近場(chǎng)擾動(dòng)項(xiàng)和遠(yuǎn)場(chǎng)波動(dòng)項(xiàng)計(jì)算方法的基礎(chǔ)上，對(duì)比分析了各種計(jì)算方法的計(jì)算精度與效率，通過(guò)研究得到以下結(jié)論:

1)在研究的場(chǎng)點(diǎn)、源點(diǎn)位置范圍內(nèi)，4種方法計(jì)算近場(chǎng)項(xiàng)N(a)及其偏導(dǎo)數(shù)的精度基本相同，但效率各異，其中切比雪夫多項(xiàng)式擬合法和制表插值法計(jì)算效率較高，在實(shí)際計(jì)算中可優(yōu)先選擇;“改進(jìn)的梯形法”在保證精度的同時(shí)效率稍低;有理函數(shù)法因偏導(dǎo)數(shù)表達(dá)式較復(fù)雜、計(jì)算效率降低，實(shí)際使用時(shí)不建議采用。

2)對(duì)比分析了波動(dòng)項(xiàng)N(a)的3種計(jì)算方法，它們的計(jì)算精度基本相當(dāng)，“變量代換法”效率較高，采用Bessel函數(shù)級(jí)數(shù)展開(kāi)的計(jì)算方法次之，Bessho型表達(dá)計(jì)算耗時(shí)較長(zhǎng)。

3)綜合采用切比雪夫多項(xiàng)式擬合法計(jì)算N(a)與變量代換法計(jì)算W(a)，能快速得到開(kāi)爾文源格林函數(shù)的計(jì)算得結(jié)果(一次計(jì)算的平均耗時(shí)約為6～9 ms)，可作為船舶水動(dòng)力求解的基本工具。

［1］李子如，陳克強(qiáng).用面元法預(yù)報(bào)船舶的興波阻力及波型［J］.武漢理工大學(xué)學(xué)報(bào)，2006，30(6):1094-1097.LI Ziru，CHEN Keqiang.Panel method for the prediction of ship wave resistance and wave pattern［J］.Journal of Wuhan University of Technology，2006，30(6):1094-1097.

［2］田超，吳有生.水線源強(qiáng)消棄法在船舶定常興波計(jì)算中的應(yīng)用［J］.船舶力學(xué)，2008，12(1):37-45.TIAN Chao，WU Yousheng.Application of the waterline source strength elimination method in the computation of the steady ship waves［J］.Journal of Ship Mechanics，2008，12(1):37-45.

［3］NOBLESSE F，DELHOMMEAU G，YANG C.Practical evaluation of steady flow due to a free-surface pressure patch［J］.Journal of Ship Research，2009，53(3):137-150.

［4］韓端鋒，黃德波.近水面潛體興波阻力的數(shù)值預(yù)報(bào)和收斂性分析［J］.船舶力學(xué)，2005，9(1):29-35.HAN Duanfeng，HUANG Debo.Numerical prediction and convergence analysis of wave resistance for submerged body near free surface［J］.Journal of Ship Mechanics，2005，9(1):29-35.

［5］NOBLESSE F.Alternative integral representations for the Green function of the theory of ship wave resistance［J］.Journal of Engineering Mathematics，1981，15(4):241-265.

［6］NOBLESSE F.The fundamental solution in the theory of steady motion of a ship［J］.Journal of Ship Research，1977，21(2):82-88.

［7］PONIZY B，BA M，GUILBAUD M，et al.Fast computation of the Green function for steady ship wave resistance［J］.Transactions on Modeling and Simulation，1993，6:89-96.

［8］PONIZY B，NOBLESSE F，BA M，et al.Numerical evaluation of free-surface Green functions［J］.Journal of Ship Research，1994，38(3):193-202.

［9］NEWMAN J N.Evaluation of the wave-resistance Green function:Part 1-The double integral［J］.Journal of Ship Research，1987，31(2):79-90.

［10］WANG H T，ROGERS J C W.Numerical evaluation of the complete wave-resistance Green’s function using Bessho’s approach［C］//Proceedings of the Fifth Conference on Numerical Ship Hydrodynamics.Washington，DC，USA，1990:133-144.

［11］URSELL F.Mathematical note on the fundamental solution(Kelvin source)in ship hydrodynamics［J］.Journal of Applied Mathematics，1984，32:335-351.

［12］NOBLESSE F，DELHOMMEAU G，HUANG Fuxin，et al.Practical mathematical representation of the flow due to a distribution of sources on a steadily advancing ship hull［J］.Journal of Engineering Mathematics，2011，71:367-392.

［13］BAAR J J M，PRICE W G.Evaluation of the wavelike disturbance in the Kelvin wave source potential［J］.Journal of Ship Research，1988，32(1):44-53.

［14］MARR G P，JACKSON P S.Some improvements and comparisons in the solution of the Neumann-Kelvin problem［J］.Journal of Ship Research，1999，43(3):170-179.

［15］MARR G P.An investigation of Neumann-Kelvin ship wave theory and its application to yacht design［D］.Auckland:University of Auckland，1995.

［16］XU Yong ，DONG Wencai.A fast integration method and characteristics research for 3-D translating-pulsating Green’s function of Havelock form in deep water［J］.China Ocean Engineering，2011，25(3):365-380.

［17］姚朝幫，董文才.開(kāi)爾文源格林函數(shù)數(shù)值積分方法研究［J］.上海交通大學(xué)學(xué)報(bào)，2014，48(1):98-105.YAO Chaobang，DONG Wencai.Evaluation method for Kelvin Source Green＇s function［J］.Journal of Shanghai Jiaotong University，2014，48(1):98-105.

［18］葉偉，陸鑫森.頻域有航速Green函數(shù)及其梯度的數(shù)值計(jì)算方法［J］.上海交通大學(xué)學(xué)報(bào)，1996，30(10):1-8.YE Wei，LU Xinsen.Numerical calculation of the Green’s function and its gradients with forward speed in the frequency domain［J］.Journal of Shanghai Jiaotong University，1996，30(10):1-8.