孫立國 江守燕 杜成斌

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

引言

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

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

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

本工作在文獻(xiàn)[24-25]的研究基礎(chǔ)上,結(jié)合四叉樹網(wǎng)格的優(yōu)勢,提出基于圖像四叉樹的改進(jìn)型比例邊界有限元法,通過結(jié)構(gòu)域圖像實現(xiàn)結(jié)構(gòu)域網(wǎng)格的全自動剖分,且結(jié)構(gòu)內(nèi)不連續(xù)處局部范圍的網(wǎng)格粗細(xì)過渡方便、快速.

1 基于圖像的四叉樹網(wǎng)格

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

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

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

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

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

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

2 改進(jìn)型比例邊界有限元法

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

2.1 裂紋的水平集描述

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

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

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

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

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

2.2 改進(jìn)子域的識別

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

圖3 四叉樹網(wǎng)格與裂紋間的切割關(guān)系Fig.3 Relationship between the quadtree cell and a crack

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

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

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

步驟4: 對于求解域內(nèi)的所有四叉樹單元,重復(fù)步驟2、步驟3,即可識別出所有改進(jìn)單元.

2.3 改進(jìn)的比例邊界有限元子域

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

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

圖4 改進(jìn)的SBFE 子域Fig.4 Enriched SBFE subdomain

2.4 改進(jìn)子域剛度矩陣的求解

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

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

SBFEM 位移模式可以表示為

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

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

線性算子矩陣L可表示為

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

根據(jù)胡克定律,應(yīng)力場可表示為

式中,D為本構(gòu)矩陣且

根據(jù)虛功原理,可得SBFE 方程為

載荷向量為

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

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

式中,哈密頓矩陣Z為

通過Schur 分解,可得

其中,Schur 矩陣S和變換V可表達(dá)為

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

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

SBFE 子域的邊界節(jié)點力為

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

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

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

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

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

2.5 改進(jìn)型比例邊界有限元法的控制方程

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

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

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

3 裂紋生長準(zhǔn)則

在SBFEM 中,I 型和II 型應(yīng)力強(qiáng)度因子KI和KII可根據(jù)裂紋尖端 θ=0o處的奇異應(yīng)力和求得,即

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

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

圖7 裂尖局部坐標(biāo)系Fig.7 Local coordinate system at the crack tip

根據(jù)關(guān)系式 ξ=r/L0,應(yīng)力強(qiáng)度因子可通過下式計算[3]最大周向應(yīng)力準(zhǔn)則用于確定裂紋擴(kuò)展方向,該準(zhǔn)則給出了以下裂紋擴(kuò)展方向

等效應(yīng)力強(qiáng)度因子Ke可通過下式計算

4 數(shù)值算例

4.1 中心斜裂紋問題

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

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

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

圖9 四叉樹網(wǎng)格Fig.9 Quad-tree cells

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

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

4.2 含預(yù)制裂紋的重力壩模型

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

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

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

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

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

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

圖13 雙孔板結(jié)構(gòu)內(nèi)多裂紋擴(kuò)展問題(單位: mm)Fig.13 Multiple crack propagation in a thin plate with two holes(unit: mm)

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

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

4.4 L 形混凝土面板

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

表1 L 型板的材料參數(shù)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 型板的細(xì)觀結(jié)構(gòu)Fig.16 Mesostructure of L-shaped panel

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

圖17 開裂路徑Fig.17 Crack path

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

5 結(jié)論

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

最后,通過三個數(shù)值算例驗證了建議方法的數(shù)值精度,提出的基于四叉樹網(wǎng)格的改進(jìn)型SBFEM可準(zhǔn)確求出裂紋的復(fù)合型應(yīng)力強(qiáng)度因子;且能夠較好地模擬出裂紋的開裂擴(kuò)展路徑,裂紋擴(kuò)展過程無需進(jìn)行網(wǎng)格重剖分.