基于徑向基函數(shù)的自適應(yīng)網(wǎng)格方法

段獻葆， 黨 妍， 秦 玲

(西安理工大學(xué)理學(xué)院，西安 710048)

1 引言

許多物理和工程實際問題都可以用偏微分方程來描述，但是只有極少數(shù)的偏微分方程可以得到精確解，所以對于一般的偏微分方程，都是借助于數(shù)值方法求解.比較成熟的數(shù)值方法中大部分是依賴于網(wǎng)格的，如有限差分法、有限元法、有限體積等等.這些方法必須先生成網(wǎng)格后才能求解，網(wǎng)格質(zhì)量的好壞直接決定了最終數(shù)值求解的精度，而網(wǎng)格生成的預(yù)處理耗費時間太大，在求解區(qū)域不規(guī)則或維數(shù)較高時，這些方法都有一定的困難.另一方面，在應(yīng)用這些方法求解大變形、斷裂問題或高維問題，特別是非定常問題，很多時候在求解的過程中都需要網(wǎng)格重構(gòu)，這大大地增加了計算量，也嚴重降低了數(shù)值解的精度[1，2].為了解決這些問題，自適應(yīng)網(wǎng)格方法是一個很好的選擇，國內(nèi)外許多學(xué)者在這方面進行了大量地研究[3-6].

同樣是為了解決有限元方法等對網(wǎng)格依賴的問題，無網(wǎng)格方法近年來受到了很大的關(guān)注[7，8].這類方法試圖徹底地或部分地消除網(wǎng)格.相對于網(wǎng)格依賴的數(shù)值方法，無網(wǎng)格方法具有一些明顯的優(yōu)點，因此在其出現(xiàn)后得到了日益眾多的關(guān)注，并且一直是偏微分方程數(shù)值方法中的一個研究熱點.一般說來，無網(wǎng)格方法基于一系列節(jié)點進行函數(shù)插值，與有限元等網(wǎng)格依賴的方法相比，避免了網(wǎng)格劃分的預(yù)處理過程，也不會出現(xiàn)網(wǎng)格畸變的問題，對間斷問題、解的變化較大的問題等具有較好的優(yōu)勢.另一方面，無網(wǎng)格方法只需各個節(jié)點的獨立信息，而不需要單元之間的相互信息，數(shù)據(jù)結(jié)構(gòu)簡單.

利用徑向基函數(shù)(Radial Basis Function， RBF)求解偏微分方程的方法，是最近幾十年來備受關(guān)注的一類無網(wǎng)格方法[9].近年來，人們已經(jīng)用RBF 方法求解了各種線性和非線性的偏微分方程，并且取得了不錯的結(jié)果.其基本思想是在求解區(qū)域內(nèi)根據(jù)所求解的問題布置節(jié)點，然后在每個節(jié)點上構(gòu)造RBF，再將RBF 代入到所求解的偏微分方程，通過求解最終得到的代數(shù)方程獲得近似解.與網(wǎng)格依賴的方法不同，RBF 方法不需要任何網(wǎng)格，對支撐域和邊界沒有要求，只需要以節(jié)點為中心的子域覆蓋所求解區(qū)域即可.同時，RBF 與空間維數(shù)無關(guān)，僅依賴于節(jié)點間的距離，低維結(jié)果可以很容易推廣到高維問題.事實上，RBF 是通過定義在[0，+∞)上的一元函數(shù)φ 與Rd上的歐幾里德范數(shù)‖·‖來表示d 元函數(shù)φ(‖x-y‖)，其中x，y ∈Rd.因此用RBF 來處理多元問題具有效率高、儲存方便、運算簡單的特點.RBF 已廣泛應(yīng)用于計算幾何、微分方程數(shù)值解、神經(jīng)網(wǎng)絡(luò)等領(lǐng)域.近年來，國內(nèi)外學(xué)者對RBF 的理論與應(yīng)用進行了系統(tǒng)的研究，隨著RBF 理論的逐漸發(fā)展，目前已成為求解偏微分方程的一個強有力的工具[10-15].但是RBF 方法的缺點也是顯而易見的，如理論方面很不完善，在逼近過程中所得到矩陣方程的系數(shù)矩陣是否可逆沒有相關(guān)理論結(jié)果，也即數(shù)值解的唯一性還沒解決；隨著中心節(jié)點增加，需要求解一個很大的方程組，并且這個方程組經(jīng)常是病態(tài)的；尚未見到開源或商業(yè)化的軟件，后處理以及所得結(jié)果與其他軟件的接口比較困難等等.為了解決這些問題，已經(jīng)有學(xué)者考慮自適應(yīng)的RBF 方法[16，17]，以及把有限元方法和RBF 方法相耦合的方法[18].

為了發(fā)揮無網(wǎng)格方法和網(wǎng)格依賴方法各自的優(yōu)勢，我們提出了一種基于徑向基函數(shù)的自適應(yīng)網(wǎng)格方法.利用有限元方法等數(shù)值結(jié)果結(jié)合徑向基函數(shù)方法在求解區(qū)域內(nèi)進行自適應(yīng)配點，克服了有限元等方法計算中網(wǎng)格畸變和重新生成帶來的困難，所得結(jié)果可以用有限元等方法的后處理技術(shù)進行分析.數(shù)值結(jié)果說明，所提方法產(chǎn)生的網(wǎng)格可以根據(jù)解的變化情況進行網(wǎng)格自適應(yīng)，從而在保證相同數(shù)值求解精度的情況下可以極大地節(jié)省計算量.

2 徑向基函數(shù)

徑向基函數(shù)插值方法以其計算格式簡單、節(jié)點配置靈活、精度高的特點而成為研究多元逼近理論的有利工具，并已經(jīng)被應(yīng)用于科學(xué)計算模擬和實際工程問題中.

徑向函數(shù)滿足以下條件[19]：如果‖x1‖ = ‖x2‖，就有φ(x1) = φ(x2)的函數(shù)φ，也即，僅依賴r = ‖x‖的函數(shù)(其中‖·‖為Euclid 范數(shù)).RBF 就是這樣的函數(shù)空間：給定一個在定義域x ∈Rd上的一元函數(shù)φ:R+→R，所有形如Φ(x-c)=φ(‖x-c‖)及其線性組合張成的函數(shù)空間稱為由函數(shù)φ 導(dǎo)出的RBF 空間.在一定的條件下，只要取{xj}兩兩不同，{Φ(x-xj)}就是線性無關(guān)的從而形成徑向基函數(shù)空間中某子空間的一組基.當(dāng){xj}幾乎充滿R 時，{Φ(x-xj)}及其線性組合可以逼近幾乎任何函數(shù)[9].

1) Kriging 方法的Gauss 分布函數(shù)(無限光滑)：φ(r)=exp(-βr2)， β ＞0.

2) Duchon 的thin plate(薄板)樣條函數(shù)(分片光滑)：

3) Sobolev 樣條函數(shù)：φ(r)=Kβ-d/2(r)rβ-d/2，其中K 為MacDonald 函數(shù).

用RBF 求解偏微分方程的一般步驟如下：

對于偏微分方程

其中L 是偏微分算子，B 是邊界算子.在Ω 內(nèi)配置N 個離散的數(shù)據(jù)點r1，r2，··· ，rN，其中r1，r2，··· ，rl是內(nèi)部節(jié)點，而rl+1，rl+2，··· ，rN是邊界節(jié)點.設(shè)方程(2)的解u 可以用RBF 表示為

其中λ1，λ2，··· ，λN為待定系數(shù).由(2)式和(3)式可得

由(4)和(5)兩式得到矩陣方程

其中

若矩陣A 是非奇異的，矩陣方程有唯一解，只要從(6)中解出λ，就可以得到方程(2)的近似解uN.

3 自適應(yīng)網(wǎng)格方法

最早由Berger 和Oliger 于1984 年提出的自適應(yīng)網(wǎng)格方法[20]，是一種高效且準確的數(shù)值方法.該方法拋棄了等距均勻的網(wǎng)格，代之以能夠自動地適應(yīng)所研究問題中解的特征的疏密程度不均的網(wǎng)格.其網(wǎng)格結(jié)構(gòu)隨著計算過程的推進而不斷的動態(tài)改變，根據(jù)計算的實際需要以及問題的特性改變計算區(qū)域內(nèi)的網(wǎng)格結(jié)構(gòu)，在物理量變化比較劇烈的地方，例如：大變形、激波面、接觸間斷面和滑移面等，采用空間尺度較小的精細網(wǎng)格，在物理量變化緩慢的地方則采用空間尺度較大的粗網(wǎng)格，這樣在保持計算高效率的同時得到高精度的數(shù)值解.

這一節(jié)我們給出一種簡單、易于實現(xiàn)的自適應(yīng)網(wǎng)格方法.該方法基于插值來調(diào)整RBF 的中心和參數(shù)，從而改變RBF 節(jié)點的分布.

圖1 網(wǎng)格自適應(yīng)過程

對于二維問題來說，首先將求解區(qū)域進行有限元分割，然后進行有限元求解，與一維問題類似，用所得到的數(shù)值解可以求出方程(3)中的系數(shù)λi；接下來，可以用徑向基和λi得到每個初始單元的中點(重心)處的值，用這個值與有限元解在這點的插值得到誤差；誤差超過給定閾值的點將成為新的節(jié)點，并且利用這個節(jié)點把所在的單元進行剖分；同樣，若誤差低于給定閾值，所在單元的節(jié)點都將被移除.

注1由于RBF 方法隨著節(jié)點的增加計算量會顯著上升，在誤差大于給定閾值的單元可以多細化幾次，但也不需要細化太多次.從我們的計算結(jié)果來看，細化兩次就可以達到非常理想的效果.

注2同一個節(jié)點可能會同時作為細化和粗化單元的節(jié)點，我們是先進行細化，然后進行粗化.

注3另一種可以采用的方法是一次增加多個節(jié)點，但在這個過程中會增加計算量，特別是RBF 方法.

注4這里給出的自適應(yīng)過程非常容易實施，并且可以推廣到更高的維度.

對于某些偏微分方程問題來說，如果初始節(jié)點不是太多，上述過程還可以簡化.如對于如下具有齊次Dirichlet 邊值的Poisson 問題

方程組(6)變?yōu)?/p>

可以用RBF 方法直接求得系數(shù)λi.此時可以用如下公式

得到在節(jié)點{y1，y2，...，yM}處的誤差，其中M 為RBF 的節(jié)點數(shù).

4 數(shù)值算例

本算例考慮Burgers 方程的求解問題.Burgers 方程是偏微分方程中的一個非常重要的方程，廣泛應(yīng)用于各個領(lǐng)域，如流體力學(xué)、非線性、氣體動力學(xué)、交通流等等.由于Burgers 方程是一個非線性方程，只有在很少的一些特殊情況下可以得到精確解，在更一般情形下，連數(shù)值求解都存在很大的困難.

我們考慮如下非定常Burgers 方程

其中Ω ?R2，?Ω 是Ω 的邊界，u =(u1，u2)是流體的速度.ν =1/Re(Re-Reynolds 數(shù))是粘性系數(shù)，f 表示體積力，g 為已知函數(shù).

設(shè)V =H1(Ω)2，則(13)-(15)的弱形式為：求u ∈V，使得

設(shè)Th是Ω 的一個三角剖分(h 是剖分參數(shù))，Vh是有限元空間，Vh?V，則(16)的有限元解為：求uh∈Vh，使得

在下面的計算中，求解區(qū)域Ω=[0，1]×[0，1]， T =1，時間步長為Δt=0.01，

其中‖·‖表示歐氏范數(shù)，x =(x1，x2)， d =(0.3，0.3)， R=0.25.

本算例求解的過程中，當(dāng)誤差大于2×10-3時進行網(wǎng)格細化，小于2×10-6時進行網(wǎng)格粗化.粘性系數(shù)ν = 0.01.圖2 至圖4 分別給出了問題在時間t = 0， 0.5， 1 時的解.其中，圖2(a)、圖3(a)、圖4(a)為RBF 節(jié)點分布；圖2(b)、圖3(b)、圖4(b)為節(jié)點對應(yīng)的有限元網(wǎng)格；圖2(c)、圖3(c)、圖4(c)為網(wǎng)格對應(yīng)的有限元解.其中為了看得清楚，有限元解的圖片旋轉(zhuǎn)了一個角度.

圖2 t=0 時的解

圖3 t=0.5 時的解

圖4 t=1 時的解

利用本文所提算法，在t = 0.5 時所用節(jié)點數(shù)和有限元網(wǎng)格單元數(shù)分別為656 個和1286 個，所用時間375.8368 秒.作為比較，我們用傳統(tǒng)有限元方法在t=0.5 時對本問題進行了求解，所得結(jié)果如圖5 所示：用相同節(jié)點數(shù)的均勻(粗)網(wǎng)格對Burgers 方程進行了求解，所得結(jié)果如圖5(a)所示，可以看出所得結(jié)果非常粗糙，根本無法接受；用均勻加密(細)網(wǎng)格可以得到與圖3(c)接近精度的解，如圖5(b)所示，采用節(jié)點和有限元網(wǎng)格單元數(shù)分別為2704 個和5202 個，所用時間1184.7050 秒，大約是本文所提算法的3.15 倍.

圖5 傳統(tǒng)有限元方法得到的解

從以上算例可以看出，本文所提算法可以在保持網(wǎng)格數(shù)量減少或不變的前提下較好的提高問題的求解精度，從而節(jié)省了求解時間.

5 結(jié)論

本文給出了一種求解偏微分方程的自適應(yīng)網(wǎng)格方法，該方法把徑向基函數(shù)計算格式簡單、節(jié)點配置靈活的優(yōu)點與網(wǎng)格依賴方法的穩(wěn)健性很好地結(jié)合起來.該算法非常容易實施.數(shù)值算例也表明所提算法可以在解變化劇烈的區(qū)域加密網(wǎng)格，而在解變化平緩的地方粗化網(wǎng)格，從而可以保持較高精度的前提下減少或不增加計算量.我們用非定常Burgers 問題對算法進行的驗證說明所提算法可以高效地求解偏微分方程問題.