無單元Galerkin方法數值求解三維泊松方程

付 偉，解紅霞

(1.太原科技大學 應用科學學院,山西 太原 030024；2.太原學院 應用數學系，山西 太原 030032)

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

無網格FEM法[3]是Idelsohn等人于2003年提出的,融合了無網格法與FEM法的優點的一類新的數值方法。

2005年,程玉民等[4]提出了復變量MLS法。與較傳統的MLS法相比,該方法減少了函數中的變量,有效實現了降維,節約了計算時間,并且有更高的計算精度。

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

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

1 無單元Galerkin方法

1.1 等效積分弱形式

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

Poisson方程:

u(x)-b(x)=0(在Ω內)

(1)

邊界條件:

(2)

(3)

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

(4)

其中,

(5)

已知

(6)

其中n是x點區域內的節點個數,Φ(x)函數，

Φ(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)是給定源函數所引起的荷載列陣

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 邊界條件的處理

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

(26)

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

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

(27)

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

(28)

其中

(29)

(30)

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

(31)

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

(32)

1.3 背景網格積分方法

定義網格為mi×mi×mi=nN,nN即為節點的總數,對于KIJ、FI轉化為如下公式進行積分。

(33)

(34)

(35)

(36)

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

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

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

(37)

邊界條件為

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

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

通過觀察解析解和數值解的誤差,得到EFG方法求解三維泊松方程的計算精度。其次,采用不同的節點分布和恰當的dmax,觀察對于相對誤差值的影響。

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

表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

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

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

圖2 節點數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

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

3 結論

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

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