節點梯度光滑有限元配點法1)

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

(廈門大學土木工程系,福建廈門 361005)

(廈門市交通基礎設施智能管養工程技術研究中心,福建廈門 361005)

引言

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

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

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

1 配點法控制方程

本文考慮穩態勢問題和彈性力學問題等典型二階問題,其控制方程的強形式可以表示為

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

其中n表示在自然邊界Γt的單位外法線向量.對于平面應力彈性力學問題,上述算子可表示為

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

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

其中K,d和f分別表示剛度矩陣,節點系數向量和力向量

由式(10)可見,配點法需要計算形函數的二階梯度,而標準有限元形函數在單元邊界處不存在二階梯度,因此無法直接用于配點法.此外,雖然通過特殊構造方式可以建立一些特殊單元的高階光滑有限元形函數,但是目前仍然缺乏構造一般單元高階光滑形函數的通用有效方法[1].本文則是通過引入無網格法中的梯度光滑技術,建立了有限元形函數的一階和二階光滑梯度,進而發展了相應的節點梯度光滑有限元配點法.

2 有限元形函數的光滑梯度構造

2.1 有限元形函數光滑梯度構造理論

考慮圖1 所示的有限元網格,?L表示第L個單元對應的空間區域,節點xI對應的有限元形函數是NI(x),NI(x)的影響域為

式中MI表示共享節點xI的單元個數.方便起見,定義下面的節點組和單元組

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

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

式中φ(x,y)為梯度光滑核函數,NI,i(x)為傳統的標準一階梯度或導數,即NI,i(x)=?NI(x)/?xi,xi為x的第i個坐標分量.由于有限元形函數的梯度在單元邊界處沒有定義,因此式(14)并不能直接用于有限元形函數的光滑梯度構造.

為了構造有限元形函數的光滑梯度,這里首先引入有限元形函數在節點處的梯度值.為此,在式(14)中選取如下的梯度光滑核函數

若進一步考慮常應變單元,例如一維兩節點單元、二維三節點三角形單元和三維四節點四面體單元,則NI,i(x)在每個單元內均為常數,式(16)可以簡化為

其中VL表示單元?L的長度、面積或體積.

在式(17)定義的有限元形函數的節點梯度基礎上,進一步選擇梯度光滑核函數φ(x,xJ)=NJ(x),則根據式(14)對應的離散形式,可得如下的有限元形函數一階光滑梯度

接下來,對式(18)進行求導,并用式(18)定義的一階光滑梯度代替傳統一階梯度,最后可得有限元形函數的二階光滑梯度

至此,已經構建有限元形函數的二階光滑梯度.下面通過一維和二維線性有限元形函數一階和二階光滑梯度的構造過程來說明有限元形函數光滑梯度構造理論.

2.2 線性有限元形函數光滑梯度構造

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

其對應的一階梯度或導數NI,x(x)為

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

另一方面,根據式(14),有限元形函數在節點處的一階光滑梯度為

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

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

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

圖2 和圖3 給出了線性有限元形函數一階和二階光滑梯度的構造過程.由圖2 和圖3 可見,雖然線性有限元形函數的梯度在節點處不連續,而且二階梯度不存在,但由本文所提有限元形函數梯度光滑理論構造所得的一階和二階光滑梯度仍然是連續函數,可以直接用于配點法求解.

對于二維和三維情況,考慮三節點三角形單元和四節點四面體單元,其對應的形函數為標準的線性形函數[1],這里不再贅述.二維和三維有限元形函數的一階和二階光滑梯度也可按照上述步驟進行構造.圖4 給出了基于平面三角形線性單元構造一階光滑梯度和和二階光滑梯度的過程.與一維情況相同,線性三角形單元的形函數在節點和邊界處的一階梯度均不連續,不存在二階梯度.另一方面,基于本文的有限元形函數光滑梯度構造理論,可以方便地構造圖5 所示的一階和二階光滑梯度.此外,圖2～圖5 同時表明,與有限元形函數和其標準梯度不同,有限元形函數的一階和二階光滑梯度與形函數本身具有相同的連續性,但光滑梯度的影響域則大于形函數的影響域.

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

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

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

3 節點梯度光滑有限元配點法及其精度分析

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

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

3.1 線性有限元形函數光滑梯度的一致性條件

對于傳統線性有限元,其形函數及其一階梯度的一致性或完備性條件為[1]

其中xI j為節點xI的第j個坐標分量,δi j為克羅內克符號.注意到有限元形函數的一致性條件取決于形函數的階次.例如,線性有限元形函數僅滿足式(26)～式(29)定義的線性一致性條件,由于其二階梯度不存在,故不滿足更高一次的二階一致性條件.但是對于有限元形函數的光滑梯度,接下來證明其不僅滿足線性一致性條件,在均布網格情況下還滿足二階一致性條件.

對于均勻有限元網格,由形函數的周期重復性可知,節點光滑梯度滿足如下條件[24-25]

其中k=0,1 或2.首先考慮k=2 的情況,利用式(32),式(33)變為

再者,對于k=1,式(33)變為

式(34)～式(36),證明了線性有限元形函數的二階光滑梯度滿足二階梯度一致性條件.與之形成鮮明對比的是,標準線性有限元形函數不存在二階梯度.

3.2 節點梯度光滑有限元配點法精度分析

為不失一般性,以平面勢問題為例,在均勻網格條件下對節點梯度光滑有限元配點法的精度進行解析研究.節點梯度光滑有限元配點法的精度可以方便地采用如下的截斷誤差進行度量[24,30]

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

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

同時注意到在均布網格條件下,有限元形函數的光滑二階導數具有周期性和偶對稱的特點,因此當n為奇數時,有

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

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

4 算例

本節通過一維、二維及三維算例驗證節點梯度光滑有限元配點法的收斂性和精度,其中FEM 和FEC 分別表示傳統的伽遼金有限元法和本文所提的節點梯度光滑有限元配點法.為了和傳統有限元法進行比較,算例中采用如下誤差形式

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

4.1 一維彈性桿問題

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

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

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

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

4.2 二維矩形區域勢問題

考慮如圖11 所示的正方形勢問題區域,其邊長L=1,對應勢問題的解析解為

圖8 采用均布網格的一維彈性桿問題收斂率對比Fig.8 Convergence comparison for the 1D elastic rod problem using uniform meshes

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

圖10 一維彈性問題應力解對比Fig.10 Comparison of stress solutions for the 1D elastic rod problem

圖11 方形區域勢問題模型Fig.11 Description of the square region potential problem

圖11 中所示自然和強制邊界條件均由式(45)給定.圖12 為該問題采用的有限元離散模型,對應的節點數和平面三角形單元數分別為289,1089,4225 和512,2048,8192.

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

4.3 懸臂梁問題

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

圖12 方形區域勢問題的有限元離散模型Fig.12 Finite element discretizations for the square region potential problem

圖13 方形區域勢問題收斂率對比Fig.13 Convergence comparison for the square region potential problem

圖14 方形區域勢問題梯度解誤差對比Fig.14 Error comparison of the gradient solution for the square region potential problem

圖15 懸臂梁問題Fig.15 Description of the cantilever beam problem

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

圖16 懸臂梁問題有限元離散模型Fig.16 Finite element discretizations for the cantilever beam problem

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

圖18 懸臂梁問題應力解誤差對比Fig.18 Error comparison of the stress solution for the cantilever beam problem

4.4 三維厚壁圓筒勢問題

為了驗證節點梯度光滑有限元配點法的三維適用性,最后一個算例是圖19 所示的厚壁圓筒勢問題,相應的幾何參數為:內徑為ri=1,外徑為ro=2,高度H=1.該問題的解析解為

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

圖19 三維厚壁圓筒勢問題模型Fig.19 Description of the 3D hollow cylinder potential problem

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

圖21 三維厚壁圓筒勢問題收斂率對比Fig.21 Convergence comparison for the 3D hollow cylinder potential problem

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

5 結論

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