基于徑向基函數(shù)逼近的波浪折射問題求解

許錫賓，束仲祎，徐績青,2

(1. 重慶交通大學 河海學院，重慶 400074；2. 重慶交通大學 國家內河航道整治工程技術研究中心 水利水運工程教育部重點實驗室，重慶 400074)

0 引 言

波浪在傳播過程中，因流場、水深和水位的不均勻變化，海洋情況極為復雜，波浪折射方程始終難以精確求解。近年來，研究波浪折射的方法一般分為3種：① 波射線法，該方法計算簡便，容易被工程所接受，但對復雜地形條件計算結果會出現(xiàn)較大誤差；② 緩坡方程，該方法適用于水下坡度較緩情況，但對緩坡方程來說，其波浪、海底高度等環(huán)境因素必須滿足某種條件[1]，其適用范圍有所欠缺；③ 有限差分法，該方法作為一種無網(wǎng)格法被用于流體力學方面[2-3]，但由于這種方法需要邊界條件，無邊界水域情況會降低其求解精度。

筆者通過對波浪折射方程分析，并根據(jù)近年來國內外廣泛應用的徑向基函數(shù)幾乎可逼近所有函數(shù)這一性質[4-6]，將其用于波浪折射方程的求解。通過求解一系列的非線性方程組，可計算出波浪折射的折射角，且求解精度高，能更好地模擬波浪折射問題。

1 徑向基函數(shù)方法及計算思路

徑向基函數(shù)利用一個影響范圍內每個點與其他所有點的距離來表達它們之間的關系，以徑向距離作為變量，是一種具有強大逼近能力的基函數(shù)。該函數(shù)具有形式簡單、各向同性、不依賴空間維數(shù)等特點。其基本形式如式(1)：

(1)

式中：y(r)為所用的插值函數(shù)；n為指所使用的徑向基函數(shù)總數(shù)目，當節(jié)點影響域與求解范圍相一致時，一般與節(jié)點數(shù)目相同，當邊界條件有多個時，還會加入輔助函數(shù)以提高求解精度；αi為與之對應的i節(jié)點的權重系數(shù)；║x-xi║為每個節(jié)點與第i節(jié)點之間距離。

由此可見，徑向基函數(shù)本質為一個一元函數(shù)，具有計算簡便，應用廣泛的優(yōu)點，能進行多階求導，幾乎可逼近所有函數(shù)，這一巨大優(yōu)點也是人們應用的重要原因。

筆者主要使用的徑向基函數(shù)有：Gauss函數(shù)和緊支柱正定徑向基函數(shù)[7]。Gauss函數(shù)基本形式為如式(2)：

φ(r)=e-c2r2

(2)

緊支柱正定徑向基函數(shù)[4]基本形式如式(3)：

(3)

將所求解的空間區(qū)域Ω用N個節(jié)點(xI,yI)(其中：I=1, 2, …,N)離散，函數(shù)α(x,y)在區(qū)域Ω中的近似函數(shù)αh(x,y)可由一組以各節(jié)點(xI,yI)為中心的徑向基函數(shù)φI(x,y)的線性組合表示，如式(4)：

(4)

式中：aJ為待定系數(shù)，a=[a1,a2,…,aN]T；Φ(t)=[Φ1(x,y),Φ2(x,y),…,ΦN(x,y)]T。

將空間區(qū)域Ω中的N個節(jié)點帶入式(4)中，可得到N個方程，聯(lián)立為線性方程組如式(5)：

Aa=α

α=[α1,α2,…,αN]T

(5)

將邊界條件帶入式(6)中，即可得到未知量系數(shù)矩陣a。

文中，根據(jù)文獻[8]提出的精密計算概念和標準，Gauss函數(shù)求逆矩陣A-1誤差大概是緊支柱正定徑向基函數(shù)的8倍。由此也可預測出以Gauss函數(shù)為基本函數(shù)的計算結果精度較低，計算結果明顯存在震蕩現(xiàn)象。

2 波浪折射方程

在實際海況中，由于波向線與海底等深線斜交，波浪傳播方向會發(fā)生變化，即產生波浪折射。在波浪由深海向近海岸傳播過程中，隨著水深變淺，波長、波高和波向等波要素都會發(fā)生變化，但是波周期始終不變。忽略風速、風阻、海底摩阻等影響因素，將波浪場視為一個穩(wěn)定的場[9-11]，其計算如式(6)：

(6)

式中：θ為波浪相位函數(shù)；(x,y)為波浪場中任意一點坐標；k為波浪的波數(shù)。

在以往對波浪折射問題的計算中，一般都是取邊界點的波浪折射角α=0作為邊界條件。但該計算方法忽略了邊界導數(shù)條件，邊界點并不符合式(5)，所求出的解并不準確。筆者對于邊界點，取α=0。則式(6)可以化簡為如式(7)：

(7)

3 算例分析

現(xiàn)有一片水域如圖1。水底等高線為同心圓，圓心處高程與靜水平面相水平。水域最大半徑為3 km，水域內最大水深為-30 m，即最外部水域的水底高程為30 m。現(xiàn)需得出波浪進入此區(qū)域1.2 km范圍內折射角變化。

由于水域對稱分布，故取如圖1中水域的一半作為計算區(qū)域。為保證邊界部分計算精度能達到要求，故取x方向為-3 000 m～-1 500 m；y方向為0～-3 000 m；最后計算結果取1.2 km以內。考慮到編程方便和計算精度，計算區(qū)域內取點采用等間距布點，布點間距為25 m。將徑向基函數(shù)帶入波浪折射方程及邊界條件方程中，可得到微分方程組，再帶入邊界條件即可求解。

3.1 算例1

將式(3)作為基本函數(shù)，不考慮輔助函數(shù)，邊界條件取邊界點為0，不考慮導數(shù)邊界條件，計算結果如圖2。

此計算結果在邊界點位置折射角為0，在邊界處存在震蕩。由于不考慮導數(shù)邊界條件，邊界處的點不一定滿足波浪折射方程導數(shù)邊界條件，結果存在誤差。故此算例證明了導數(shù)邊界條件的重要性，后續(xù)算例都會考慮導數(shù)邊界條件。

3.2 算例2

將式(3)作為基本函數(shù)和輔助函數(shù)，邊界條件取邊界點為0和導數(shù)邊界條件，計算結果如圖3。

此計算結果滿足邊界點位置折射角為0，同時也滿足波浪折射方程導數(shù)形式，且計算出的曲面形態(tài)較好，無明顯震蕩現(xiàn)象。該計算方法計算出來的結果是符合波浪折射方程的，為筆者所建議采用。

3.3 算例3

將式(3)作為基本函數(shù)和輔助函數(shù)，邊界條件取邊界點為0和導數(shù)邊界條件，計算結果如圖4。

由圖4可看出，由于Gauss函數(shù)求逆矩陣A-1誤差較大，計算結果震蕩現(xiàn)象比較嚴重，存在較大誤差，因此這種Gauss函數(shù)作為基本函數(shù)形式不被采用。

3.4 算例4

將式(3)作為基本函數(shù)，式(3)的導數(shù)作為輔助函數(shù)，邊界條件取邊界點為0和導數(shù)邊界條件，計算結果如圖5。

此計算結果在局部位置起伏，存在震蕩現(xiàn)象，雖然不是很明顯，但對計算結果仍然存在影響。計算結果與算例2相比，計算精度有所降低，由此證明該算例所取的基本函數(shù)和輔助函數(shù)對計算結果存在影響，不建議使用。

4 對比分析

由于算例1中沒有考慮導數(shù)邊界條件，將其邊界點計算結果帶入式(7)中，求得誤差如圖6。

為對比分析其他算例計算結果，選擇一處相同的計算范圍(x軸坐標為-1 940～-2 060；y軸坐標為-1 600～-1 850)，分別截取各算例計算結果局部放大圖，如圖7。

由理論分析角度來看：在此計算范圍中，折射角方向應向右偏轉。由圖7(a)可看出：算例2計算結果符合理論分析，效果較好；從圖7(b)中可看出：算例3計算結果存在較明顯震蕩現(xiàn)象，計算結果誤差較大；從圖7(c)中可看出：算例4計算結果趨近于0，即在此范圍內其折射角為0，顯然不符合理論分析的結果。故筆者選取算例2計算為最終結果，將其與Snell折射定律進行對比分析。

取x軸坐標為-2 500～-1 500 m，間隔取100 m，y軸坐標為-1 500 m，共11個點，運用Snell折射定律計算，如式(8)：

(8)

式中：αi為折射角；ci為折射點處的波速；α0為入射角；c0為入射波速；k為折射點處波數(shù)；d為折射點處水深。

圖8為Snell定律和算例計算結果的對比。由圖8可看出：由于Snell折射定律是特殊情況下的波浪折射，它要求水域內各點水深相同，等深線平行，這種情況與算例實際情況不符合，故其計算結果呈線性。

5 結 論

筆者通過徑向基函數(shù)和等間距布點方法計算波浪折射時存在的問題。在運用Gauss函數(shù)時參考了緊支柱思想，將其中的參數(shù)c取值與全局支撐域相關聯(lián)，為今后Gauss函數(shù)運用提供了一種參數(shù)的取值方法。但Gauss函數(shù)的求逆矩陣A-1誤差較大，效果上存在震蕩現(xiàn)象，在使用上仍存在一定問題。故筆者的基本函數(shù)取緊支柱正定徑向基函數(shù)和輔助函數(shù)取式(3)。

筆者提供了一種計算波浪折射問題方法。在邊界點上符合導數(shù)邊界條件，在理論方面符合波浪折射原理，計算結果更精確，方法更合理。由于Snell折射定律的適用條件特殊，不適用于復雜水深下的海洋環(huán)境；而筆者提出的方法則適用于水深變化的波浪情況，為求解波浪折射數(shù)值解提供了一種新思路。