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

2020-11-04 03:06段獻(xiàn)葆

工程數(shù)學(xué)學(xué)報(bào) 2020年5期

段獻(xiàn)葆，黨妍，秦玲

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

1 引言

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

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

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

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

2 徑向基函數(shù)

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

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

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

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

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

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

對(duì)于偏微分方程

其中L 是偏微分算子，B 是邊界算子.在Ω 內(nèi)配置N 個(gè)離散的數(shù)據(jù)點(diǎn)r1，r2，··· ，rN，其中r1，r2，··· ，rl是內(nèi)部節(jié)點(diǎn)，而rl+1，rl+2，··· ，rN是邊界節(jié)點(diǎn).設(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]，是一種高效且準(zhǔn)確的數(shù)值方法.該方法拋棄了等距均勻的網(wǎng)格，代之以能夠自動(dòng)地適應(yīng)所研究問(wèn)題中解的特征的疏密程度不均的網(wǎng)格.其網(wǎng)格結(jié)構(gòu)隨著計(jì)算過(guò)程的推進(jìn)而不斷的動(dòng)態(tài)改變，根據(jù)計(jì)算的實(shí)際需要以及問(wèn)題的特性改變計(jì)算區(qū)域內(nèi)的網(wǎng)格結(jié)構(gòu)，在物理量變化比較劇烈的地方，例如：大變形、激波面、接觸間斷面和滑移面等，采用空間尺度較小的精細(xì)網(wǎng)格，在物理量變化緩慢的地方則采用空間尺度較大的粗網(wǎng)格，這樣在保持計(jì)算高效率的同時(shí)得到高精度的數(shù)值解.

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

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

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

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

注2同一個(gè)節(jié)點(diǎn)可能會(huì)同時(shí)作為細(xì)化和粗化單元的節(jié)點(diǎn)，我們是先進(jìn)行細(xì)化，然后進(jìn)行粗化.

注3另一種可以采用的方法是一次增加多個(gè)節(jié)點(diǎn)，但在這個(gè)過(guò)程中會(huì)增加計(jì)算量，特別是RBF 方法.

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

對(duì)于某些偏微分方程問(wèn)題來(lái)說(shuō)，如果初始節(jié)點(diǎn)不是太多，上述過(guò)程還可以簡(jiǎn)化.如對(duì)于如下具有齊次Dirichlet 邊值的Poisson 問(wèn)題

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

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

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

4 數(shù)值算例

本算例考慮Burgers 方程的求解問(wèn)題.Burgers 方程是偏微分方程中的一個(gè)非常重要的方程，廣泛應(yīng)用于各個(gè)領(lǐng)域，如流體力學(xué)、非線(xiàn)性、氣體動(dòng)力學(xué)、交通流等等.由于Burgers 方程是一個(gè)非線(xiàn)性方程，只有在很少的一些特殊情況下可以得到精確解，在更一般情形下，連數(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是Ω 的一個(gè)三角剖分(h 是剖分參數(shù))，Vh是有限元空間，Vh?V，則(16)的有限元解為：求uh∈Vh，使得

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

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

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

圖2 t=0 時(shí)的解

圖3 t=0.5 時(shí)的解

圖4 t=1 時(shí)的解

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

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

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

5 結(jié)論