節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法1)

樊禮恒 王東東 劉宇翔 杜洪輝

(廈門(mén)大學(xué)土木工程系,福建廈門(mén) 361005)

(廈門(mén)市交通基礎(chǔ)設(shè)施智能管養(yǎng)工程技術(shù)研究中心,福建廈門(mén) 361005)

引言

基于變分法或等效積分弱形式的伽遼金法和基于強(qiáng)形式的配點(diǎn)法是計(jì)算力學(xué)領(lǐng)域兩種常用的數(shù)值計(jì)算方法[1].與伽遼金法相比,配點(diǎn)法具有構(gòu)造簡(jiǎn)單、計(jì)算高效的特點(diǎn),但是配點(diǎn)法需要計(jì)算數(shù)值離散形函數(shù)的高階梯度,例如勢(shì)問(wèn)題或彈性力學(xué)問(wèn)題等二階問(wèn)題需要用到形函數(shù)的二階導(dǎo)數(shù).同時(shí)配點(diǎn)法的剛度矩陣不對(duì)稱(chēng),穩(wěn)定性相對(duì)較差.有限元法是目前應(yīng)用最為廣泛的一種典型伽遼金法[1],其形函數(shù)定義在單元上,一般網(wǎng)格情況下僅具有C0連續(xù)性[2],在相鄰單元邊界處會(huì)出現(xiàn)梯度跳躍,因而不能直接應(yīng)用于需要計(jì)算形函數(shù)高階梯度的配點(diǎn)法[1].近年來(lái),隨著采用高階光滑形函數(shù)的無(wú)網(wǎng)格法和等幾何分析的快速發(fā)展[3-7],配點(diǎn)法逐漸得到了愈來(lái)愈廣泛的應(yīng)用,例如徑向基函數(shù)配點(diǎn)法[8-13],移動(dòng)最小二乘與再生核無(wú)網(wǎng)格配點(diǎn)法[14-18],等幾何配點(diǎn)法[19-21]等.值得指出的是,與無(wú)網(wǎng)格和等幾何形函數(shù)相比[1],有限元形函數(shù)形式簡(jiǎn)單,計(jì)算快捷.最近,高效偉等[22-23]利用有限元形函數(shù)在單元內(nèi)部高階連續(xù)的特點(diǎn),通過(guò)在選定配點(diǎn)鄰域內(nèi)構(gòu)建局部單元,提出了自由單元配點(diǎn)法.

就配點(diǎn)法的精度而言,一個(gè)突出的問(wèn)題就是其計(jì)算精度依賴(lài)于形函數(shù)所采用的基函數(shù)階次的奇偶性,即奇數(shù)階次基函數(shù)的配點(diǎn)法收斂率較其基函數(shù)階次下降一次[19,24].王東東等[24-26]系統(tǒng)研究了二階和四階問(wèn)題無(wú)網(wǎng)格配點(diǎn)法奇數(shù)階次基函數(shù)精度掉階的原因.研究表明,配點(diǎn)法的精度與形函數(shù)及其梯度的一致性或完備性條件密切相關(guān).因此,王東東等[24-25]通過(guò)引入無(wú)網(wǎng)格形函數(shù)的梯度光滑方法,建立了超收斂無(wú)網(wǎng)格配點(diǎn)法,并在不引入待定配點(diǎn)位置的情況下有效解決了奇數(shù)階次基函數(shù)配點(diǎn)法的精度掉階問(wèn)題,實(shí)現(xiàn)了采用奇數(shù)次基函數(shù)配點(diǎn)法精度的超收斂提升.基于梯度光滑方法,鄧立克等[27]通過(guò)線(xiàn)性基函數(shù)構(gòu)造了二階光滑梯度,使得采用線(xiàn)性基函數(shù)就可以進(jìn)行四階薄板問(wèn)題的伽遼金無(wú)網(wǎng)格分析.然而,注意到上述工作中所討論的形函數(shù)高階光滑梯度構(gòu)造的基礎(chǔ)是形函數(shù)的標(biāo)準(zhǔn)一階梯度,這對(duì)于本身具有高階光滑特性的無(wú)網(wǎng)格形函數(shù)沒(méi)有問(wèn)題,但是由于有限元形函數(shù)的一階梯度在單元節(jié)點(diǎn)和邊界處并不存在,因此難以直接構(gòu)造有限元形函數(shù)的高階光滑梯度和發(fā)展相應(yīng)的配點(diǎn)法.

為了構(gòu)造有限元形函數(shù)的高階光滑梯度,本文在無(wú)網(wǎng)格法的梯度光滑理論框架內(nèi)[24,28],首先引入了有限元形函數(shù)在節(jié)點(diǎn)處和全域內(nèi)的一階光滑梯度定義,然后在此基礎(chǔ)上建立了有限元形函數(shù)的二階光滑梯度,進(jìn)而采用有限元形函數(shù)的光滑梯度建立了一種新型的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法.與傳統(tǒng)有限元形函數(shù)梯度不同,有限元形函數(shù)的一階和二階光滑梯度具有全域連續(xù)的特點(diǎn).不失一般性,本文以線(xiàn)性形函數(shù)為例,詳細(xì)研究了有限元形函數(shù)一階和二階光滑梯度的構(gòu)造特點(diǎn)和一致性條件,并對(duì)節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法進(jìn)行了精度分析.分析表明,線(xiàn)性有限元形函數(shù)的二階光滑梯度能夠滿(mǎn)足二階一致性條件,因而對(duì)應(yīng)的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法具有二次精度和超收斂特性.注意到基于梯度光滑理論[28],Liu 等[29]構(gòu)造了有限元形函數(shù)的一階光滑梯度,提出了伽遼金光滑有限元法.但本文所討論的一階與二階有限元形函數(shù)光滑梯度,其光滑核函數(shù)選取和構(gòu)造方式均與伽遼金光滑有限元法不同,同時(shí)著重研究的是配點(diǎn)法計(jì)算列式.文中通過(guò)典型算例系統(tǒng)驗(yàn)證了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的精度和收斂性.

1 配點(diǎn)法控制方程

本文考慮穩(wěn)態(tài)勢(shì)問(wèn)題和彈性力學(xué)問(wèn)題等典型二階問(wèn)題,其控制方程的強(qiáng)形式可以表示為

其中u(x)為場(chǎng)變量,例如彈性力學(xué)問(wèn)題的位移或熱傳導(dǎo)問(wèn)題的溫度,L 為空間域? 內(nèi)的偏微分算子,Bg和Bt為強(qiáng)制邊界Γg和自然邊界Γt上的算子,b為源項(xiàng),g和t分別為Γg和Γt上給定的邊界值.方便起見(jiàn),這里給出二維情況下的各個(gè)算子的定義.例如,對(duì)于穩(wěn)態(tài)勢(shì)問(wèn)題,L,Bg和Bt為

其中n表示在自然邊界Γt的單位外法線(xiàn)向量.對(duì)于平面應(yīng)力彈性力學(xué)問(wèn)題,上述算子可表示為

式(4)和式(6)中D為平面應(yīng)力彈性本構(gòu)矩陣,I為二階單位陣,E和ν 分別代表材料的楊氏模量和泊松比.

將式(7)代入控制方程(1),并令其離散節(jié)點(diǎn)處的殘值為零,便可以得到配點(diǎn)法的離散方程

其中K,d和f分別表示剛度矩陣,節(jié)點(diǎn)系數(shù)向量和力向量

由式(10)可見(jiàn),配點(diǎn)法需要計(jì)算形函數(shù)的二階梯度,而標(biāo)準(zhǔn)有限元形函數(shù)在單元邊界處不存在二階梯度,因此無(wú)法直接用于配點(diǎn)法.此外,雖然通過(guò)特殊構(gòu)造方式可以建立一些特殊單元的高階光滑有限元形函數(shù),但是目前仍然缺乏構(gòu)造一般單元高階光滑形函數(shù)的通用有效方法[1].本文則是通過(guò)引入無(wú)網(wǎng)格法中的梯度光滑技術(shù),建立了有限元形函數(shù)的一階和二階光滑梯度,進(jìn)而發(fā)展了相應(yīng)的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法.

2 有限元形函數(shù)的光滑梯度構(gòu)造

2.1 有限元形函數(shù)光滑梯度構(gòu)造理論

考慮圖1 所示的有限元網(wǎng)格,?L表示第L個(gè)單元對(duì)應(yīng)的空間區(qū)域,節(jié)點(diǎn)xI對(duì)應(yīng)的有限元形函數(shù)是NI(x),NI(x)的影響域?yàn)?/p>

式中MI表示共享節(jié)點(diǎn)xI的單元個(gè)數(shù).方便起見(jiàn),定義下面的節(jié)點(diǎn)組和單元組

根據(jù)廣義梯度光滑理論[4,24],有限元形函數(shù)NI(x)的一階光滑梯度的定義為

圖1 有限元形函數(shù)影響域示意圖Fig.1 Illustration of the influence domains of finite element shape functions

式中φ(x,y)為梯度光滑核函數(shù),NI,i(x)為傳統(tǒng)的標(biāo)準(zhǔn)一階梯度或?qū)?shù),即NI,i(x)=?NI(x)/?xi,xi為x的第i個(gè)坐標(biāo)分量.由于有限元形函數(shù)的梯度在單元邊界處沒(méi)有定義,因此式(14)并不能直接用于有限元形函數(shù)的光滑梯度構(gòu)造.

為了構(gòu)造有限元形函數(shù)的光滑梯度,這里首先引入有限元形函數(shù)在節(jié)點(diǎn)處的梯度值.為此,在式(14)中選取如下的梯度光滑核函數(shù)

若進(jìn)一步考慮常應(yīng)變單元,例如一維兩節(jié)點(diǎn)單元、二維三節(jié)點(diǎn)三角形單元和三維四節(jié)點(diǎn)四面體單元,則NI,i(x)在每個(gè)單元內(nèi)均為常數(shù),式(16)可以簡(jiǎn)化為

其中VL表示單元?L的長(zhǎng)度、面積或體積.

在式(17)定義的有限元形函數(shù)的節(jié)點(diǎn)梯度基礎(chǔ)上,進(jìn)一步選擇梯度光滑核函數(shù)φ(x,xJ)=NJ(x),則根據(jù)式(14)對(duì)應(yīng)的離散形式,可得如下的有限元形函數(shù)一階光滑梯度

接下來(lái),對(duì)式(18)進(jìn)行求導(dǎo),并用式(18)定義的一階光滑梯度代替?zhèn)鹘y(tǒng)一階梯度,最后可得有限元形函數(shù)的二階光滑梯度

至此,已經(jīng)構(gòu)建有限元形函數(shù)的二階光滑梯度.下面通過(guò)一維和二維線(xiàn)性有限元形函數(shù)一階和二階光滑梯度的構(gòu)造過(guò)程來(lái)說(shuō)明有限元形函數(shù)光滑梯度構(gòu)造理論.

2.2 線(xiàn)性有限元形函數(shù)光滑梯度構(gòu)造

清晰起見(jiàn),首先考慮圖2 所示的一維均勻有限元網(wǎng)格,其中h表示單元的長(zhǎng)度,在此情況下,?I=,因此VI=h,=2h.有限元節(jié)點(diǎn)xI的形函數(shù)NI(x)可表示為

其對(duì)應(yīng)的一階梯度或?qū)?shù)NI,x(x)為

由式(20)、式(21)和圖2 可見(jiàn),有限元形函數(shù)NI(x)為C0連續(xù),其一階梯度NI,x(x)在單元節(jié)點(diǎn)xI?1,xI及xI+1 處不連續(xù),無(wú)法直接計(jì)算二階梯度,因而無(wú)法用于式(8)定義的配點(diǎn)法求解.

另一方面,根據(jù)式(14),有限元形函數(shù)在節(jié)點(diǎn)處的一階光滑梯度為

圖2 一維有限元形函數(shù)一階光滑梯度構(gòu)造過(guò)程Fig.2 Construction of the first order smoothed gradient of 1D finite element shape function

在此基礎(chǔ)上,利用式(18),可得一維線(xiàn)性有限元形函數(shù)的一階光滑梯度

進(jìn)一步將式(22)和式(23)代入式(19),可得一維線(xiàn)性有限元形函數(shù)的二階光滑梯度

圖2 和圖3 給出了線(xiàn)性有限元形函數(shù)一階和二階光滑梯度的構(gòu)造過(guò)程.由圖2 和圖3 可見(jiàn),雖然線(xiàn)性有限元形函數(shù)的梯度在節(jié)點(diǎn)處不連續(xù),而且二階梯度不存在,但由本文所提有限元形函數(shù)梯度光滑理論構(gòu)造所得的一階和二階光滑梯度仍然是連續(xù)函數(shù),可以直接用于配點(diǎn)法求解.

對(duì)于二維和三維情況,考慮三節(jié)點(diǎn)三角形單元和四節(jié)點(diǎn)四面體單元,其對(duì)應(yīng)的形函數(shù)為標(biāo)準(zhǔn)的線(xiàn)性形函數(shù)[1],這里不再贅述.二維和三維有限元形函數(shù)的一階和二階光滑梯度也可按照上述步驟進(jìn)行構(gòu)造.圖4 給出了基于平面三角形線(xiàn)性單元構(gòu)造一階光滑梯度和和二階光滑梯度的過(guò)程.與一維情況相同,線(xiàn)性三角形單元的形函數(shù)在節(jié)點(diǎn)和邊界處的一階梯度均不連續(xù),不存在二階梯度.另一方面,基于本文的有限元形函數(shù)光滑梯度構(gòu)造理論,可以方便地構(gòu)造圖5 所示的一階和二階光滑梯度.此外,圖2～圖5 同時(shí)表明,與有限元形函數(shù)和其標(biāo)準(zhǔn)梯度不同,有限元形函數(shù)的一階和二階光滑梯度與形函數(shù)本身具有相同的連續(xù)性,但光滑梯度的影響域則大于形函數(shù)的影響域.

圖3 一維有限元形函數(shù)二階光滑梯度構(gòu)造過(guò)程Fig.3 Construction of the second order smoothed gradient of 1D finite element shape function

圖4 二維有限元形函數(shù)一階光滑梯度構(gòu)造過(guò)程Fig.4 Construction of the first order smoothed gradient of 2D finite element shape function

圖5 二維有限元形函數(shù)二階光滑梯度構(gòu)造過(guò)程Fig.5 Construction of the second order smoothed gradient of 2D finite element shape function

3 節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法及其精度分析

將式(18)和式(19)定義的有限元一階和二階節(jié)點(diǎn)光滑梯度代入配點(diǎn)法列式(8),便形成了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法.以勢(shì)問(wèn)題為例,節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法離散方程剛度矩陣的各個(gè)元素為

節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法力向量的各個(gè)元素與式(10)相同.為不失一般性,下面首先分析線(xiàn)性有限元形函數(shù)光滑梯度的一致性條件,然后以勢(shì)問(wèn)題為例對(duì)基于線(xiàn)性單元的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的計(jì)算精度進(jìn)行理論研究.

3.1 線(xiàn)性有限元形函數(shù)光滑梯度的一致性條件

對(duì)于傳統(tǒng)線(xiàn)性有限元,其形函數(shù)及其一階梯度的一致性或完備性條件為[1]

其中xI j為節(jié)點(diǎn)xI的第j個(gè)坐標(biāo)分量,δi j為克羅內(nèi)克符號(hào).注意到有限元形函數(shù)的一致性條件取決于形函數(shù)的階次.例如,線(xiàn)性有限元形函數(shù)僅滿(mǎn)足式(26)～式(29)定義的線(xiàn)性一致性條件,由于其二階梯度不存在,故不滿(mǎn)足更高一次的二階一致性條件.但是對(duì)于有限元形函數(shù)的光滑梯度,接下來(lái)證明其不僅滿(mǎn)足線(xiàn)性一致性條件,在均布網(wǎng)格情況下還滿(mǎn)足二階一致性條件.

對(duì)于均勻有限元網(wǎng)格,由形函數(shù)的周期重復(fù)性可知,節(jié)點(diǎn)光滑梯度滿(mǎn)足如下條件[24-25]

其中k=0,1 或2.首先考慮k=2 的情況,利用式(32),式(33)變?yōu)?/p>

再者,對(duì)于k=1,式(33)變?yōu)?/p>

式(34)～式(36),證明了線(xiàn)性有限元形函數(shù)的二階光滑梯度滿(mǎn)足二階梯度一致性條件.與之形成鮮明對(duì)比的是,標(biāo)準(zhǔn)線(xiàn)性有限元形函數(shù)不存在二階梯度.

3.2 節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法精度分析

為不失一般性,以平面勢(shì)問(wèn)題為例,在均勻網(wǎng)格條件下對(duì)節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的精度進(jìn)行解析研究.節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的精度可以方便地采用如下的截?cái)嗾`差進(jìn)行度量[24,30]

其中uI為節(jié)點(diǎn)xI的解析常變量,即uI=u(xI),h表示單元的特征長(zhǎng)度,k為精度的階次.

在式(37)中引入uI的泰勒展開(kāi)形式

同時(shí)注意到在均布網(wǎng)格條件下,有限元形函數(shù)的光滑二階導(dǎo)數(shù)具有周期性和偶對(duì)稱(chēng)的特點(diǎn),因此當(dāng)n為奇數(shù)時(shí),有

將式(31),式(33)～式(36)及式(40)代入式(39),可得

由式(41)可以看出,節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法具有二階精度,其L2誤差和H1或能量誤差均為二階收斂.與線(xiàn)性有限元法相比,雖然兩者L2誤差收斂特性相同,但節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的H1或能量誤差具有高一階次的超收斂特性.

4 算例

本節(jié)通過(guò)一維、二維及三維算例驗(yàn)證節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的收斂性和精度,其中FEM 和FEC 分別表示傳統(tǒng)的伽遼金有限元法和本文所提的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法.為了和傳統(tǒng)有限元法進(jìn)行比較,算例中采用如下誤差形式

其中nd表示算例的空間維度,σij和εi j分別為應(yīng)力和應(yīng)變分量.

4.1 一維彈性桿問(wèn)題

考慮兩端固定的一維彈性桿問(wèn)題,其幾何和材料參數(shù)為:桿長(zhǎng)為L(zhǎng)=1,抗拉剛度EA=1,體力b(x)=?24x2.該問(wèn)題的解析解為u(x)=2x4.計(jì)算中采用圖6 和圖7 所示包含20,40 和80 個(gè)單元的均布和非均布有限元網(wǎng)格進(jìn)行分析.

圖6 一維彈性桿問(wèn)題均布有限元離散模型Fig.6 Uniform finite element meshes for the 1D elastic rod problem

圖7 一維彈性桿問(wèn)題非均布有限元離散模型Fig.7 Non-uniform finite element meshes for the 1D elastic rod problem

圖8 和圖9 給出了一維彈性桿問(wèn)題的L2和能量誤差收斂情況.結(jié)果顯示,節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的位移誤差收斂率與有限元法相同,精度略低,但是其能量誤差收斂率比有限元法高一階,對(duì)應(yīng)的精度也明顯優(yōu)于有限元法,很好地驗(yàn)證了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的精度和收斂率.圖10 為采用20 個(gè)均布和非均布單元得到的應(yīng)力結(jié)果,可見(jiàn)節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法(FEC)的應(yīng)力精度顯著高于傳統(tǒng)有限元法(FEM).

4.2 二維矩形區(qū)域勢(shì)問(wèn)題

考慮如圖11 所示的正方形勢(shì)問(wèn)題區(qū)域,其邊長(zhǎng)L=1,對(duì)應(yīng)勢(shì)問(wèn)題的解析解為

圖8 采用均布網(wǎng)格的一維彈性桿問(wèn)題收斂率對(duì)比Fig.8 Convergence comparison for the 1D elastic rod problem using uniform meshes

圖9 采用非均布網(wǎng)格的一維彈性桿問(wèn)題收斂率對(duì)比Fig.9 Convergence comparison for the 1D elastic rod problem using non-uniform meshes

圖10 一維彈性問(wèn)題應(yīng)力解對(duì)比Fig.10 Comparison of stress solutions for the 1D elastic rod problem

圖11 方形區(qū)域勢(shì)問(wèn)題模型Fig.11 Description of the square region potential problem

圖11 中所示自然和強(qiáng)制邊界條件均由式(45)給定.圖12 為該問(wèn)題采用的有限元離散模型,對(duì)應(yīng)的節(jié)點(diǎn)數(shù)和平面三角形單元數(shù)分別為289,1089,4225 和512,2048,8192.

圖13 給出了該勢(shì)問(wèn)題的L2和H1誤差收斂結(jié)果,其中節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的L2和H1誤差收斂率分別為2 和1.9,與上節(jié)的理論收斂率基本一致,尤其是H1誤差的收斂率和精度都明顯高于有限元法,具有超收斂特性.圖14 為采用512 個(gè)三角形單元計(jì)算得到的梯度解的相對(duì)誤差對(duì)比,即,其中圖14(a)和圖14(b)分別表示有限元法和節(jié)點(diǎn)梯度光滑有限元法的誤差分布圖.由圖14 的對(duì)比結(jié)果清晰可見(jiàn),傳統(tǒng)線(xiàn)性三角形單元的梯度解為常數(shù)且在單元邊界處不連續(xù),因此節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的梯度解誤差遠(yuǎn)小于傳統(tǒng)有限元法,表明了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法梯度解的精度優(yōu)勢(shì).

4.3 懸臂梁?jiǎn)栴}

圖15 所示為一平面應(yīng)力懸臂梁?jiǎn)栴},其左端為強(qiáng)制位移邊界,右端受豎向力P=1000 作用,幾何與材料參數(shù)為:長(zhǎng)度L=4,梁截面高H=1,寬W=1,楊氏模量E=2.0×107,泊松比ν=0.3,慣性矩I=WH3/12.該問(wèn)題的解析解為

圖12 方形區(qū)域勢(shì)問(wèn)題的有限元離散模型Fig.12 Finite element discretizations for the square region potential problem

圖13 方形區(qū)域勢(shì)問(wèn)題收斂率對(duì)比Fig.13 Convergence comparison for the square region potential problem

圖14 方形區(qū)域勢(shì)問(wèn)題梯度解誤差對(duì)比Fig.14 Error comparison of the gradient solution for the square region potential problem

圖15 懸臂梁?jiǎn)栴}Fig.15 Description of the cantilever beam problem

該懸臂梁?jiǎn)栴}的有限元離散模型如圖16 所示,三種有限元離散對(duì)應(yīng)的節(jié)點(diǎn)數(shù)225,833 和3201,三角形單元數(shù)為384,1536 和6144.圖17 為該算例的位移和能量誤差收斂對(duì)比結(jié)果.與一維彈性桿問(wèn)題的結(jié)果一致,傳統(tǒng)線(xiàn)性有限元法的位移和能量誤差收斂率分別為2 和1,而節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的位移誤差收率也為2,但其能量誤差收斂率為2,比傳統(tǒng)有限元法高一階,具有更好的收斂特性.此外,圖17也表明節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的能量誤差明顯小于傳統(tǒng)有限元法.圖18 給出了采用384 個(gè)單元的剪應(yīng)力解的相對(duì)誤差云圖,結(jié)果進(jìn)一步顯示出節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的高精度應(yīng)力計(jì)算特點(diǎn).

圖16 懸臂梁?jiǎn)栴}有限元離散模型Fig.16 Finite element discretizations for the cantilever beam problem

圖17 懸臂梁誤差收斂率對(duì)比Fig.17 Convergence comparison for the cantilever beam problem

圖18 懸臂梁?jiǎn)栴}應(yīng)力解誤差對(duì)比Fig.18 Error comparison of the stress solution for the cantilever beam problem

4.4 三維厚壁圓筒勢(shì)問(wèn)題

為了驗(yàn)證節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的三維適用性,最后一個(gè)算例是圖19 所示的厚壁圓筒勢(shì)問(wèn)題,相應(yīng)的幾何參數(shù)為:內(nèi)徑為ri=1,外徑為ro=2,高度H=1.該問(wèn)題的解析解為

利用該問(wèn)題的對(duì)稱(chēng)性,計(jì)算中取四分之一模型進(jìn)行離散分析.圖20 為該三維勢(shì)問(wèn)題的有限元離散情況,對(duì)應(yīng)的節(jié)點(diǎn)數(shù)為225,1377 和9537,四節(jié)點(diǎn)四面體單元數(shù)為768,6144 和49 152.由于空間區(qū)域的構(gòu)造特點(diǎn),該問(wèn)題的有限元網(wǎng)格具有非均布特性.圖21 為采用圖20 中三種有限元離散網(wǎng)格計(jì)算得到的L2和H1誤差收斂率結(jié)果.與一維和二維算例的結(jié)果類(lèi)似,有限元法和節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的L2誤差收斂率都為2,節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法的H1誤差收斂率受非均布有限元網(wǎng)格影響,略低于理論值2,但其H1誤差的收斂率和精度仍然明顯高于傳統(tǒng)有限元法,驗(yàn)證了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法求解三維問(wèn)題的精度和有效性.圖22 給出了768個(gè)單元對(duì)應(yīng)的梯度解的相對(duì)誤差對(duì)比.結(jié)果再次證明,即使對(duì)于非均布有限元離散,節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法仍然在梯度解方面具有突出的精度優(yōu)勢(shì).

圖19 三維厚壁圓筒勢(shì)問(wèn)題模型Fig.19 Description of the 3D hollow cylinder potential problem

圖20 三維厚壁圓筒勢(shì)問(wèn)題有限元離散模型Fig.20 Non-uniform finite element discretizations for the 3D hollow cylinder potential problem

圖21 三維厚壁圓筒勢(shì)問(wèn)題收斂率對(duì)比Fig.21 Convergence comparison for the 3D hollow cylinder potential problem

圖22 三維厚壁圓筒梯度解誤差對(duì)比Fig.22 Error comparison of the gradient solution for the 3D hollow cylinder potential problem

5 結(jié)論

傳統(tǒng)有限元形函數(shù)在離散區(qū)域內(nèi)一般僅有C0連續(xù)性,在單元邊界和節(jié)點(diǎn)處上不存在一階和二階梯度,因而不能直接用于配點(diǎn)法.本文提出了有限元形函數(shù)的光滑梯度構(gòu)造理論,建立了具有與有限元形函數(shù)同樣連續(xù)性的一階和二階光滑梯度,進(jìn)而發(fā)展了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法.有限元形函數(shù)的一階光滑梯度依賴(lài)于在廣義梯度光滑理論框架內(nèi)的一階光滑梯度節(jié)點(diǎn)值的定義,以及選擇有限元形函數(shù)作為核函數(shù)對(duì)一階光滑梯度節(jié)點(diǎn)值進(jìn)行梯度光滑.對(duì)有限元一階光滑梯度進(jìn)行求導(dǎo),并用一階光滑梯度替換有限元形函數(shù)的標(biāo)準(zhǔn)梯度,即可得到有限元形函數(shù)的二階光滑梯度.文中詳細(xì)證明了線(xiàn)性有限元形函數(shù)的光滑梯度不僅滿(mǎn)足常規(guī)的一階梯度一致性條件,而且在均布網(wǎng)格條件下還滿(mǎn)足二階梯度一致性條件.與之對(duì)應(yīng)的是,線(xiàn)性有限元形函數(shù)的標(biāo)準(zhǔn)一階梯度只滿(mǎn)足一階梯度一致性條件.正是由于有限元形函數(shù)的光滑梯度滿(mǎn)足二階梯度一致性條件,理論分析表明,基于有限元形函數(shù)光滑梯度的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法具有二階精度,即L2誤差和H1或能量誤差的收斂率均為2,與傳統(tǒng)伽遼金有限元法相比呈現(xiàn)超收斂特性.文中通過(guò)典型算例系統(tǒng)驗(yàn)證了節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法在梯度或應(yīng)力求解方面的超收斂特性和顯著精度優(yōu)勢(shì).本文的節(jié)點(diǎn)梯度光滑有限元配點(diǎn)法采用了線(xiàn)性有限元形函數(shù),對(duì)于高次有限元形函數(shù)及高階光滑梯度還需要進(jìn)一步研究.