開爾文源格林函數(shù)數(shù)值計算方法對比研究

姚朝幫，董文才

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

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

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

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

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

1 Kelvin源格林函數(shù)

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

G*對場點的偏導數(shù)的表達式可參見文獻［8］，其和G*的計算相比并未引入計算難度。

2 近場擾動項N(a)的計算方法

2.1 制表法

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

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

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

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

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

2.2 切比雪夫多項式擬合方法

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

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

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

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

式中:初值U0=

式中，s=

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

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

2.3 有理函數(shù)法

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

2.4 改進的梯形法

圖1 改進梯形法流程Fig.1 Solution strategy for improved trapezoidal rule

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

3 遠場波動項W(a)的計算方法

3.1 Bessel函數(shù)法

Bessho通過引入Dawson型積分，將開爾文源格林函數(shù)的波動項轉(zhuǎn)換成了Bessel函數(shù)的級數(shù)形式［13］為

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

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

3.2 變量代換法

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

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

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

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

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

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

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

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

需要指出的是，式(18)及(19)中N(a)的積分區(qū)間是(-∞，+∞)，不利于數(shù)值積分的實施，實際數(shù)值計算過程中需要對積分區(qū)間進行適當截斷，具體截斷的方法可參見文獻［17］。

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

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

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

圖3 移動源點與場點的位置Fig.3 The location of translating source and field points

4.1 N(a)計算結(jié)果分析及效率評估

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

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

表1 N(a)計算耗時對比Table 1 Computation time comparison of N(a)

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

4.2 W(a)計算結(jié)果分析及效率評估

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

表2 W(a)計算耗時對比Table 2 Computation time comparison of W(a)

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

4.3 G*的計算結(jié)果及驗證

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

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

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

5 結(jié)論

本文在分析開爾文源格林函數(shù)近場擾動項和遠場波動項計算方法的基礎(chǔ)上，對比分析了各種計算方法的計算精度與效率，通過研究得到以下結(jié)論:

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

2)對比分析了波動項N(a)的3種計算方法，它們的計算精度基本相當，“變量代換法”效率較高，采用Bessel函數(shù)級數(shù)展開的計算方法次之，Bessho型表達計算耗時較長。

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

［1］李子如，陳克強.用面元法預(yù)報船舶的興波阻力及波型［J］.武漢理工大學學報，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］田超，吳有生.水線源強消棄法在船舶定常興波計算中的應(yīng)用［J］.船舶力學，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ù)報和收斂性分析［J］.船舶力學，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］姚朝幫，董文才.開爾文源格林函數(shù)數(shù)值積分方法研究［J］.上海交通大學學報，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］.上海交通大學學報，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.