賈亞茹,付 偉,曹 健

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

無網(wǎng)格方法是近年來興起的一種數(shù)值計(jì)算方法。1994年,Belyschko等人[1]通過采用MLS法計(jì)算形函數(shù),并利用Lagrange乘子法引入本質(zhì)邊界條件,提出了無單元伽遼金法(EFG),改進(jìn)了之前Nayroles等[2]提出的彌散單元法(DEM),提高了算法質(zhì)量,為無網(wǎng)格方法的發(fā)展邁出了關(guān)鍵一步,掀起了無網(wǎng)格方法研究的熱潮。使用MLS方法構(gòu)造形函數(shù)時(shí),所得到的最后的代數(shù)方程組有可能是病態(tài)的,導(dǎo)致得不到理想的數(shù)值解。2013年,顧天奇、張雷等人[3]提出了改進(jìn)的MLS方法,解決了這一問題。2013年,張贊等人[4]提出了關(guān)于三維瞬態(tài)熱傳導(dǎo)問題的改進(jìn)無網(wǎng)格伽遼金法。2014年,王聚豐、程玉民等人[5]對插值型移動最小二乘法進(jìn)行了誤差估計(jì),得出結(jié)論,如果多項(xiàng)式基函數(shù)的階足夠大且原始函數(shù)足夠光滑,則逼近函數(shù)及其偏導(dǎo)數(shù)在影響域的最大半徑上收斂到精確值。2015年,孫鳳欣等人[6]在n維空間內(nèi)對插值型最小二乘法進(jìn)行了誤差估計(jì),在王聚豐、程玉民等人的基礎(chǔ)上提出了基于IMLS方法的插值無單元Galerkin(IEFG)方法。2018年,蔡小杰等人[7]建立了彈塑性大變形問題的改進(jìn)的無單元Galerkin方法,此方法是基于改進(jìn)的移動最小二乘法建立形函數(shù),根據(jù)彈塑性大變形問題的Galerkin弱形式建立離散方程,利用罰函數(shù)法施加位移邊界條件,推導(dǎo)了彈塑性大變形問題的改進(jìn)的無單元Galerkin方法的公式,采用Newton-Raphson迭代法進(jìn)行求解。

無網(wǎng)格法發(fā)展短短數(shù)年,雖為工程計(jì)算問題帶來一些新方法所具有的變革,但事實(shí)上是遠(yuǎn)遠(yuǎn)不夠的,其中還有一些不足和有待研究者來解決的問題。

本文介紹了無單元伽遼金方法中構(gòu)造形函數(shù)的方法——移動最小二乘法,重點(diǎn)探究了權(quán)函數(shù)的影響半徑對無網(wǎng)格方法精度的影響。結(jié)合權(quán)函數(shù)的性質(zhì)和形式,分析了權(quán)函數(shù)的影響半徑對于數(shù)值結(jié)果的影響,給出了結(jié)合參數(shù)構(gòu)造的新的權(quán)函數(shù)的影響半徑的確定方法和基于能量范數(shù)的誤差模型。通過對比不同參數(shù)值下的誤差的能量范數(shù),得出最優(yōu)的參數(shù)值,從而確定權(quán)函數(shù)的影響半徑的取值,并對實(shí)際的數(shù)值算例進(jìn)行了研究,通過設(shè)定不同積分網(wǎng)格分布下的Laplace方程和Poisson方程對數(shù)值結(jié)果進(jìn)行分析,結(jié)合數(shù)值算例證明了給出的新的權(quán)函數(shù)的影響半徑的確定方法的有效性。

1 移動最小二乘法(MLS)

無單元伽遼金方法中運(yùn)用移動最小二乘法確定形函數(shù)。移動最小二乘法的優(yōu)勢在于可以通過不同階數(shù)的基函數(shù)和不同的權(quán)函數(shù)提高移動最小二乘近似的計(jì)算精度[8]。

在移動最小二乘法(MLS)中,函數(shù)u(X)近似為:

(1)

這里m是基函數(shù)的個數(shù),PT(X)是基函數(shù),a(X)是相應(yīng)的系數(shù)。

基函數(shù)在二維空間的形式為:

線性基:

PT(X)=[1,x,y]

(2)

二次基:

PT(X)=[1,x,y,x2,xy,y2]

(3)

定義如下函數(shù):

(4)

對加權(quán)殘量泛函J關(guān)于a(x)求偏導(dǎo),

(5)

令式(5)中,

(6)

(7)

從而有

A(X)a(X)=B(X)Us,

(8)

其中,Us={u1,u2,…,un}T為積分域中所有節(jié)點(diǎn)形成的節(jié)點(diǎn)參數(shù)向量。

針對一個二維問題,選取線性基作為基函數(shù),

(9)

(10)

根據(jù)移動最小二乘法(MLS)可得到逼近函數(shù)uh(X)的表達(dá)式為:

(11)

其中，ΦT(X)定義為點(diǎn)x所在的積分單元內(nèi)所有節(jié)點(diǎn)的MLS形函數(shù),具體表示為:

(12)

在移動最小二乘法中,求解的關(guān)鍵之一在于需要通過合理控制支持域中節(jié)點(diǎn)的數(shù)量來保證A的可逆性,從而獲得MLS法的穩(wěn)定性。另外,待定系數(shù)a(X)是一個關(guān)于X的函數(shù),這就使得整個未知函數(shù)的擬合近似成為一個動態(tài)過程,同時(shí)在定義域內(nèi)也具有了連續(xù)性。以上結(jié)論都是建立在選取了合適的權(quán)函數(shù),由此說明,在MLS方法中權(quán)函數(shù)對于數(shù)值結(jié)果的質(zhì)量有著極為重要的影響。

1.1 權(quán)函數(shù)的性質(zhì)

權(quán)函數(shù)的選取應(yīng)根據(jù)以下幾方面的性質(zhì)來提升整個算法的動態(tài)性,支持MLS法的相容性[9]。

在這幾個性質(zhì)中,緊支性是最重要的,它對于局部區(qū)域數(shù)值的控制使得近似成為可能。

1.2 權(quán)函數(shù)的形式

權(quán)函數(shù)的形式多種多樣,目前常用的權(quán)函數(shù)有指數(shù)型、Gauss型、樣條函數(shù)等。在解決實(shí)際問題時(shí)應(yīng)注重權(quán)函數(shù)的改變?yōu)樗惴◣淼木扔绊憽?/p>

1)指數(shù)函數(shù)

(13)

2)Gauss函數(shù)

(14)

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

(15)

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

(16)

MLS近似法可以非常接近于內(nèi)插方法(局部近似)或類似于傳統(tǒng)最小二乘法近似(全局近似),這都要取決于所采用的權(quán)函數(shù)的影響半徑。要使得MLS近似更加精確,對每個權(quán)函數(shù)影響范圍的重疊區(qū)域要有足夠的要求,同時(shí)要保證影響域內(nèi)有足夠的節(jié)點(diǎn)數(shù)量。權(quán)函數(shù)的影響半徑選取要使得影響域內(nèi)節(jié)點(diǎn)數(shù)量不能出現(xiàn)太多或太少的情況,一方面要求影響域內(nèi)的節(jié)點(diǎn)數(shù)量應(yīng)足夠多,這是為了保證矩陣A(x)的可逆性,另一方面保證影響半徑足夠小,MLS方法是一種注重局部特性的研究方法,過大的影響域會使得擬合效果失真,計(jì)算效率降低,為此保證影響半徑足夠小也是提高擬合效果的重要方法。

2 誤差估計(jì)

本文使用相同的權(quán)函數(shù),在不同的影響半徑下使用無網(wǎng)格伽遼金法來計(jì)算得到的不同的近似結(jié)果。

定義是三個影響半徑rinf-min,rinf,rinf-max:

1)rinf:rinf=α×dI,其中dI為距離節(jié)點(diǎn)I第n個最近節(jié)點(diǎn)的距離,α為正實(shí)數(shù)。

2)rinf-min:rinf-min=prod×dI,其中prod∈∪(1.5,δ)。

3)rinf-max:rinf-max=β×rinf-min,其中β∈∪(2,δ)。

2.1 誤差公式的定義

誤差,定義為近似解與精確解之間的差異。無網(wǎng)格伽遼金法大多被應(yīng)用于解決工程力學(xué)方面的實(shí)際問題,其中會產(chǎn)生位移和應(yīng)變兩個值,而應(yīng)變的誤差研究比位移更重要,因?yàn)橥ㄟ^大量的數(shù)值計(jì)算發(fā)現(xiàn)位移所表現(xiàn)的誤差沒有應(yīng)變反應(yīng)的誤差更靈敏,所以這里采用能量范數(shù)誤差對每個積分單元進(jìn)行誤差估計(jì)。

本文的重點(diǎn)在于,通過誤差的變化來確定參數(shù)prod、β的取值,從而獲得最優(yōu)誤差下的影響半徑。

2.2 誤差范數(shù)

(17)

誤差的能量范數(shù)為:

(18)

對于每個單元,有

(19)

結(jié)構(gòu)的總體能量范數(shù)

(20)

2.3 誤差指標(biāo)

總體誤差指標(biāo)以η表示,則

(21)

對于每個單元,誤差指標(biāo)

(22)

其中,‖e‖a是單元應(yīng)力誤差能量范數(shù)‖e‖e的允許值。

3 數(shù)值算例

3.1 數(shù)值算例1——矩形域上的Laplace方程

選取Laplace方程作為算例,通過判斷誤差的大小對新的影響辦法選取方法進(jìn)行可行性檢驗(yàn)。

方程為:

(23)

邊界條件:

(24)

解析解:

u(x,y)=-x3-y3+3x2y+3xy2

(25)

在本算例中,選取四次樣條權(quán)函數(shù),懲罰因子為1015,矩形域[0,1]×[0,1]內(nèi)分別均勻分布了三種節(jié)點(diǎn):15×15,25×25,35×35,采用的背景網(wǎng)格與節(jié)點(diǎn)分布相同。在影響半徑rinf-min=prod×dI中,令dI=1。

為了確定prod和β對于誤差的影響,本文分別對這兩個參數(shù)進(jìn)行了研究:

3.1.1參數(shù)prod的確定

表1 prod取不同值,β=2.0時(shí)的能量范數(shù)誤差

圖1 prod取不同值,β=2.0時(shí),不同節(jié)點(diǎn)下的能量范數(shù)誤差

由圖1可以看出,在三種不同的網(wǎng)格劃分下,參數(shù)prod在1.2處都處于曲線的最低點(diǎn),能量范數(shù)誤差在prod=1.2時(shí)取得最優(yōu)值。同時(shí),也可以看出網(wǎng)格劃分越多,節(jié)點(diǎn)密度越大,能量范數(shù)誤差越小,數(shù)值的近似結(jié)果越精準(zhǔn)。

3.1.2參數(shù)β的確定

表2 β取不同值,prod=1.2時(shí)的能量范數(shù)誤差

圖2 β取不同值,prod=1.2時(shí),不同節(jié)點(diǎn)下的能量范數(shù)誤差

由圖2可以看出,能量范數(shù)誤差在β=2.0時(shí)達(dá)到最優(yōu),同時(shí)網(wǎng)格劃分密度越大,能量范數(shù)誤差越小,數(shù)值近似結(jié)果越精確。

3.2 數(shù)值算例2——矩形域上的Poisson方程

方程為:

(26)

邊界條件為:

(27)

解析解為:

u(x,y)=x2+y2

(28)

在本算例中,選取四次樣條權(quán)函數(shù),懲罰因子為1015,矩形域[0,1]×[0,1]內(nèi)分別均勻分布了三種節(jié)點(diǎn):15×15,25×25,35×35,采用的背景網(wǎng)格與節(jié)點(diǎn)分布相同。在影響半徑rinf-min=prod×dI中,令dI=1。

3.2.1參數(shù)prod的確定

表3 prod取不同值,β=2.2時(shí)的能量范數(shù)誤差

圖3 prod取不同值,β=2.2時(shí),不同節(jié)點(diǎn)下的能量范數(shù)誤差

由圖3可以看出,在三種不同的網(wǎng)格劃分下,參數(shù)prod在1.4處都處于曲線的最低點(diǎn),能量范數(shù)誤差在prod=1.4時(shí)取得最優(yōu)值。同時(shí),也可以看出網(wǎng)格劃分越多,節(jié)點(diǎn)密度越大,能量范數(shù)誤差越小,數(shù)值的近似結(jié)果越精準(zhǔn)。

3.2.2參數(shù)β的確定

表4 β取不同值,prod=1.4時(shí)的能量范數(shù)誤差

圖4 β取不同值,prod=1.4時(shí),不同節(jié)點(diǎn)下的能量范數(shù)誤差

由圖4可以看出,在三種不同的網(wǎng)格劃分下,參數(shù)β在2.2處都處于曲線的最低點(diǎn),能量范數(shù)誤差在β=2.2時(shí)取得最優(yōu)值。同時(shí),也可以看出網(wǎng)格劃分越多,節(jié)點(diǎn)密度越大,能量范數(shù)誤差越小,數(shù)值的近似結(jié)果越精準(zhǔn)。

4 結(jié)論

此誤差估計(jì)方法是以每個積分單元為單位進(jìn)行的,可以提高準(zhǔn)確性和有效性。

通過以上兩個數(shù)值算例的結(jié)果所示,當(dāng)β的取值在[2.0,2.2]之間,prod的取值在[1.2,1.4]之間時(shí),數(shù)值算例中在不同節(jié)點(diǎn)劃分下的能量范數(shù)誤差達(dá)到最小值。同時(shí),通過數(shù)值算例的分析研究可以看出,合理調(diào)整節(jié)點(diǎn)和幾分單元的劃分可以增強(qiáng)數(shù)值結(jié)果的精度?？偟膩碚f,適當(dāng)增加節(jié)點(diǎn)的密度劃分獲得的誤差結(jié)果越小,近似結(jié)果的準(zhǔn)確度越高。

選取合適的權(quán)函數(shù)的影響半徑對EFG方法的精度有重要影響,同時(shí)對誤差的估計(jì)有助于調(diào)整節(jié)點(diǎn)和積分單元的劃分,從而增強(qiáng)結(jié)果的準(zhǔn)確性。本文提出的確定權(quán)函數(shù)影響半徑的方法,通過了實(shí)際的數(shù)值算例驗(yàn)證了此方法的合理性,證明了使用此方法可以在問題域中得到較為理想的結(jié)果。