三維脈動(dòng)源格林函數(shù)的一種簡(jiǎn)化算法

劉 磊，賈敬蓓，劉大路，黃亞南

(1．中國(guó)船級(jí)社天津分社，天津300457;2．大連海洋大學(xué)航海與船舶工程學(xué)院，遼寧大連116023)

0 引言

在水動(dòng)力學(xué)研究中，經(jīng)常采用格林函數(shù)建立方程求解速度勢(shì)函數(shù)，進(jìn)而分析流體對(duì)結(jié)構(gòu)的作用力以及結(jié)構(gòu)的運(yùn)動(dòng)響應(yīng)。針對(duì)不同的實(shí)際問(wèn)題，需要采用不同的格林函數(shù)。很多學(xué)者從不同角度、采用不同的方法對(duì)各種形式格林函數(shù)的數(shù)值計(jì)算方法進(jìn)行了研究，以期能夠更加快速、準(zhǔn)確的計(jì)算格林函數(shù)。

許多學(xué)者對(duì)三維移動(dòng)脈動(dòng)源的數(shù)值計(jì)算方法進(jìn)行了廣泛深入的研究［1～3］。關(guān)于頻域格林函數(shù)和時(shí)域格林函數(shù)的數(shù)值計(jì)算方法，也有不少學(xué)者進(jìn)行了許多的研究工作［4～6］。然而，目前為止，針對(duì)三維脈動(dòng)源格林函數(shù)的數(shù)值研究并不多。紀(jì)剛［7］和陳崧等［8］采用零階 Bessel函數(shù)的 Laplace變換，將無(wú)航速具有自由液面的三維脈動(dòng)源格林函數(shù)變形為易于截?cái)嗵幚淼男问?，?shí)現(xiàn)了無(wú)窮積分向有限積分的轉(zhuǎn)變，并采用了樣條差值技術(shù)和自適應(yīng)方法解決了高頻振蕩函數(shù)的數(shù)值積分的速度和精度問(wèn)題。

本文主要針對(duì)三維脈動(dòng)源，在分析其格林函數(shù)形式特點(diǎn)的基礎(chǔ)上，對(duì)三維脈動(dòng)源格林函數(shù)的數(shù)值計(jì)算問(wèn)題進(jìn)行了簡(jiǎn)化研究。

1 Bessel函數(shù)的三個(gè)積分

在 John V．Wehausen和 Edmund V．Laitone于1960年發(fā)表的文章《Surface Waves》中，給出了三維脈動(dòng)源的一種格林函數(shù)形式:式中:k為波數(shù)空間的自變量;v為無(wú)限水深波數(shù)，，其中，ω 為波動(dòng)角速度，g為重力加速度;J0(kR)和J0(vR)分別是kR和vR的零階Bessel函數(shù);r為場(chǎng)點(diǎn)和源點(diǎn)之間的距離，r=R為場(chǎng)點(diǎn)和源點(diǎn)之間的水平距離y，z)是場(chǎng)點(diǎn)坐標(biāo)，(ξ，η，ζ)是源點(diǎn)坐標(biāo)。

在靠近自由表面的位置z=ζ=0，因此式(1)可以變化為下面形式:

由式(2)可以看出，式(1)最終可以歸結(jié)為求以下三種形式的Bessel函數(shù)積分，即:

如果直接在這種格林函數(shù)表達(dá)形式的基礎(chǔ)上進(jìn)行計(jì)算，過(guò)程比較繁瑣。因此，本文從上述格林函數(shù)所包含的三種Bessel函數(shù)積分形式入手，對(duì)這種三維脈動(dòng)源格林函數(shù)進(jìn)行分析和簡(jiǎn)化，進(jìn)而提出了一種簡(jiǎn)化算法。

編制FORTRAN程序即可求得使計(jì)算量比較少，又能在μ的取值比較小的時(shí)候滿足精度要求。以μ為自變量，積分值為因變量，得到如圖1的曲線。

表1 系數(shù)值

圖1 積分曲線

2 積分曲線的擬合

在得到圖1的曲線后，就期望通過(guò)已有的數(shù)據(jù)獲得一個(gè)函數(shù)f(μ)，要求這個(gè)函數(shù)可以比較好的符合得到的曲線。那么，在以后的計(jì)算中就可以在誤差允許的范圍內(nèi)使用擬合出的函數(shù)來(lái)求出需要的數(shù)值。除了在精度上對(duì)擬合函數(shù)有要求外，還要求所得函數(shù)的形式相對(duì)簡(jiǎn)單。相反如果函數(shù)形式很復(fù)雜，計(jì)算步驟很繁瑣甚至比直接計(jì)算還費(fèi)時(shí)，就失去了進(jìn)行這項(xiàng)工作的意義，也與初衷不符。

2．1 擬合方法及結(jié)果

擬合函數(shù)的優(yōu)化方法有很多中可供選擇，但是為了使得到的函數(shù)形式盡量簡(jiǎn)單，所以選擇使用G-A算法和可變誤差多面體算法。這兩種算法的缺點(diǎn)在于需要先大致確定出函數(shù)的形式，然后才能通過(guò)程序求出各項(xiàng)系數(shù)的值。優(yōu)點(diǎn)是操作比較簡(jiǎn)單，而且在函數(shù)形式相對(duì)簡(jiǎn)單的前提下精度基本符合要求。前期使用的是G-A算法，但是所得到的結(jié)果并不理想。而使用可變誤差多面體算法，在初始函數(shù)形式以及系數(shù)個(gè)數(shù)相同的前提下，后者明顯比G-A算法的精度要好。這并不是說(shuō)G-A算法不好，而是說(shuō)明在這個(gè)函數(shù)的擬合過(guò)程中，可變誤差多面體法更合適?？勺冋`差多面體法的計(jì)算過(guò)程實(shí)際上是一個(gè)可變多面體的尋找過(guò)程，也就是通過(guò)反射、膨脹、壓縮和收縮來(lái)改變多面體。但是在尋找過(guò)程中得到的每一個(gè)點(diǎn)若不在近似能行區(qū)域內(nèi)，就要經(jīng)過(guò)一個(gè)恢復(fù)過(guò)程，用過(guò)程中得到的一個(gè)屬于近似能行區(qū)域中的點(diǎn)來(lái)代換。這樣可以保證可變多面體算法是在能行區(qū)域中進(jìn)行的。

通過(guò)在［0．000 1，100］上得到的曲線，即圖 1，可以看到:在［0．000 1，1)上，函數(shù)基本以對(duì)數(shù)形式衰減;在(10，100］上，函數(shù)則以冪函數(shù)1/x的形式衰減;在［1，10］上，函數(shù)出現(xiàn)了波動(dòng)。由此推斷，該擬合函數(shù)是對(duì)數(shù)函數(shù)和冪函數(shù)的復(fù)合函數(shù)。確定函數(shù)的形式后，系數(shù)的個(gè)數(shù)由少到多進(jìn)行試驗(yàn)，嘗試的范圍最少為3個(gè)，最多為15個(gè)。在這個(gè)過(guò)程中，除了要改變系數(shù)的個(gè)數(shù)，還要對(duì)給出的函數(shù)形式及目標(biāo)函數(shù)進(jìn)行調(diào)整。函數(shù)的形式是基于對(duì)數(shù)函數(shù)和冪函數(shù)的基本性質(zhì)，而在程序中目標(biāo)函數(shù)有2種形式可供選擇。

當(dāng)對(duì)數(shù)函數(shù)的系數(shù)包含二次以上的項(xiàng)時(shí)，在［0．000 1，1)上會(huì)產(chǎn)生波動(dòng)，所以選擇以關(guān)于x的一次函數(shù)作為對(duì)數(shù)函數(shù)的系數(shù)。此外分子部分還有x2項(xiàng)和x3項(xiàng)，而分母部分的最高次數(shù)為四次，這種形式既可以保證在區(qū)間［1，10］上會(huì)產(chǎn)生波動(dòng)，又能使當(dāng)x的值比較大時(shí)候函數(shù)以1/x的形式衰減。在式(9)的基礎(chǔ)上，經(jīng)過(guò)比較采用式(7)作為目標(biāo)函數(shù)，得到了表2的9個(gè)系數(shù)值。將系數(shù)值代入式(8)，所得結(jié)果如圖2所示。

由圖2可以看出，在區(qū)間［10，100］上2條曲線符合的比較好。而在區(qū)間［1，10］上，尤其是波動(dòng)區(qū)域的附近，2條曲線出入比較大。此外，在區(qū)間［0．000 1，1)上2條曲線也有一定的差別。因此，為了達(dá)到更好的精度以滿足工程需要，決定對(duì)得到的9個(gè)系數(shù)值進(jìn)行多次優(yōu)化。以表2中9個(gè)數(shù)值為初始值，仍然使用可變誤差多面體算法進(jìn)行優(yōu)化。優(yōu)化的次數(shù)則以擬合的精度為準(zhǔn)，一共進(jìn)行了5次優(yōu)化，其中以第3次優(yōu)化后的結(jié)果最好。其系數(shù)值見(jiàn)

式中:Λ(μi)為原數(shù)據(jù)值，f(μi)為擬合函數(shù)值。Sum的值越小，結(jié)果就越好。

經(jīng)過(guò)多次試驗(yàn)后發(fā)現(xiàn)，系數(shù)為9個(gè)時(shí)得到的結(jié)果最理想:表3。

表2 系數(shù)值

圖2 擬合曲線與積分曲線的比較

表3 優(yōu)化系數(shù)值

根據(jù)表3所得的優(yōu)化系數(shù)值繪制的曲線，較之未經(jīng)優(yōu)化的曲線在精度上有所提高。優(yōu)化曲線與積分曲線的比較如圖3所示。

與最初擬合函數(shù)的曲線相比，優(yōu)化后的曲線無(wú)論在原本較好的區(qū)間(10，100］上，還是在出入比較大的區(qū)間［0．000 1，1)上，甚至在誤差最大的波動(dòng)區(qū)域，都符合的比較好了。

因此，最終的擬合函數(shù)確定為:

圖3 優(yōu)化曲線與積分曲線的比較

2．2 誤差曲線及分析

在前面的工作里，雖然我們已經(jīng)得到了精度比較好的結(jié)果，但是與原數(shù)據(jù)相比還有一定的誤差。出現(xiàn)這些誤差的根源，主要是簡(jiǎn)化、截?cái)嗪蜕崛胍鸬?。在最初的積分計(jì)算中，采用了式(6)的Hess-Smith近似逼近式，這是模型簡(jiǎn)化誤差的主要來(lái)源。而在積分程序以及優(yōu)化程序的運(yùn)算過(guò)程中，不可避免的又引入了截?cái)嗾`差和舍入誤差。因此，式(10)只是一個(gè)基本上能滿足工程問(wèn)題要求的近似式，它在每一個(gè)計(jì)算點(diǎn)的誤差可以用對(duì)應(yīng)的曲線表示，如圖4所示，圖中是對(duì)應(yīng)最優(yōu)結(jié)果的誤差曲線。

圖4 誤差曲線

從圖 4 得知，在區(qū)間［0．000 1，18］上誤差變化比較大，其值為|ε|＜0．3;而在區(qū)間(18，100］上，誤差值|ε|＜0．05。這樣可以計(jì)算出擬合函數(shù)的相對(duì)誤差|η|＜6%。這個(gè)誤差值在所有四次優(yōu)化結(jié)果中是最小的。初始擬合函數(shù)的誤差與優(yōu)化后的誤差相比是比較大的，而只有第三次優(yōu)化后的誤差亦即最佳擬合函數(shù)的誤差可以控制在0．3以內(nèi)。

對(duì)于普通的工程計(jì)算，這個(gè)誤差已經(jīng)可以滿足要求。但是，為了可以更加精確，我們希望可以將誤差曲線也擬合出來(lái)，這樣在式(10)后只要減掉誤差函數(shù)，就可以使擬合函數(shù)的誤差進(jìn)一步降低，這在理論上是完全可以實(shí)現(xiàn)的。

通過(guò)對(duì)圖4的研究，我們可以初步判斷這個(gè)要擬合出來(lái)的誤差函數(shù)應(yīng)該是一個(gè)多項(xiàng)式、一個(gè)余弦函數(shù)和一個(gè)指數(shù)函數(shù)的乘積，其形式應(yīng)大致為:

當(dāng)然上式還有很多地方值得商榷，比如說(shuō)x的次數(shù)、項(xiàng)數(shù)等，又或其中可能還有其他形式的函數(shù)可以替換，這都很有研究?jī)r(jià)值。我們只在式(11)的基礎(chǔ)上，使用可變誤差多面體法，得到了一個(gè)精度很一般的結(jié)果，其系數(shù)值見(jiàn)表4。

表4 系數(shù)值

由表4中的數(shù)值繪制出的誤差擬合函數(shù)的曲線如圖5所示。從圖5中可以看出，精度是比較差的，但是擬合曲線的大體趨勢(shì)是對(duì)的，這說(shuō)明擬合函數(shù)的形式基本正確。

圖5 擬合曲線與誤差曲線的比較

3 三維脈動(dòng)源格林函數(shù)的簡(jiǎn)化

在靠近自由表面的位置z=ζ=0，則可得如下形式:

4 結(jié)論

為了簡(jiǎn)化計(jì)算，本文提出了擬和簡(jiǎn)單函數(shù)的思想。采用可變誤差多面體算法來(lái)擬和形式復(fù)雜的函數(shù)，將一個(gè)表示三維脈動(dòng)源的格林函數(shù)簡(jiǎn)化為形式簡(jiǎn)單計(jì)算方便的復(fù)合函數(shù)。在滿足精度要求的前提下，大大減少了計(jì)算量。盡管原函數(shù)的擬和曲線比較理想，但誤差曲線的擬和還有 一些不足之處。不過(guò)這項(xiàng)探索性的研究說(shuō)明這種思想是可行的，為今后的工作提供了可以使用的簡(jiǎn)單函數(shù)。

［1］ 宗智，黃鼎良．三維移動(dòng)脈動(dòng)源速度勢(shì)的數(shù)值研究［J］．水動(dòng)力學(xué)研究與進(jìn)展，1991，(6):55-63．

［2］ 盧曉平，葉恒奎，張緯康，等．Haskind源格林函數(shù)的奇異性研究與數(shù)值積分方法［J］．水動(dòng)力學(xué)研究與進(jìn)展，1999，14(4):444-452．

［3］ XU Y，DONG W C．Study on characteristics of 3-D translating-pulsating source Green Function of deep-water Havelock form and its fast integration method［J］．China Ocean Engineering，2011，25(3):365-380．

［4］ YAO X L，SUN S L，WANG S P，YANG S T．The research on the highly efficient calculation method of 3-D frequency-domain Green function［J］．Journal of Marine Science and Application，2009，8:196-203．

［5］ 戴愚志，余建星，郭海濤．時(shí)域有限水深格林函數(shù)及其導(dǎo)數(shù)的數(shù)值計(jì)算［J］．船舶力學(xué)，2006，10(1):15-21．

［6］ 劉昌鳳，勾瑩，滕斌．時(shí)域格林函數(shù)卷積的新算法［J］．水動(dòng)力學(xué)研究與進(jìn)展，2010，25(4):524-533．

［7］ 紀(jì)剛，張緯康，盧曉平．無(wú)航速具有自由液面的三維脈動(dòng)源Green函數(shù)的數(shù)值積分方法［J］．海軍工程大學(xué)學(xué)報(bào)，2004，16(6):89-92．

［8］陳崧，趙留平，侯磊．無(wú)航速具有自由液面的三維Green脈動(dòng)源函數(shù)的數(shù)值積分方法［J］．中國(guó)艦船研究，2006，1(2):50-52．