無(wú)單元Galerkin方法數(shù)值求解三維泊松方程

付 偉，解紅霞

(1.太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院,山西 太原 030024；2.太原學(xué)院 應(yīng)用數(shù)學(xué)系，山西 太原 030032)

無(wú)網(wǎng)格方法是近20年來(lái)熱門的一種數(shù)值計(jì)算方法,其中一種方法是無(wú)單元Galerkin方法[1-2],簡(jiǎn)稱EFG方法。無(wú)網(wǎng)格方法拋開(kāi)了網(wǎng)格的限制,這是它較有限元算法的優(yōu)勢(shì)。

無(wú)網(wǎng)格FEM法[3]是Idelsohn等人于2003年提出的,融合了無(wú)網(wǎng)格法與FEM法的優(yōu)點(diǎn)的一類新的數(shù)值方法。

2005年,程玉民等[4]提出了復(fù)變量MLS法。與較傳統(tǒng)的MLS法相比,該方法減少了函數(shù)中的變量,有效實(shí)現(xiàn)了降維,節(jié)約了計(jì)算時(shí)間,并且有更高的計(jì)算精度。

程玉民課題組于2009年提出了二維彈性動(dòng)力學(xué)的復(fù)變量邊界無(wú)單元(CVBEFM)法[5]。與EFG法相結(jié)合,先后于2011、2012年提出了針對(duì)二維彈性、彈塑性問(wèn)題的復(fù)變量EFG法[6]。MLS法的形函數(shù)不具Kronecker δ函數(shù)特性,采取措施來(lái)施加邊界條件就會(huì)使Galerkin弱形式更加復(fù)雜。

前人用EFG方法數(shù)值求解二維泊松方程。本文在此基礎(chǔ)上,推導(dǎo)了無(wú)單元Galerkin方法數(shù)值求解三維泊松方程時(shí)的計(jì)算公式,改進(jìn)算法流程,對(duì)比解析解和數(shù)值解,說(shuō)明該方法有較好的精度。

1 無(wú)單元Galerkin方法

1.1 等效積分弱形式

這里引入一個(gè)泊松方程,來(lái)探討無(wú)單元Galerkin方法。

Poisson方程:

u(x)-b(x)=0(在Ω內(nèi))

(1)

邊界條件:

(2)

(3)

式(1)、式(2)、式(3)的等效積分弱形式為:

(4)

其中,

(5)

已知

(6)

其中n是x點(diǎn)區(qū)域內(nèi)的節(jié)點(diǎn)個(gè)數(shù),Φ(x)函數(shù)，

Φ(x)=(φ1(x),φ2(x),…,φn(x))=pTA-1(x)B(x)

(7)

u(x)=(u1,u2,…,un)T=(u(x1),u(x2),…,u(xn))T

(8)

(9)

其中

B(x)=(B1(x),B2(x),…,Bn(x))

(10)

(11)

將(6)和(9)帶入(4)得到:

(12)

然后我們對(duì)(1.12)的各項(xiàng)進(jìn)行逐一積分,首先對(duì)第一項(xiàng)進(jìn)行積分:

(13)

其中，

(14)

K=[KIJ]為nt×nt階的矩陣,

(15)

下面是對(duì)(12)第二項(xiàng)進(jìn)行積分:

(16)

其中F(1)是給定源函數(shù)所引起的荷載列陣

F(1)=(f(1)(x1),f(1)(x2),…,f(1)(xn))T

(17)

(18)

最后對(duì)(12)的第三項(xiàng)進(jìn)行積分:

(19)

F(2)=(f(2)(x1),f(2)(x2),…,f(2)(xn))T

(20)

(21)

將式(1.13),(1.16)和(1.19)代入式(1.12)中,得到:

δuT·K·u-δuTF(1)-δuTF(2)=0

(22)

即

δuT·(Ku-F)=0

(23)

其中

F=F(1)+F(2)

(24)

由于δuT的任意性,使我們可以得到線性方程:

Ku=F

(25)

1.2 邊界條件的處理

根據(jù)式(1)、式(2)和式(3),引入罰因子我們可以得到:

(26)

其中α=(α1,α2,…αi),這里α是罰因子的對(duì)角矩陣,針對(duì)于二維問(wèn)題i=2,針對(duì)于三維問(wèn)題時(shí)i=3,這個(gè)罰因子αk(k=1,2,…i),可能是一個(gè)坐標(biāo)函數(shù),通常情況它是由極大的正數(shù)組成。

α=1.0×104～13×max(剛度矩陣K的對(duì)象元素)

(27)

[K+Kα]u=F+Fα

(28)

其中

(29)

(30)

其中額外的矩陣Kα是由下式(31)所定義的節(jié)點(diǎn)矩陣組成的全局懲罰矩陣:

(31)

向量Fα來(lái)自于初始邊界條件,它的節(jié)點(diǎn)矩陣采用下列形式:

(32)

1.3 背景網(wǎng)格積分方法

定義網(wǎng)格為mi×mi×mi=nN,nN即為節(jié)點(diǎn)的總數(shù),對(duì)于KIJ、FI轉(zhuǎn)化為如下公式進(jìn)行積分。

(33)

(34)

(35)

(36)

1.4 無(wú)單元Galerkin方法程序求解流程

2 無(wú)單元Galerkin方法解決三維算例

解決三維泊松方程的數(shù)值算例,通過(guò)上文中提到的無(wú)單元Galerkin方法具體解決一個(gè)三維算例。

(37)

邊界條件為

圖1 三維數(shù)值算例選取節(jié)點(diǎn)5×6×20的節(jié)點(diǎn)分布Fig.1 Node selection for three-dimensional numerical examples node distribution of 5×6×20

(38)

x∈[0,1],y∈[0,1],z∈[0,1]

(39)

解析解為:

u=sinπxsinπysinπz

(40)

例如對(duì)于節(jié)點(diǎn)5×6×20，節(jié)點(diǎn)分布如圖1所示。

通過(guò)觀察解析解和數(shù)值解的誤差,得到EFG方法求解三維泊松方程的計(jì)算精度。其次,采用不同的節(jié)點(diǎn)分布和恰當(dāng)?shù)膁max,觀察對(duì)于相對(duì)誤差值的影響。

表1 不同布點(diǎn)方式下EFG方法的 相對(duì)誤差值和所用時(shí)間(一)Table 1 The relative error of EFG method and the time used (一)

表2 不同布點(diǎn)方式下EFG方法的 相對(duì)誤差值和所用時(shí)間(二)Table 2 The relative error of EFG method and the time used (二)

表1和表2的第一行節(jié)點(diǎn)數(shù)5×5×5的數(shù)據(jù)為校準(zhǔn)數(shù)據(jù),其目的是保證上面兩個(gè)表的測(cè)試條件相同,以免因?yàn)閮?nèi)部因素或外部因素的差別而影響誤差。從表1的第二行看起,我們是固定x=5,y=8,依次增加z的取值,可以看出隨著z的增加,即節(jié)點(diǎn)數(shù)增多,計(jì)算時(shí)間也隨之增加,誤差越來(lái)越小,慢慢趨于穩(wěn)定。單獨(dú)看表2,從第二行看起同樣是固定x和y的值,單純?cè)黾觶,與上述結(jié)論一致。以表1為比較組,從兩個(gè)表第二行看起,表2較表1相當(dāng)于固定了x和z,在對(duì)應(yīng)的同一行中表2中的y等于表1中的y+1,顯然表2中的誤差更小,同樣證明節(jié)點(diǎn)數(shù)越多,誤差越小,擬合效果更佳,但需要的計(jì)算時(shí)間也更久。

表3 不同布點(diǎn)方式下EFG方法的 相對(duì)誤差值和所用時(shí)間(三)Table 3 The relative error of EFG method and the time used (三)

表4 不同的dmax下EFG方法的 相對(duì)誤差值和所用時(shí)間Table 4 Different dmaxthe relative error value of the EFG method and the time used

為了保證該方法的科學(xué)性和有效性,我們首先要做的就是實(shí)驗(yàn)的完整性,表3是對(duì)于固定y,z,增加x,觀察EFG方法所造成的相對(duì)誤差值和時(shí)間。顯然隨節(jié)點(diǎn)數(shù)增加誤差會(huì)減小,計(jì)算時(shí)間增加,與上述結(jié)論一致。綜上3個(gè)表我們可以看出,隨著節(jié)點(diǎn)數(shù)的增加,誤差會(huì)逐漸降低。

對(duì)于表1、表2和表3是在改變節(jié)點(diǎn)數(shù)并且固定了dmax=1.22情況下運(yùn)用EFG的方法,通過(guò)觀察相對(duì)誤差和運(yùn)算時(shí)間,首先我們可以看出隨著節(jié)點(diǎn)數(shù)的增加,計(jì)算的精確度也隨之增加,相對(duì)誤差值減小,但是同時(shí)也增加了計(jì)算量,所以所用的時(shí)間也增加了。當(dāng)節(jié)點(diǎn)數(shù)量增加時(shí)基函數(shù)的更高完備性階數(shù)實(shí)現(xiàn)比低階有更好的收斂性特征。相反地,從表4可以看出,在節(jié)點(diǎn)數(shù)不改變的情況,緩慢增加dmax的數(shù)值,相對(duì)誤差值并不是一直減小或增大,最終會(huì)圍繞著解析解擺動(dòng)；其中dmax=3.50,是所能測(cè)試到的最大值,超過(guò)它將無(wú)法產(chǎn)生運(yùn)行結(jié)果。

圖2 節(jié)點(diǎn)數(shù)x=8,y=20,dmax=1.60, z取不同值時(shí)的曲線擬合圖Fig.2 Node number x=8,y=20 dmax=1.60 curve fitting graphs with different value z

下面選擇上述測(cè)試結(jié)果的最優(yōu)條件,即在x=8,y=20,dmax=1.60,z取不同值的情況下的擬合點(diǎn)和解析解的圖像關(guān)系如圖2所示。

3 結(jié)論

從解決三維泊松方程可以看出,對(duì)于一個(gè)恰當(dāng)?shù)膁max,當(dāng)節(jié)點(diǎn)數(shù)目增加時(shí)擬合效果更好。而當(dāng)節(jié)點(diǎn)數(shù)保持不變,單獨(dú)增加dmax的值時(shí),相對(duì)誤差值并不是一直減小,最終會(huì)圍繞著解析解擺動(dòng)。

總之,增加節(jié)點(diǎn)數(shù)目對(duì)于減小相對(duì)誤差是有效果的,也就是增加節(jié)點(diǎn)數(shù)目,細(xì)化節(jié)點(diǎn)單元可達(dá)到更好的擬合效果。其次使用EFG方法解決具有特殊邊界的三維問(wèn)題,應(yīng)當(dāng)用一個(gè)恰當(dāng)?shù)膁max和較多的節(jié)點(diǎn)。