基于圖像四叉樹的改進型比例邊界有限元法研究1)

孫立國 江守燕 杜成斌

(河海大學工程力學系,南京 210098)

引言

比例邊界有限元法(scaled boundary finite element method,SBFEM)具有僅需離散結構域邊界使得計算問題的維度降低一維、可存在懸掛節點而不會帶來額外的計算問題等優點,該方法已廣泛應用于結構的計算分析中[1-2],特別是在斷裂損傷問題領域[3-11],SBFEM 具有其獨特的優勢,可自動地解析求出縫端的應力強度因子.

由于SBFEM 本身的固有優勢,SBFE 網格可為任意的多邊形單元[12],特別地,SBFEM 與四叉樹網格的結合[13-17]使得前處理工作具有快速、高效的特點,粗細網格過渡十分方便,并能實現自動化的網格剖分,大大地減輕了網格剖分的負擔.然而,上述這些研究在模擬界面演化問題(如裂紋擴展問題)時,雖然能夠使得網格重剖分工作達到最小化,但局部的網格重剖分工作仍不可避免.

Natarajan 等[18]和大連理工大學李建波等[19-23]結合擴展有限元法和SBFEM 的優點,提出了擴展比例邊界有限元法(X-SBFEM)的概念,該方法采用擴展有限單元模擬裂紋貫穿單元、SBFE 模擬裂尖單元,從而避免了傳統擴展有限元法需要構造裂尖單元復雜改進函數以及需要采用特殊的數值積分技術進行裂尖單元剛度矩陣求解的弊端.江守燕等[24-25]進一步完善和改進了該方法,借助于擴展有限元法的基本思想,將擴展有限元法中所有改進單元轉化為SBFE 的子域,被裂紋切割的改進單元無需引入特殊的富集函數,取得了較好的模擬效果.Huynh 等[26]通過在擴展有限元法中耦合SBFEM 的基本思想,緩解了裂紋生長模擬中網格細化過程的懸掛節點問題.

本工作在文獻[24-25]的研究基礎上,結合四叉樹網格的優勢,提出基于圖像四叉樹的改進型比例邊界有限元法,通過結構域圖像實現結構域網格的全自動剖分,且結構內不連續處局部范圍的網格粗細過渡方便、快速.

1 基于圖像的四叉樹網格

四叉樹生成的網格可以直接轉換為SBFEM 子域,而不需要像傳統有限元那樣區分網格的角節點和懸掛節點.由于SBFEM 的半解析性質,它只需對結構邊界進行離散化,徑向具有解析解,只需比例中心對所有邊界都有良好的可視性,即使邊界大小不同,也不會影響數值精度.因此,四叉樹網格可以完美地用作SBFEM 子域,基于圖像的四叉樹網格剖分技術可實現SBFEM 網格的快速、全自動剖分[27].下面具體介紹基于圖像的四叉樹網格剖分的關鍵步驟.

(1)獲取待四叉樹網格剖分的結構圖像.圖1(a)顯示了一個像素為N×M的圖像,圖像內部包含一個夾雜,基體、夾雜以及兩者之間的界面通過不同顏色區分.使用四叉樹算法對其進行網格剖分,圖像顏色矩陣QN×M用于存儲像素的所有顏色.

(2)初始化四叉樹網格.如圖1(b)所示,通過將圖像分割為多個邊長為smax的單元來初始化四叉樹.

(3)如果單元中最大和最小顏色強度之間的差異大于顏色閾值,則該單元被遞歸地進一步劃分為四個大小相等的單元.該過程遞歸執行,直到所有單元滿足均勻性標準或達到最小邊緣長度.然后得到最終的四叉樹網格,見圖1(c).

在四叉樹分解過程中使用了平衡四叉樹分解,該分解僅產生16 種可能的單元模式,通過考慮旋轉對稱性,只有六種類型的唯一節點排列是可能的,見圖1(d).在SBFEM 中,這些懸掛節點可以直接地視為新多邊形單元的節點,無需任何特殊處理.

圖1 基于圖像的四叉樹網格剖分Fig.1 Image-based quadtree cell generation

2 改進型比例邊界有限元法

改進型SBFEM 的主要思想是: 在SBFEM 的理論框架下,基于水平集法表征材料的內部裂紋面,含不連續裂紋面的子域可通過節點水平集函數識別,裂紋擴展時無需進行網格重剖分,界面的幾何特征通過SBFEM 子域的附加自由度表征.

2.1 裂紋的水平集描述

φ(x)Φk(x)(k=1,2)

如圖2 所示,兩個水平集和用于描述裂紋[28].通常假設裂紋尖端水平集函數Φk(x)(k=1,2)與裂紋面水平集函數 φ(x) 正交.函數φ(x)可以用符號距離函數表示為

圖2 裂紋的水平集表征Fig.2 Crack description based on two level set functions

式中,x是考察點的坐標;x*是考察點在裂紋表面上的投影點;n是垂直于裂紋表面的外法向向量;sign(x)是符號函數.函數 Φk(x)(k=1,2) 可以定義為

式中,xk是第k個裂紋尖端的坐標,t是第k個裂紋尖端的單位切向矢量.

2.2 改進子域的識別

基于水平集函數,可以確定四叉樹網格與裂紋的切割關系.如圖3 所示,可以觀察到兩種類型的關系: (1)四叉樹單元被裂紋完全切割,即域 Ωcc;(2)四叉樹單元被裂紋部分切割,單元包含一個裂紋尖端,即域 Ωcp.改進子域識別的具體流程如下.

圖3 四叉樹網格與裂紋間的切割關系Fig.3 Relationship between the quadtree cell and a crack

步驟1: 遍歷求解域內所有節點,根據式(1)和式(2)求出所有節點的兩個水平集函數值.

步驟2: 針對某一個含有N個節點的四叉樹單元,根據式(2)計算的每個節點的波前水平集函數值,判斷裂尖點是否位于該單元內.若該單元包含裂尖點,則屬于域 Ωcp.若不包含裂尖點,則轉至步驟3 進一步判斷該單元是否屬于域 Ωcc.

步驟3: 若單元屬于域 Ωcc,需要滿足條件: 其中M(0 ＜M＜N)個節點的裂紋面水平集函數 φ(x)＜0,則剩余的 (N-M) 個節點的裂紋面水平集函數φ(x)≥0 .若不滿足前述條件,則該單元為常規單元.

步驟4: 對于求解域內的所有四叉樹單元,重復步驟2、步驟3,即可識別出所有改進單元.

2.3 改進的比例邊界有限元子域

如圖4(a) 所示,完全由裂紋切割的單元視為SBFEM 的兩個獨立子域.與常規SBFEM 不同的是,通過在裂紋與單元邊界的交點處引入額外的節點,將原四叉樹單元一分為二.在本研究中,SBFE 子域將通過在交叉點處引入額外的自由度而得到富集,單元并沒有物理分區,這樣就無需改變網格以遵從幾何界面,尤其是當裂紋擴展時無需進行局部的網格重剖分.

如圖4(b)所示,對于包含裂紋尖端的單元,選擇裂紋尖端周圍的若干單元形成超級單元.根據裂尖所在單元找出裂尖單元周邊的單元,如周邊單元無法形成一個方形域,則繼續尋找,直至包圍裂尖的單元可形成一個方形域為止,超級單元作為一個SBFEM子域.超級單元內部的節點位移可通過SBFE 位移場近似獲得.被裂紋切割的單元邊界并沒有物理分開,而是增加了額外的自由度來表示不連續性.

圖4 改進的SBFE 子域Fig.4 Enriched SBFE subdomain

2.4 改進子域剛度矩陣的求解

如圖5 所示,具有4 個角節點和2 個懸掛節點的四叉樹單元(SBFE 子域)的邊被劃分為多個線單元,即L1,L2,L3,L4,L5,L6.單元局部坐標系ξ-η的 ξ 方向由比例中心指向子域邊界,η 方向沿著子域環向邊界.比例中心O位于四邊形單元的中心,對所有邊界均具有較好的可見性.值得注意的是,與有限元不同,在SBFEM 中,懸掛節點不會帶來額外的計算困難,它只是將一條線單元分成兩條.

圖5 比例邊界有限元的四叉樹網格Fig.5 SBFE quadtree cell

SBFEM 位移模式可以表示為

式中,N(η) 為形函數矩陣,u(ξ) 為徑向的位移函數.

由位移模式式(3),可得應變場為

線性算子矩陣L可表示為

式中,|J| 為雅克比矩陣行列式,可表示為

根據胡克定律,應力場可表示為

式中,D為本構矩陣且

根據虛功原理,可得SBFE 方程為

載荷向量為

系數矩陣 (Ei,i=0,1,2) 為

為求解式(12),首先變換該式為一個一階微分方程

式中,哈密頓矩陣Z為

通過Schur 分解,可得

其中,Schur 矩陣S和變換V可表達為

式中,Sn和Sp為與哈密頓矩陣Z的負特征值和正特征值對應的上三角矩陣;Vu和分別為有限域和無限域的模態位移;Vq和分別為與有限域和無限域相應的模態力.對于有限域問題,僅Sn在比例中心處產生有限位移,u(ξ) 可以表示為

式中,c為依賴于邊界條件的積分常數.根據節點位移ub=u(ξ=1),積分常數可表達為

SBFE 子域的邊界節點力為

因此,SBFE 子域的剛度矩陣可通過下式求解得到

如圖6(a)所示,四叉樹單元被裂紋完全切割,切割面為位移不連續界面,在局部坐標系下,該單元包含四個角節點(節點編號: 1,2,3,4)和兩個懸掛節點(節點編號: 5,6),由于被裂紋切割,該單元還包括四個虛節點(節點編號: 7,8,9,10),將四叉樹單元(子域)分割成兩個次子域SS1 和SS2.在兩個次子域SS1 和SS2 內,分別采用式(25)計算這兩個次子域的剛度矩陣和,然后集成四叉樹單元(子域)的剛度矩陣,其形式如下

圖6 改進單元節點編號示意圖(*表示比例中心)Fig.6 Schematic diagram of node number of enriched element(* means scale center)

式中,kuu表示與真實節點(即角節點和懸掛節點)相應的剛度矩陣貢獻;kaa表示與虛節點(裂紋與四叉樹單元邊界的交叉點)相應的剛度矩陣貢獻;kua和kau為真實節點與虛節點的耦合貢獻.

如圖6(b)所示,四叉樹單元被裂紋面部分切割,取圍繞裂尖一圈的若干單元組成一個超級單元,超級單元的比例中心取為裂尖位置,該超級單元包括兩個虛節點(節點編號: 17,18),采用式(25)計算超級單元的剛度矩陣.

2.5 改進型比例邊界有限元法的控制方程

由2.4 節的討論可知,無論何種改進類型的單元,單元剛度矩陣均具有式(26)的形式,因此,可得到改進型SBFEM 的控制方程為

約78%的學生表示遇到過此類困難。隨機抽樣訪談中，部分受訪學生反映曾因文化差異在涉外場合導致過溝通障礙，也有為此鬧過笑話的經歷。

式中,u為常規節點(角節點和懸掛節點)的位移未知量,a為附加的虛節點的位移未知量.使用這種方法引入SBFEM 的附加節點未知量a可以解釋為裂紋與單元邊界交點處的實際位移.

3 裂紋生長準則

在SBFEM 中,I 型和II 型應力強度因子KI和KII可根據裂紋尖端 θ=0o處的奇異應力和求得,即

在包含裂紋尖端的超級單元中,當r→0 時,式(20)包含對角塊,奇異應力可表達為

圖7 為裂尖局部坐標系,通過插值可以得到裂紋尖端 θ=0o處的奇異應力模態值.奇異應力模態(η(θ))在沿裂紋子域邊界的高斯積分點使用式(30)計算. θ=0o時的應力模態可以通過插值得到.本研究中使用樣條插值來獲得更精確的應力.

圖7 裂尖局部坐標系Fig.7 Local coordinate system at the crack tip

根據關系式 ξ=r/L0,應力強度因子可通過下式計算[3]最大周向應力準則用于確定裂紋擴展方向,該準則給出了以下裂紋擴展方向

等效應力強度因子Ke可通過下式計算

4 數值算例

4.1 中心斜裂紋問題

本算例為復合型應力強度因子的求解問題,如圖8 所示,雙向拉伸板中心含一傾斜裂紋,裂紋傾角為 θ,裂紋半長度為a,該板在無限遠處受單軸拉伸載荷 σ 作用.該問題I 和II 型應力強度因子的解析解為[29]

圖8 含中心斜裂紋的拉伸板Fig.8 Tension plate with an inclined centre crack

數值計算時,取拉伸板邊長L=16 m,裂紋半長度a=0.5 m,拉伸載荷 σ=1.0 kPa,按平面應變問題考慮,為消除剛體位移,板的下邊緣豎向約束,板的左下角水平向約束,如圖9 所示,四叉樹網格單元總數為1276,節點總數為1365,彈性模量取為1.0 MPa,泊松比取為0.2,這里取的彈性模量和泊松比數值僅為了驗證論文建議方法計算復合型應力強度因子的正確性,彈性模量和泊松比的取值不同對計算結果影響很小.

圖9 四叉樹網格Fig.9 Quad-tree cells

圖10 給出了采用建議的改進型SBFEM 計算得到的裂紋不同傾斜角度下拉伸板的I 型和II 型應力強度因子值,并給出該問題的解析解作為比較,從圖中可看出,數值解與解析解吻合較好,從而驗證了數值方法的有效性.四叉樹網格可實現網格的快速過渡,在含有裂紋的局部區域可采用較密的網格.

圖10 拉伸板應力強度因子隨傾角變化曲線Fig.10 Stress intensity factors of the tensile plate with different crack inclination angles

4.2 含預制裂紋的重力壩模型

如圖11,某混凝土重力壩比例模型含1 條初始裂紋,試件的尺寸和邊界條件如圖所示.數值計算時,該重力壩模型被離散成四叉樹網格,包括4140 個單元和4392 個節點,按平面應變問題考慮.重力壩模型的材料參數取值參照文獻[30],彈性模量E=35.7 GPa、泊松比 ν=0.1、密度 ρ=2400 kg/m3、斷裂能G=0.184 N/mm、臨界應力強度因子KIC=81.456 N·mm1/2、裂紋生長長度取為150 mm.

重力壩模型所受載荷包括自重和水載荷.水載荷的模擬簡化成在圖11 中四個加載點處施加一定大小的集中力,施加的集中力大小比例如圖所示.初始裂紋的位置如圖11 所示,深度為150 mm.圖12給出了裂紋的最終擴展路徑,同時與文獻[30]的數值和試驗結果進行比較,二者一致性較好.

圖11 含初始裂紋的重力壩模型(單位: mm)Fig.11 Gravity dam model with initial crack (unit: mm)

圖12 重力壩模型的裂紋開裂擴展路徑Fig.12 Crack propagation path of gravity dam model

4.3 含雙裂紋雙孔洞的薄板模型

如圖13 所示,一尺寸為 20 mm×10 mm 帶有半徑r=2 mm 的雙孔洞薄板受雙向集中拉伸載荷作用,板的左、右邊緣處分別有長a=1 mm 的裂紋.數值計算時,該薄板模型被離散成四叉樹網格,包括2456 個單元和3000 個節點.按平面應力問題考慮,模型的材料參數取值參照文獻[13],彈性模量E=200 GPa、泊松比ν=0.3、密度ρ=7800 kg/m3、臨界應力強度因子KIC=1500 N·mm1/2、裂紋生長長度取為0.5 mm.

圖13 雙孔板結構內多裂紋擴展問題(單位: mm)Fig.13 Multiple crack propagation in a thin plate with two holes(unit: mm)

圖14 給出了預測的多裂紋擴展路徑,當開始施加外載荷后,兩條裂紋都向鄰近的孔洞位置擴展,隨后裂紋逐漸改變擴展角度,以幾乎水平的方向向薄板中部擴展,當兩裂紋都擴展到試件中部以后,由于裂尖的應力集中效應,兩條裂紋彼此相向擴展,模擬結果與文獻[13]采用的自適應網格的SBFEM 數值結果進行了比較,二者一致性較好.

圖14 含雙孔板問題的預測多裂紋開裂擴展路徑Fig.14 Predicted multiple crack paths for the plate with two holes

4.4 L 形混凝土面板

Winkler 等[31]對L 形混凝土面板進行了實驗測試,許多學者已經對其進行了模擬,以驗證數值模型.面板的幾何尺寸如圖15 所示.面板承受垂直集中載荷F.混凝土面板的細觀結構,使用具有代表性的直徑為45 mm、32.5 mm 和17.5 mm 的三段式圓形骨料級配,骨料面積分數為40%.L 形面板的細觀結構如圖16 所示.要生成四叉樹網格,具體圖像的大小為1024 × 1024 像素,最小邊緣長度設置為4 × 4 像素.L 形面板的材料特性見表1.初始裂紋生長長度取為5 mm.

表1 L 型板的材料參數Table 1 Material properties for the L-shaped panel

圖15 L 型板的示意圖(單位: mm)Fig.15 Schematic diagram of an L-shaped panel (unit: mm)

圖16 L 型板的細觀結構Fig.16 Mesostructure of L-shaped panel

圖17 給出了當前數值結果計算得到的裂紋開裂路徑并與試驗結果進行了對比,圖18 給出了數值計算得到的載荷-加載點位移曲線并與試驗結果進行了對比,從圖中可以看出,當前建議的基于四叉樹網格和改進型比例邊界有限元法的模擬結果均與試驗結果吻合較好.

圖17 開裂路徑Fig.17 Crack path

圖18 載荷-加載點位移曲線Fig.18 Load-displacement curves at loading point

5 結論

SBFEM 是一種半解析的數值方法,由于其本身的固有優勢,可以方便使用懸掛節點而無需任何特殊處理,且SBFEM 可以在求解過程中自動解析地求出裂尖的應力強度因子,SBFEM 與四叉樹可完美結合并實現網格的全自動剖分.論文提出的改進型SBFEM 繼承了傳統SBFEM 的所有優點,通過引入虛節點的思想,將裂紋與四叉樹單元邊界交叉點作為虛節點,虛節點的自由度作為附加自由度處理,并采用水平集函數表征裂紋面,在模擬裂紋擴展問題時無需進行裂尖的局部網格重剖分工作.

最后,通過三個數值算例驗證了建議方法的數值精度,提出的基于四叉樹網格的改進型SBFEM可準確求出裂紋的復合型應力強度因子;且能夠較好地模擬出裂紋的開裂擴展路徑,裂紋擴展過程無需進行網格重剖分.