關(guān)于配點(diǎn)型無(wú)網(wǎng)格法邊界條件處理技術(shù)

張宏偉, 李美香, 李衛(wèi)國(guó)

(大連理工大學(xué)數(shù)學(xué)科學(xué)學(xué)院,遼寧大連 116024)

0 引 言

無(wú)網(wǎng)格方法是過(guò)去十多年興起的數(shù)值方法,該方法是利用一組散布在問(wèn)題域及其邊界上的節(jié)點(diǎn)表示該問(wèn)題域和其邊界,通過(guò)幾個(gè)互不相關(guān)節(jié)點(diǎn)上的值擬合出一個(gè)逼近函數(shù),該函數(shù)具有較好的光滑性而且導(dǎo)數(shù)連續(xù),這樣不僅擺脫了網(wǎng)格的束縛,避免了復(fù)雜的網(wǎng)格生成及重新劃分工作,而且提供了連續(xù)性好、形式靈活的形函數(shù).該方法在函數(shù)的逼近、邊界條件的引入和能量泛函的積分等方面具有一定的靈活性,此外,還具有精度高、前后處理比有限元簡(jiǎn)便等特性.因此,無(wú)網(wǎng)格法可以用來(lái)處理大量的基于網(wǎng)格單元的計(jì)算方法難以求解的問(wèn)題,具有有限元法和有限差分法不可比擬的優(yōu)點(diǎn).現(xiàn)有的無(wú)網(wǎng)格法基本上可以分為Galerkin型(如 EFGM[1]、RKPM[2]等)和配點(diǎn)型(如 SPH[3]、RBF[4]等)兩大類.Galerkin法具有良好的穩(wěn)定性和精確性,通過(guò)積分運(yùn)算降低了對(duì)形函數(shù)的連續(xù)性要求,可方便地施加導(dǎo)數(shù)邊界條件.但Galerkin法需要布置背景無(wú)網(wǎng)格進(jìn)行數(shù)值積分,計(jì)算量大.另外,Galerkin法要求形函數(shù)具有全局相容性,處理本質(zhì)邊界條件困難.配點(diǎn)型無(wú)網(wǎng)格法是真正的無(wú)網(wǎng)格法,它應(yīng)用于實(shí)際的關(guān)鍵問(wèn)題是導(dǎo)數(shù)邊界條件的存在.邊界條件對(duì)微分方程至關(guān)重要,但導(dǎo)數(shù)邊界條件才是造成基于任意節(jié)點(diǎn)的配點(diǎn)型無(wú)網(wǎng)格方法精度偏低且穩(wěn)定性差的真正原因.本文將介紹配點(diǎn)型無(wú)網(wǎng)格法及其特點(diǎn),對(duì)配點(diǎn)法中幾種流行的導(dǎo)數(shù)邊界條件處理技術(shù)進(jìn)行詳細(xì)的探討.著重利用有限差分法中的積分插值法,對(duì)一維問(wèn)題的導(dǎo)數(shù)邊界條件提出操作簡(jiǎn)單、計(jì)算精度高、不增加自由度的較理想的解決方案,并通過(guò)具體的數(shù)值例子說(shuō)明它的優(yōu)越性.

1 配點(diǎn)格式及特點(diǎn)

設(shè)在問(wèn)題域Ψ(其邊界為Γ)上未知函數(shù)u(x)滿足的微分方程為

和邊界條件

其中F和G是變量的微分算子.

配點(diǎn)法是使微分方程在問(wèn)題域上的一組特定的點(diǎn)上成立.設(shè)函數(shù)u(x)的無(wú)網(wǎng)格近似函數(shù)為uh(x)(各近似函數(shù)的構(gòu)造方法可參看文獻(xiàn)[5]),r1和r2分別為計(jì)算微分方程和邊界條件的配點(diǎn)個(gè)數(shù) ,且r=r1+r2,則有配點(diǎn)法的實(shí)質(zhì)是強(qiáng)迫殘量在r個(gè)配點(diǎn)上等于零.r大于或等于待求參量的個(gè)數(shù).配點(diǎn)型無(wú)網(wǎng)格法有以下特點(diǎn):

(1)離散微分方程的過(guò)程直接,實(shí)現(xiàn)其離散方程的算法簡(jiǎn)潔;

(2)是真正的無(wú)網(wǎng)格法,只需要在各節(jié)點(diǎn)處計(jì)算形函數(shù)及其導(dǎo)數(shù),不需要進(jìn)行數(shù)值積分,易于實(shí)現(xiàn),計(jì)算效率高;

(3)位移邊界條件很容易實(shí)現(xiàn),不需要特殊處理;

(4)需要在各節(jié)點(diǎn)處計(jì)算形函數(shù)的二階導(dǎo)數(shù);

(5)系數(shù)矩陣一般是稀疏非對(duì)稱的;

(6)在某些情況下,離散代數(shù)方程是病態(tài)的,其解不穩(wěn)定,且精度差.

2 導(dǎo)數(shù)邊界條件處理技術(shù)

配點(diǎn)法應(yīng)用于實(shí)際的關(guān)鍵問(wèn)題是導(dǎo)數(shù)邊界條件的存在.當(dāng)邊界條件均為Dirichlet型時(shí),用配點(diǎn)型無(wú)網(wǎng)格法可得到微分方程的精確結(jié)果,但存在任意導(dǎo)數(shù)邊界條件時(shí),其結(jié)果將顯著惡化且不穩(wěn)定.這是由于采用無(wú)網(wǎng)格近似函數(shù)進(jìn)行函數(shù)逼近時(shí),在邊界處對(duì)函數(shù)的二階導(dǎo)數(shù)逼近誤差很大[6].有許多用于在配點(diǎn)型無(wú)網(wǎng)格法中施加導(dǎo)數(shù)邊界條件的方法,其中常用的5種為

(1)直接配點(diǎn)法

通過(guò)簡(jiǎn)單的配點(diǎn)對(duì)導(dǎo)數(shù)邊界進(jìn)行離散而得到不同于控制系統(tǒng)方程的一組單獨(dú)方程,即對(duì)導(dǎo)數(shù)邊界條件不作特殊處理.

(2)虛點(diǎn)法

沿導(dǎo)數(shù)邊界在問(wèn)題域的外側(cè)增設(shè)一組虛點(diǎn).在此條件下對(duì)于每個(gè)導(dǎo)數(shù)邊界節(jié)點(diǎn)可得到兩個(gè)方程:一個(gè)方程對(duì)應(yīng)導(dǎo)數(shù)邊界條件,另一個(gè)對(duì)應(yīng)微分方程.

(3)Hermite型配點(diǎn)法

對(duì)導(dǎo)數(shù)邊界節(jié)點(diǎn)利用附加的導(dǎo)數(shù)變量施加導(dǎo)數(shù)邊界條件.

(4)規(guī)則網(wǎng)格法

該方法沿問(wèn)題域的導(dǎo)數(shù)邊界內(nèi)側(cè)布置一層或數(shù)層規(guī)則分布的節(jié)點(diǎn).對(duì)這些節(jié)點(diǎn)采用FDM的標(biāo)準(zhǔn)差分方案,即采用與標(biāo)準(zhǔn)FDM一樣的處理方法施加導(dǎo)數(shù)邊界條件.

(5)加密導(dǎo)數(shù)邊界附近的節(jié)點(diǎn)

本文針對(duì)求解Helmholtz問(wèn)題的基于點(diǎn)插值的配點(diǎn)型無(wú)網(wǎng)格法[6],對(duì)具有導(dǎo)數(shù)邊界條件的Helmholtz問(wèn)題采用積分插值的方法處理導(dǎo)數(shù)邊界條件.對(duì)3個(gè)節(jié)點(diǎn)的插值法構(gòu)造形函數(shù)的配點(diǎn)法,采用積分插值法1,其過(guò)程為設(shè)Helmholtz問(wèn)題的導(dǎo)數(shù)邊界條件

對(duì)u′(x1/2)采用中心差商采用左矩形公式得

于是有

那么對(duì)式(3)可以用下式離散:

其中u1、u2分別為u(x1)、u(x2)的近似值.

對(duì)邊界條件(4)在區(qū)間[xN-1/2,xN]上構(gòu)造積分守恒式,采用同樣的方法處理.

對(duì)于用5個(gè)節(jié)點(diǎn)插值構(gòu)造形函數(shù)的配點(diǎn)法,采用的積分插值法2,其過(guò)程如下:在[xN-4,xN]上Helmholtz方程滿足積分方程

從而有

即得

對(duì)其用New ton-Cotes公式近似,即

因此,導(dǎo)數(shù)邊界條件(4)可以近似表示為

其中dc為節(jié)點(diǎn)間距,對(duì)u′(xN-4)用已求得的包含節(jié)點(diǎn)xN-4的函數(shù)近似式uh(x)的導(dǎo)數(shù)近似.對(duì)邊界條件(3),在區(qū)間[x1,x3]上構(gòu)造積分方程,采用同樣的方法處理(但對(duì)散亂節(jié)點(diǎn)的情形,導(dǎo)數(shù)邊界節(jié)點(diǎn)要進(jìn)行規(guī)則化處理).

規(guī)則節(jié)點(diǎn)法是在導(dǎo)數(shù)邊界處安排3個(gè)規(guī)則分布的節(jié)點(diǎn)xN-2、xN-1、xN.使用如下有限差分算法計(jì)算導(dǎo)數(shù)邊界節(jié)點(diǎn)的一階導(dǎo)數(shù)

利用上式代替導(dǎo)數(shù)邊界條件得到離散代數(shù)方程組.

在Hermite型近似中,將增設(shè)導(dǎo)數(shù)變量作為附加的自由度,對(duì)于內(nèi)部節(jié)點(diǎn),如果其局部支持域中不含導(dǎo)數(shù)邊界節(jié)點(diǎn),則采用傳統(tǒng)的無(wú)網(wǎng)格形函數(shù),形成配點(diǎn)方程.如果其局部支持域中包含導(dǎo)數(shù)邊界節(jié)點(diǎn),則采用如下的Herm ite型近似式

式中:N為由Hermite型近似得到的形函數(shù)向量,NH(x)為與導(dǎo)數(shù)自由度u′N有關(guān)的形函數(shù).對(duì)采用基于三角函數(shù)點(diǎn)插值的配點(diǎn)型無(wú)網(wǎng)格法,Hermite型近似式的形成過(guò)程如下:

對(duì)局部支持域中含有導(dǎo)數(shù)邊界節(jié)點(diǎn)的計(jì)算點(diǎn)x處的函數(shù)近似為

如采用5個(gè)節(jié)點(diǎn)插值構(gòu)造形函數(shù),則為

為確定式中的系數(shù)a1、a2、a3、a4,可用通過(guò)支持域中的3個(gè)節(jié)點(diǎn)的函數(shù)值的插值以及等于導(dǎo)數(shù)邊界節(jié)點(diǎn)的導(dǎo)數(shù)值而確定.即對(duì)局部支持域中的3個(gè)節(jié)點(diǎn)(包括導(dǎo)數(shù)邊界節(jié)點(diǎn))可得到

對(duì)式(8)求導(dǎo)后可得到邊界方程

聯(lián)立式(8)、(9)求出系數(shù)a1、a2、a3、a4就可以得到式(7)的形式.

3 數(shù)值算例

下面通過(guò)數(shù)值算例揭示各種不同導(dǎo)數(shù)邊界條件處理技術(shù)對(duì)計(jì)算精度的影響.設(shè)帶有Neum ann條件的一維Helmholtz問(wèn)題為

為揭示導(dǎo)數(shù)邊界條件對(duì)計(jì)算精度的影響,平均相對(duì)誤差e定義如下:

其中N為節(jié)點(diǎn)數(shù);ui為理論解;uhi為數(shù)值解.相對(duì)數(shù)值誤差范數(shù)的h收斂率R(e)定義為

其中為平均節(jié)點(diǎn)間距.

表1、2列出了不同的插值方案和幾種不同的處理導(dǎo)數(shù)邊界條件的方法求解Helm ho ltz問(wèn)題所得到的誤差結(jié)果.各種不同方法的h收斂率情況如表3所示.

由表1、2、3可得到如下結(jié)論:

(1)如果Helmholtz問(wèn)題僅有Dirichlet邊界條件,用點(diǎn)插值配點(diǎn)法可得到非常理想的結(jié)果,其誤差很小.用11個(gè)規(guī)則節(jié)點(diǎn)離散問(wèn)題域,對(duì)三節(jié)點(diǎn)插值的模型e≈0.080600%,對(duì)五節(jié)點(diǎn)插值的模型e≈0.000027%.

(2)Neumann條件的存在將導(dǎo)致較大的計(jì)算誤差.如果不對(duì)Neum ann條件采用特殊的處理方法(直接配點(diǎn)法),對(duì)三節(jié)點(diǎn)插值的模型e≈1.580000%,對(duì)五節(jié)點(diǎn)插值的模型e≈0.000833%,即誤差分別放大約20和31倍.因此,直接配點(diǎn)法雖然既簡(jiǎn)單又直接,但是它僅對(duì)Dirichlet邊界條件的問(wèn)題有效,對(duì)存在Neum ann條件的問(wèn)題通常不穩(wěn)定或不精確.

(3)對(duì)Neumann條件進(jìn)行特殊處理有助于改善解的精度.

(4)Hermite型配點(diǎn)法和積分插值法對(duì)兩種模型均可得到精確而穩(wěn)定的解.對(duì)三節(jié)點(diǎn)插值的模型兩種方法誤差放大約為4倍.對(duì)五節(jié)點(diǎn)插值的模型兩種方法的數(shù)值誤差幾乎與僅有Dirichlet條件產(chǎn)生的誤差相等.但Hermite法增加了自由度從而增加了計(jì)算成本.

(5)虛點(diǎn)法對(duì)兩種模型的處理也能得到較好的結(jié)果.然而它使得方程數(shù)增加并需要額外的網(wǎng)格劃分.

(6)規(guī)則網(wǎng)格法對(duì)五節(jié)點(diǎn)插值模型的效果比直接配點(diǎn)法差.這是因?yàn)槭?8)采用接近導(dǎo)數(shù)邊界的3個(gè)節(jié)點(diǎn)處理其導(dǎo)數(shù)邊界,而其誤差大體上由邊界條件所引起的較大誤差所決定,盡管增加了內(nèi)節(jié)點(diǎn)的插值階數(shù),但是其精度將不會(huì)改善.

(7)雖然各種方法的精度各不相同,但它們的收斂率幾乎相等.三節(jié)點(diǎn)插值方案的收斂率約為2,五節(jié)點(diǎn)為4～5.

表1 三節(jié)點(diǎn)插值不同方法所得到的誤差結(jié)果Tab.1 The resu lt errors of differentmethods with 3 point interpo lation

表2 五節(jié)點(diǎn)插值不同方法所得到的誤差結(jié)果Tab.2 The resu lt errors of differentmethods with 5 point interpo lation

表3 不同方法所得的h收斂率Tab.3 The h rate of convergence of differentmethods

4 結(jié) 語(yǔ)

本文在對(duì)基于三角函數(shù)點(diǎn)插值的配點(diǎn)型無(wú)網(wǎng)格法解Helmholtz問(wèn)題研究的基礎(chǔ)上,分析了導(dǎo)數(shù)邊界條件是影響配點(diǎn)型無(wú)網(wǎng)格方法求解精度的真正元兇,對(duì)幾種流行的導(dǎo)數(shù)邊界條件處理技術(shù)進(jìn)行了詳細(xì)的探討.對(duì)一維問(wèn)題的導(dǎo)數(shù)邊界條件提出了操作簡(jiǎn)單、計(jì)算精度高,不增加自由度的較理想的解決方案——積分插值法.并通過(guò)具體的數(shù)值例子說(shuō)明了它的優(yōu)越性.

[1]BELYTSCHKO T,LU Y Y,GU L.Element free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2):229-256

[2]LIU W K,JUN S,ZHANG Y F.Reproducing kernel particlemethods[J].International Journal for Numerical M ethods in Fluids,1995,20(8-9):1081-1106

[3]LUCY L B.A numerical app roach to the testing o f the fission hypothesis[J].The Astronom ical Journal,1977,8(12):1013-1024

[4]CHEN W.Symmetric boundary knot method[J].Engineering Ana lysis with Boundary Elements,2002,26(6):333-343

[5]張 雄,劉 巖.無(wú)網(wǎng)格法[M].北京:清華大學(xué)出版社,2004

[6]LIU G R,GU Y T.無(wú)網(wǎng)格法理論及程序設(shè)計(jì)[M].3版.王建明,周學(xué)軍,譯.濟(jì)南:山東大學(xué)出版社,2007