付 偉，解紅霞

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

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

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

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

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

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

1 無單元Galerkin方法

1.1 等效積分弱形式

這里引入一個泊松方程,來探討無單元Galerkin方法。

Poisson方程:

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

(1)

邊界條件:

(2)

(3)

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

(4)

其中,

(5)

已知

(6)

其中n是x點區(qū)域內(nèi)的節(jié)點個數(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)

然后我們對(1.12)的各項進行逐一積分,首先對第一項進行積分:

(13)

其中，

(14)

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

(15)

下面是對(12)第二項進行積分:

(16)

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

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

(17)

(18)

最后對(12)的第三項進行積分:

(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),這里α是罰因子的對角矩陣,針對于二維問題i=2,針對于三維問題時i=3,這個罰因子αk(k=1,2,…i),可能是一個坐標(biāo)函數(shù),通常情況它是由極大的正數(shù)組成。

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

(27)

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

(28)

其中

(29)

(30)

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

(31)

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

(32)

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

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

(33)

(34)

(35)

(36)

1.4 無單元Galerkin方法程序求解流程

2 無單元Galerkin方法解決三維算例

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

(37)

邊界條件為

圖1 三維數(shù)值算例選取節(jié)點5×6×20的節(jié)點分布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)

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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