三維脈動源格林函數的一種簡化算法

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

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

0 引言

在水動力學研究中，經常采用格林函數建立方程求解速度勢函數，進而分析流體對結構的作用力以及結構的運動響應。針對不同的實際問題，需要采用不同的格林函數。很多學者從不同角度、采用不同的方法對各種形式格林函數的數值計算方法進行了研究，以期能夠更加快速、準確的計算格林函數。

許多學者對三維移動脈動源的數值計算方法進行了廣泛深入的研究［1～3］。關于頻域格林函數和時域格林函數的數值計算方法，也有不少學者進行了許多的研究工作［4～6］。然而，目前為止，針對三維脈動源格林函數的數值研究并不多。紀剛［7］和陳崧等［8］采用零階 Bessel函數的 Laplace變換，將無航速具有自由液面的三維脈動源格林函數變形為易于截斷處理的形式，實現了無窮積分向有限積分的轉變，并采用了樣條差值技術和自適應方法解決了高頻振蕩函數的數值積分的速度和精度問題。

本文主要針對三維脈動源，在分析其格林函數形式特點的基礎上，對三維脈動源格林函數的數值計算問題進行了簡化研究。

1 Bessel函數的三個積分

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

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

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

如果直接在這種格林函數表達形式的基礎上進行計算，過程比較繁瑣。因此，本文從上述格林函數所包含的三種Bessel函數積分形式入手，對這種三維脈動源格林函數進行分析和簡化，進而提出了一種簡化算法。

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

表1 系數值

圖1 積分曲線

2 積分曲線的擬合

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

2．1 擬合方法及結果

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

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

當對數函數的系數包含二次以上的項時，在［0．000 1，1)上會產生波動，所以選擇以關于x的一次函數作為對數函數的系數。此外分子部分還有x2項和x3項，而分母部分的最高次數為四次，這種形式既可以保證在區間［1，10］上會產生波動，又能使當x的值比較大時候函數以1/x的形式衰減。在式(9)的基礎上，經過比較采用式(7)作為目標函數，得到了表2的9個系數值。將系數值代入式(8)，所得結果如圖2所示。

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

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

經過多次試驗后發現，系數為9個時得到的結果最理想:表3。

表2 系數值

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

表3 優化系數值

根據表3所得的優化系數值繪制的曲線，較之未經優化的曲線在精度上有所提高。優化曲線與積分曲線的比較如圖3所示。

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

因此，最終的擬合函數確定為:

圖3 優化曲線與積分曲線的比較

2．2 誤差曲線及分析

在前面的工作里，雖然我們已經得到了精度比較好的結果，但是與原數據相比還有一定的誤差。出現這些誤差的根源，主要是簡化、截斷和舍入引起的。在最初的積分計算中，采用了式(6)的Hess-Smith近似逼近式，這是模型簡化誤差的主要來源。而在積分程序以及優化程序的運算過程中，不可避免的又引入了截斷誤差和舍入誤差。因此，式(10)只是一個基本上能滿足工程問題要求的近似式，它在每一個計算點的誤差可以用對應的曲線表示，如圖4所示，圖中是對應最優結果的誤差曲線。

圖4 誤差曲線

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

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

通過對圖4的研究，我們可以初步判斷這個要擬合出來的誤差函數應該是一個多項式、一個余弦函數和一個指數函數的乘積，其形式應大致為:

當然上式還有很多地方值得商榷，比如說x的次數、項數等，又或其中可能還有其他形式的函數可以替換，這都很有研究價值。我們只在式(11)的基礎上，使用可變誤差多面體法，得到了一個精度很一般的結果，其系數值見表4。

表4 系數值

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

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

3 三維脈動源格林函數的簡化

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

4 結論

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

［1］ 宗智，黃鼎良．三維移動脈動源速度勢的數值研究［J］．水動力學研究與進展，1991，(6):55-63．

［2］ 盧曉平，葉恒奎，張緯康，等．Haskind源格林函數的奇異性研究與數值積分方法［J］．水動力學研究與進展，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］ 戴愚志，余建星，郭海濤．時域有限水深格林函數及其導數的數值計算［J］．船舶力學，2006，10(1):15-21．

［6］ 劉昌鳳，勾瑩，滕斌．時域格林函數卷積的新算法［J］．水動力學研究與進展，2010，25(4):524-533．

［7］ 紀剛，張緯康，盧曉平．無航速具有自由液面的三維脈動源Green函數的數值積分方法［J］．海軍工程大學學報，2004，16(6):89-92．

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