劉衛(wèi)桃，孟智娟，李興國，任紅萍

(太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院，太原 030024)

無網(wǎng)格方法是一種只需要節(jié)點(diǎn)信息，不需要?jiǎng)澐志W(wǎng)格的方法。在處理不連續(xù)、大變形等數(shù)值問題時(shí)很有優(yōu)勢，因此在計(jì)算科學(xué)中的研究已很熱門[1-3]。1992年，Nayroles等將Galerkin方法與移動(dòng)最小二乘法(MLS)結(jié)合，提出了擴(kuò)散單元法(DEM)[4].但DEM方法邊界條件不易引入，對(duì)此Belystschko等人將DEM法進(jìn)行了改進(jìn)，提出了無單元Galerkin方法(EFG)[5]，是無網(wǎng)格方法中最具有發(fā)展前景的方法之一，并且利用該方法處理了動(dòng)態(tài)裂紋擴(kuò)展和大變形等復(fù)雜的三維問題。Zhang等人提出了改進(jìn)的無單元Galerkin方法(IEFG)，并求解了熱傳導(dǎo)等問題[6]。

為提高無網(wǎng)格方法的計(jì)算速度，Ren等人提出了插值型EFG方法，且應(yīng)用于求解彈性力學(xué)問題、三維彈塑性等問題[7-8]。張國達(dá)等人[9]研究了瞬態(tài)熱傳導(dǎo)問題的插值型EFG方法及誤差。針對(duì)求解三維問題時(shí)計(jì)算量大、計(jì)算效率低的問題，孟智娟提出的維數(shù)分裂無單元Galerkin方法(DSEFG)[10]解決了這一問題，并將其應(yīng)用于求解三維勢問題等。

由于DSEFG方法中構(gòu)造形函數(shù)[11]時(shí)權(quán)函數(shù)[12]的選擇會(huì)直接影響計(jì)算的精度和效率，因此研究權(quán)函數(shù)的選取具有關(guān)鍵意義。

1 三維勢問題的分維處理

三維勢問題的控制方程為：

(1)

本質(zhì)邊界條件為：

(2)

自然邊界條件為：

(3)

取x3為分裂方向，將問題域Ω在x3方向上分裂為L層。便得到L+1個(gè)二維子域Ω(k)，(k=0，

1，…，L)，則有：

(4)

(5)

相應(yīng)的邊界條件為：

(6)

(7)

對(duì)于求解問題(1)-(3)的數(shù)值解，選用的DSEFG方法是用IEFG方法求解二維邊值問題，在分裂方向x3上采用有限差分法。

采用罰函數(shù)法施加于Γu，則有：

(8)

(9)

2 三維勢問題的維數(shù)分裂無單元Galerkin方法(DSEFG)

(10)

2.1 權(quán)函數(shù)類型

(1)三次樣條函數(shù)

(11)

(2)四次樣條函數(shù)

(12)

(3)正定緊支徑向基函數(shù)(CSBRF2)

(13)

(4)指數(shù)型函數(shù)

(14)

2.2 三維勢問題的維數(shù)分裂無單元Galerkin方法

由方程(10)，可得:

Φ*(x(k))u″

(15)

(16)

這里

(17)

B(x(k))=

(B1(x(k))，B2(x(k))，…，Bn(x(k)))

(18)

(19)

將式(10)、式(15)-式(16)代入式(8)可得:

(20)

方程(20)的各項(xiàng)積分可討論為：

δuT·C·u″

(21)

δuT·K·u

(22)

(23)

(24)

δuT·Kα·u

(25)

(26)

其中:

(27)

(28)

(29)

(30)

(31)

(32)

將式(21)-式(26)代入方程式(20)得:

δuT·(Cu″-Ku+Kαu-F1-F2-Fα)=0

(33)

由δuT的任意性，可得:

(34)

其中

(35)

(36)

?

(37)

在x3上采用有限差分法，得:

(k=1，2，…，L-1)

(38)

從而，方程組(34)表示為：

(39)

?

(40)

其矩陣形式為：

(41)

這里

(42)

令

(43)

(44)

(45)

這樣，方程(41)簡化為：

EU=R

(46)

(47)

3 數(shù)值算例

在算例中，通過對(duì)三維勢問題計(jì)算結(jié)果收斂性的研究，說明了權(quán)函數(shù)在維數(shù)分裂無單元Galerkin方法中的重要性。算例中的每個(gè)積分單元上，選擇4×4高斯積分。

相對(duì)誤差定義為：

(48)

3.1 算例1

三維立方體內(nèi)具有Dirichlet邊界條件的Laplace方程

0(x∈Ω)

(49)

邊界條件為：

u=sin(πx2)sin(πx3)(x1=0)

(50)

u=2sin(πx2)sin(πx3)(x1=1)

(51)

u=0(x2=0，x2=1，x3=0，x3=1)

(52)

其中求解域Ω=[0，1]×[0，1]×[0，1].

對(duì)應(yīng)的解析解為：

(53)

將問題域Ω沿x1方向均勻地分裂為L個(gè)子域，在平面Ox2x3上節(jié)點(diǎn)均勻分布為10×10，即整個(gè)問題域Ω上的節(jié)點(diǎn)分布為(L+1)×10×10.

表1為指數(shù)型權(quán)函數(shù)在α不同時(shí)，分別對(duì)應(yīng)于最佳節(jié)點(diǎn)分布處，DSEFG方法取得的相對(duì)誤差最低和CPU運(yùn)行時(shí)間的結(jié)果。圖1顯示了L=10，α=0.3時(shí)，指數(shù)型權(quán)函數(shù)的DSEFG方法取得的相對(duì)誤差最低。

表1 以指數(shù)型權(quán)函數(shù)的DSEFG在不同參數(shù)下與相對(duì)應(yīng)的節(jié)點(diǎn)分布取到的最低相對(duì)誤差的比較Tab.1 Comparison of the lowest relative error between DSEFG with exponential weight function and corresponding node distribution under different parameters

圖1 指數(shù)型權(quán)函數(shù)下的DSEFG方法在α不同時(shí)的相對(duì)誤差Fig.1 Relative error of DSEFG method under exponential weight function with different α

圖2-圖5是DSEFG方法在L不同時(shí)的相對(duì)誤差。當(dāng)問題域Ω分裂的子域L不同時(shí)，隨著L的增加，以三次樣條，四次樣條為權(quán)函數(shù)時(shí)計(jì)算值逐漸收斂，以CSBRF2函數(shù)和指數(shù)型函數(shù)為權(quán)函數(shù)時(shí)計(jì)算值均是在一點(diǎn)處取到最低，且在以指數(shù)型函數(shù)為權(quán)函數(shù)時(shí)，DSEFG方法的計(jì)算精度最高。

圖2 三次樣條函數(shù)下的相對(duì)誤差與L的關(guān)系Fig.2 The relationship between relative error and L under cubic spline function

圖3 四次樣條函數(shù)下的相對(duì)誤差與L的關(guān)系Fig.3 The relationship between relative error and L under quartic spline weight function

圖4 CSBRF2權(quán)函數(shù)下的相對(duì)誤差與L的關(guān)系Fig.4 The relationship between relative error and L under CSBRF2 weight function

圖5 指數(shù)型權(quán)函數(shù)下的相對(duì)誤差與L的關(guān)系Fig.5 The relationship between relative error and L under exponential weight function

表2顯示了采用不同權(quán)函數(shù)得到的最低相對(duì)誤差和運(yùn)行時(shí)間，且通過比較得出在以指數(shù)型函數(shù)為權(quán)函數(shù)時(shí)，維數(shù)分裂無單元Galerkin方法的計(jì)算精度和計(jì)算時(shí)間都相對(duì)較低。

表2 不同權(quán)函數(shù)下DSEFG方法的相對(duì)誤差與CPU時(shí)間Tab.2 Relative error and CPU time of DSEFG method under different weight functions

圖6為以不同函數(shù)為權(quán)函數(shù)時(shí)，DSEFG方法在(x1，7/9，7/9)處的值，可以看出都能很好的接近于解析解。并且在以指數(shù)型函數(shù)為權(quán)函數(shù)時(shí)，DSEFG方法的運(yùn)行速度相對(duì)快一點(diǎn)。

圖6 不同權(quán)函數(shù)下的DSEFG方法在u(x1，7/9，7/9)處的計(jì)算結(jié)果Fig.6 Calculation results of DSEFG method at u(x1，7/9，7/9) under different weight functions

3.2 算例2

三維立方體內(nèi)具有Neumann邊界條件的Laplace方程

?2u(x)=0(x∈Ω)

(54)

邊界條件為：

(55)

(56)

(57)

(58)

其中求解域Ω=[0，1]×[0，1]×[0，1].

對(duì)應(yīng)的解析解

(59)

對(duì)于該問題，選取x3為分裂方向，節(jié)點(diǎn)分布為11×11×11，dmax=1.19.圖7顯示了在沿x3方向不同權(quán)函數(shù)下DSEFG方法的取值結(jié)果，都能較好的落在解析解上。表3顯示了在以指數(shù)型函數(shù)為權(quán)函數(shù)時(shí)，維數(shù)分裂無單元Galerkin方法具有較好的精度，且計(jì)算速度較快。

表3 不同權(quán)函數(shù)下DSEFG方法的相對(duì)誤差與時(shí)間Tab.3 Relative error and time of DSEFG method under different weight functions

圖7 不同權(quán)函數(shù)下DSEFG方法的比較結(jié)果Fig.7 Comparison results of DSEFG method under different weight functions

4 結(jié)論

本文中，通過數(shù)值算例說明了利用DSEFG方法求解三維勢問題時(shí)，權(quán)函數(shù)及權(quán)函數(shù)中的參數(shù)對(duì)數(shù)值計(jì)算結(jié)果的精度和速率都有重要的影響。并且得出，求解三維Laplace方程時(shí)，權(quán)函數(shù)選擇指數(shù)型函數(shù)，維數(shù)分裂無單元Galerkin方法可以得到較好的結(jié)果。