基于赫林格-賴斯納變分原理的一致高效無網格本質邊界條件施加方法1)

吳俊超 吳新瑜 趙珧冰 王東東

*(華僑大學土木工程學院,福建省智慧基礎設施與監測重點實驗室,福建廈門 361021)

?(廈門大學土木工程系,福建省濱海土木工程數字仿真重點實驗室,福建廈門 361005)

引言

無網格法[1-7]是指基于節點信息構造離散形函數的一類數值方法的總稱.該類方法通常具有高階光滑、全域協調的形函數,形函數構造過程不依賴于網格單元的拓撲信息,適用于大變形分析[8]、薄板殼高階問題[9]及裂紋擴展模擬[10]等.然而,與傳統有限元法相比,高階連續光滑的特點導致無網格形函數在離散節點上通常不具有插值性,在求解過程中難以直接施加本質邊界條件[11-12].其次,無網格形函數通常為有理式,并且形函數影響域高度重疊,導致形函數在背景積分單元上為分段的有理式.在伽遼金法的求解過程中,傳統高斯積分法無法精確數值積分由形函數組成的剛度矩陣和力向量,導致伽遼金弱形式不滿足積分約束條件或稱變分一致性條件[13-14],無法保證計算精度和最優誤差收斂率.

為了使無網格法能夠直接施加本質邊界條件,許多學者構造了諸多具有插值性的無網格近似方法,例如奇異權函數法[15]、插值最小二乘法[16-17]、復變量移動最小二乘法[18]、廣義移動最小二乘法[19]、變換法[20]等.然而,這類方法不是建立在變分原理基礎上,并不能保證節點之間位移邊界條件施加精度和無網格法的變分一致性.此外,對于滿足積分約束條件的無網格數值積分方法,例如穩定節點積分法[13]、一致性積分法[21-22]、變分一致積分法[23]、嵌套子域積分法[24]、再生光滑梯度積分法[25]等,在計算過程中采用形函數的光滑梯度替代傳統無網格形函數的導數,在保證無網格法的計算精度和最優誤差收斂率的同時提高了計算效率.但其本質邊界條件仍需要具有變分一致性的方法進行施加[26-27].

在伽遼金無網格法中,拉格朗日乘子法[1]是常采用的一種變分一致本質邊界條件施加方法.該方法需要在整體剛度矩陣上增加額外的自由度離散拉格朗日乘子.當采用滿足變分一致無網格數值積分方法時,拉格朗日乘子的自由度需要與光滑梯度構造過程中的數值積分點相一致.但是,過多的拉格朗日乘子自由度將增加整體剛度矩陣奇異性,不適用于高階基函數的無網格法.Lu等[28]根據拉格朗日乘子的物理意義,采用位移自由度離散拉格朗日乘子,無需增加額外自由度,并稱之為修正變分原理法.但該方法施加本質邊界條件過程中的修正變分項降低了剛度矩陣的正定性,計算精度不高.Zhu和Atluri[29]將罰函數法引入無網格法中施加本質邊界條件,該方法格式簡潔易實現,但其形式并不滿足變分的一致性,求解精度依賴人工參數,穩定性較低.Fernández-Méndez和Huerta[11]采用尼茲法(Nitsche method)施加無網格法本質邊界條件,該方法可視為修正變分原理與罰函數法的結合,修正變分項保證了變分一致性,而罰函數項恢復了剛度矩陣的正定性,提高了求解的穩定性,但該方法仍然需要選擇合適的人工參數來保證數值解的精度.

本文針對滿足積分約束條件的伽遼金無網格法,提出了一種具有變分一致且高效的本質邊界條件施加方法.該方法以赫林格-賴斯納(Hellinger-Reissner,HR)變分原理[30]為基礎,采用混合離散的方式近似弱形式中的位移和應力.其中,位移采用傳統無網格形函數進行離散,而應力則在每個背景積分單元中近似為局部多項式.此時,近似應力中多項式的系數可在背景積分單元中采用位移自由度表示,無需額外增加整體自由度.不僅如此,應力的系數表達式中包含再生光滑梯度表達式,在再生光滑梯度的理論框架下,光滑梯度構造過程中可采用優化的數值積分點以減少全域無網格形函數的計算量,計算效率高.整體的離散平衡方程具有與尼茲法相類似的形式,其中尼茲法中的修正變分項分別引入了再生光滑梯度和無網格形函數離散應力和位移,整體求解過程無需計算傳統無網格形函數的導數.而在赫林格-賴斯納變分原理的框架下,弱形式中自帶穩定項,無需增加額外的穩定項,且穩定項中不存在人工參數,消除了對人工參數的依賴性.文中通過二維彈性力學算例驗證了所提方法的精度、效率和誤差收斂性.

1 赫林格-賴斯納變分原理

不失一般性,在求解域 Ω 內考慮如下彈性力學控制方程

式中,W為彈性體的余能密度函數,其與應力之間的關系為

其中Cijkl為四階彈性張量.

對上式進行變分可得與之相對應的弱形式

2 位移離散與再生核近似

在赫林格-賴斯納變分原理的弱形式中,位移與應力可分別采用不同的近似方案進行離散.這里,位移采用基于再生核近似[29]的無網格形函數進行離散.首先,在求解域 Ω 上布置一系列無網格節點.此時,位移分量ui的近似表達式可表示為

式中,diI為無網格節點xI上的位移分量節點系數,ΨI為與之相對應的無網格形函數,根據再生核近似理論[2],無網格形函數具有如下表達式

其中,p為基函數向量,以二維p階基函數為例,此時基函數向量為

φ為核函數,其影響域也決定了無網格形函數的影響域.在二維情況下,核函數的影響域通常為圓形或矩形,本文采用基于三次樣條函數、具有矩形影響域的核函數

式中,rx=|xI-x|/sx,ry=|yI-y|/sy,sx和sy分別為x和y方向的影響域尺寸.φ為三次樣條函數[1]

無網格形函數表達式(8)中的c為修正函數向量,該表達式可通過滿足再生條件[2]確定

將式無網格形函數表達式(8) 代入再生條件(13)中,可得

式中M稱為矩量矩陣

將式(14)、式(15)代入式(8)中,可得到再生核無網格形函數的表達式

圖1為二維二次基函數無網格形函數,從圖中可以看出,無網格形函數在全域上連續光滑.但在無網格節點處,無論節點位于域內或邊界處,形函數都不具有插值性,這是伽遼金無網格法難以施加本質邊界的主要原因.同時,相較于有限元形函數,無網格形函數并不具有顯示表達式,形函數計算復雜耗時,無網格形函數的導數更是如此.

圖1 無網格形函數Fig.1 Meshfree shape functions

3 應力離散與再生光滑梯度近似

表述方便起見,引入沃伊特標記(Voigt notation),將近似應力 σh采用矩陣向量形式改寫為

4 赫林格-賴斯納變分原理下的本質邊界條件施加方法

基于赫林格-賴斯納變分原理的弱形式(4)中已考慮了本質邊界條件.前兩節已分別采用再生核無網格形函數和背景積分單元上的分段多項式近似式(4)中的位移和應力,進一步將近似的位移分量表達式(7)和應力分量表達式(17)代入弱形式(6)中,有

其中,δd為位移變分節點系數向量,f為外力向量.將式(19)代入式(33)中,式(33)等式左邊可化簡為

式中最后一個等式引入了關系式(20)～式(26)進行化簡,詳細的推導過程可參考文末附錄A.將簡化后的式(34)代入式(33),即可得到離散的平衡方程

從式(35)和式(36)可以看出,在赫林格-賴斯納變分原理的框架下,本質邊界條件的施加過程與尼茲法[11]具有相類似的格式.清晰起見,這里給出基于尼茲法的無網格法離散平衡方程

其中,K和f與式(36) 中的剛度矩陣K和力向量f相同.Kv和fv是基于變分原理施加本質邊界條件產生的剛度和力向量,兩者保證了邊界條件施加方法的變分一致性,其具體表達式為

再者,式(40)中的Kp和fp對應于尼茲法中采用罰函數法施加穩定項的剛度矩陣和力向量

式中 1為單位對角矩陣

該部分通過設置合適的罰因子 α[31]來保證剛度矩陣的正定性,提高求解的穩定性.

對照本文所提方法和尼茲法的無網格離散方程,即式(35)和式(40),可見式(35)中的對應于尼茲法中利用變分原理施加本質邊界條件,但其應變離散過程采用式(25)中的光滑梯度替代傳統無網格形函數梯度,而剛度矩陣的正定性則是通過施加得以保證.與尼茲法相比,本文所提方法的穩定項表達式相對較為復雜,但在求解全過程中,無需計算復雜耗時的無網格形函數梯度,整體過程計算效率優于尼茲法,且穩定項也不包括任何人工參數.

在引入數值積分求解式(35)的過程中,為保證離散弱形式的變分一致性,數值積分在求解的各個過程也需保持一致.如圖2 所示,所提方法需要兩套數值積分點,和中的數值積分方案需要與,,f,和保持一致,為了保證理論誤差收斂率,該數值積分方案的積分精度在文獻[25]有詳細討論.另一方面,變分一致性要求式(20)中G與式(35)中的K采用相同數值積分方案.在所提方法中,由于再生光滑梯度已滿足積分約束條件,再生光滑梯度也為低階多項式,此時數值積分方案采用傳統有限元中對應階次的高斯積分即可保證數值精度.

圖2 背景積分單元示意圖和優化的數值積分方案Fig.2 Illustration of background integration cells and optimized quadrature rules

此外,構造再生光滑梯度需要計算積分點處無網格形函數.為了減小形函數的計算量、提高計算效率,這里采用無網格再生光滑梯度積分法[25],利用積分點在背景積分單元間的共享特性,從全域上優化整體求解過程中積分點的數量,提高計算效率.本文采用二維二次、三次基函數的優化數值積分方案,具體的數值積分點的位置與權重詳見附錄B.

5 數值算例

本節首先通過分片試驗,驗證了采用傳統高斯積分法和再生光滑梯度積分法的不同本質邊界條件施加方法下是否滿足積分約束條件.然后,通過典型的懸臂梁問題和帶孔無限大平板問題來驗證采用二次和三次基函數時,所提算法的計算精度和效率.在精度分析中,采用如下的位移誤差(L2-Err)和能量誤差(He-Err)

表述清晰起見,下面的討論中用“GI”和“RKGSI”表示高斯積分法和再生光滑梯度積分法,用“penalty”、“LM”和“Nitsche”表示罰函數法、拉格朗日乘子法和尼茲法這三種施加本質邊界條件的方法.例如,本文所提的基于赫林格-賴斯納變分原理的本質邊界條件施加方法簡寫為“RKGSI-HR”.為了保證拉格朗日乘子法求解的穩定性,拉格朗日乘子均采用線性有限元形函數進行離散.當采用二次基函數時,GI 的域內積分采用13 點高斯積分,邊界積分采用3 點高斯積分;當采用三次基函數時,GI 采用16 點高斯積分,邊界積分采用5 點高斯積分.

5.1 分片試驗

首先,采用線性、二次和三次彈性力學分片試驗驗證不同無網格本質邊界條件施加方案的變分一致性.分片試驗的求解域為邊長等于1 的正方形,求解域的四邊施加本質邊界條件,邊界條件的預設值根據如下所示的精確解施加

其中,n=1,2,3 代表線性、二次和三次分片試驗.分片試驗采用如圖3 所示的非均布節點離散求解域,針對二次基函數的無網格近似,采用線性和二次分片試驗進行測試,核函數相對影響域為2.5;而三次基函數的無網格近似采用二次和三次分片試驗,核函數的相對影響域為3.5.

圖3 分片試驗無網格離散模型Fig.3 Meshfree discretization for the patch test

表1、表2 分別為具有二次、三次基函數的無網格法分片試驗結果,從表中可以看出,傳統高斯積分法由于不滿足積分約束條件,即使采用高階高斯積分也不能通過分片試驗.采用滿足積分約束條件的再生光滑梯度積分法時,罰函數法不具有變分一致性,無法通過分片試驗.而拉格朗日乘子法由于其拉格朗日乘子采用線性形函數進行離散,無法與再生光滑梯度相匹配,也不能通過分片試驗.當采用尼茲法和本文所提基于赫林格-賴斯納變分原理的本質邊界施加方法時,RKGSI 可通過分片試驗.

表1 二次基函數無網格法分片試驗結果Table 1 The results of patch test with quadratic basis functions

表2 三次基函數無網格法分片試驗結果Table 2 The results of patch test with cubic basis functions

5.2 懸臂梁問題

考慮圖4 所示懸臂梁問題,懸臂梁求解域的長和寬分別為L=48,D=12,懸臂梁楊氏模量E=3×106、泊松比 ν=0.3.懸臂梁的右端沿y軸正方向施加荷載P,而左端為固定支座.根據平面應力假設和圣維南原理,該問題的解析解[32]為

圖4 懸臂梁問題模型Fig.4 Description of the cantilever beam problem

與之相應的應力分量為

為保證一致性,該問題的本質邊界條件和自然邊界條件分別根據解析解式(47)和式(48)進行施加.

懸臂梁求解域采用圖5 所示的四個疏密程度不同的節點離散研究計算誤差收斂特性,算例采用具有二次基函數的無網格近似,相應的核函數相對影響域為2.5.圖6為該問題的位移誤差和能量誤差收斂率結果,結果表明所提RKGSI-HR 方法的精度與RKGSI-Nitsche 相當,均達到理論收斂率,但RKGSIHR 方法無需人工參數.傳統高斯積分法、RKGSILM和RKGSI-penalty 均不具有變分一致性,不能達到理論誤差收斂率.

圖5 懸臂梁問題節點離散Fig.5 Meshfree discretizations of the cantilever beam problem

圖6 懸臂梁問題誤差對比Fig.6 Error comparison for the cantilever beam problem

圖7為該問題的效率對比,其中圖7(a)為整體求解過程中的效率分析.從圖7(a) 中可以看出,RKGSI 整體的計算效率要優于GI,GI 所用CPU 耗時大概為RKGSI 的2.5 倍.這是由于RKGSI 整體無需計算傳統無網格形函數導數,同時采用優化的數值積分點減小形函數的計算量,兩方面共同作用的結果.此外,從該圖可以看出拉格朗日乘子法施加邊界條件時,整體的計算效率要略高于其他方法,這是由于拉格朗日乘子法需要增加自由度離散拉格朗日乘子,導致整體線性代數的計算量增加.為了更好說明本文所提RKGSI-HR 的計算效率,圖7(b)對比了不同本質邊界方法施加過程的效率對比.該圖對比了施加過程計算形函數及梯度和組裝相應的剛度矩陣和力向量所用時間,從圖中可以看出,罰函數法和拉格朗日乘子法由于計算過程僅需要無網格形函數本身,兩種方法的計算形函數耗時相同且優于其余兩種方法.相比于傳統尼茲法,基于赫林格-賴斯納變分原理施加本質邊界條件過程也無需計算形函數梯度,但需要構造相應的光滑梯度,此部分所用施加大約為罰函數法和拉格朗日乘子法的1.6 倍.而在組裝力向量方面,赫林格-賴斯納變分原理法和尼茲法具有相當的效率.拉格朗日乘子法由于采用有限元形函數離散拉格朗日乘子,其計算效率最高.需要指出的是,罰函數法和拉格朗日乘子法并不能保證理論誤差收斂率,且從整體看拉格朗日乘子法的計算效率并不優于其他方法.綜合本算例的精度分析和效率分析可得,RKGSI-HR 不僅保證計算精度和最優理論收斂性,并且相較于傳統尼茲法具有更高的計算效率.

圖7 懸臂梁問題效率對比Fig.7 Efficiency comparison for the cantilever beam problem

5.3 帶孔無限大平板問題

考慮圖8 所示帶孔無限大平板問題,其中板的中心有一半徑為a=1 的小孔,平板的無窮遠處沿x軸方向施加均布荷載T=1000.平板材料參數為:楊氏模量E=3×106、泊松比 ν=0.3.根據Michell解[32]可得該問題的解析解為

其中 κ和μ為

與之相對應的應力表達式為

根據帶孔無限大平板的對稱性,如圖8 所示,取邊長為b=5 的四分之一方形域進行研究,求解域的左端與下端(圖中紅線邊界)約束法向位移,求解域上端、右端和圓孔邊界(圖中藍線邊界)根據解析應力表達式(51)施加自然邊界條件.該問題的誤差收斂率測試采用圖9 所示的112 個、403 個、1525 個、5929 個節點進行離散,無網格近似采用三次基函數,核函數的相對影響域為3.5.

圖8 帶孔無限大平板問題模型Fig.8 Description of the plate with a hole problem

圖9 帶孔無限大平板問題無網格離散Fig.9 Meshfree discretizations of the plate with a hole problem

圖10為位移和能量誤差收斂的結果,結果再次驗證了所提赫林格-賴斯納變分原理法能夠保證理論誤差收斂率.同樣地,具有變分一致性的RKGSINitsche 也達到了理論誤差收斂率,而傳統高斯積分法即使在節點較為稀疏時,也難以保證收斂效果.該問題中,雖然RKGSI-LM 不能滿足變分一致性,但拉格朗日乘子法保持著較高的計算精度,其誤差并未影響該方法達到理論誤差收斂率.圖11為該問題的效率對比,與同樣具有變分一致性的RKGSINitsche 方法相比,RKGSI-HR 無需計算傳統形函數導數,具有更高的計算效率.最后,圖12為該問題的應力云圖,圖中清晰表明具有RKGSI 的計算精度要高于GI,GI 的應力云圖均出現不同程度的振蕩.

圖10 帶孔無限大平板問題誤差對比Fig.10 Error comparison for the plate with a hole problem

圖11 帶孔無限大平板問題效率對比Fig.11 Efficiency comparison for the plate with a hole problem

圖12 帶孔無限大平板問題 σxx應力云圖Fig.12 Contour plot of stress σxxfor the plate with a hole problem

6 結論

本文以赫林格-賴斯納變分原理為基礎,提出了一種滿足積分約束條件的變分一致高效本質邊界條件施加方法.該方法采用混合離散近似赫林格-賴斯納變分原理弱形式中位移和應力,其中位移采用傳統無網格形函數進行離散,而應力則在每個背景積分單元上近似為對應階次的多項式.在赫林格-賴斯納變分原理的框架下,該方法的離散平衡方程具有與傳統尼茲法相類似的格式,可視為與再生光滑梯度積分法相配套的新型尼茲法.與傳統尼茲法相比,所提方法的修正變分項采用無網格形函數和再生光滑梯度進行混合離散,在保證了變分一致性的同時避免了復雜耗時的形函數導數計算,明顯提高了計算效率;而對應于尼茲法中的穩定項則直接源于赫林格-賴斯納變分原理弱形式,無需額外增加穩定項,更重要的是穩定項中不包含任何人工參數,有效消除了尼茲法中的人工參數依賴性問題.文中通過典型算例系統地驗證了所提基于赫林格-賴斯納變分原理施加本質邊界條件方法的變分一致性、計算精度和計算效率.

附錄A

本附錄中詳細推導了第4 小節式(34)中最后一個等式的推導過程,該等式為

引入關系式(20)～式(26),式(A1)中等式的左端各項簡化過程如下所示

綜合式(A2)～式(A6),可得等式(A1)成立.需要注意的是,從式(A4)的化簡過程可以看出,的表達式雖然直觀上看不對稱,但從原表達式中可知其實是對稱矩陣.

附錄 B

本附錄中列出了本文所提方法在數值實現過程中所采用的積分點位置和權重,其中附表B1和附表B2 分別為二次基函數情況和三次基函數情況的優化數值積分方案.表中“”為KIJ,G所采用的數值積分點,所采用的數值積分點.并且 ξ,η和γ為三角形參數空間坐標,w和wB分別為三角形域內積分權重和邊界積分權重.

表B1 二次基函數無網格法優化的數值積分方案Table B1 The optimized quadrature rules for quadratic basis function

表B2 三次基函數無網格法優化的數值積分方案Table B2 The optimized quadrature rules for cubic basis function