非連續(xù)問題中單元分割的模板方法1)

王理想 文龍飛 肖桂仲,?? 田榮,2)

?(中物院高性能數(shù)值模擬軟件中心,北京 100088)

?(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100088)

??(南京理工大學(xué),南京 210094)

引言

非連續(xù)問題廣泛存在于基礎(chǔ)研究[1-4]和工程應(yīng)用[5-8]中,如裂紋、孔洞、夾雜、多材料、流固耦合、相變、剪切帶等.研究這類問題,具有重要的科學(xué)意義和廣泛的應(yīng)用價值.

傳統(tǒng)有限元法模擬這類問題面臨諸多困難,如網(wǎng)格劃分需與非連續(xù)面對齊、非連續(xù)面更新(如裂紋擴展)時需網(wǎng)格重分重映、空間梯度變化較大處(如裂尖處)需進行網(wǎng)格加密等.

鑒于傳統(tǒng)有限元法局限性,1990 年代開始,發(fā)展出一批基于單位分解(PU)[9]的數(shù)值方法,如單位分解有限元(PUFEM)[10]、廣義有限元(GFEM)[11-15]、數(shù)值流形法(NMM)[16-20]、擴展有限元(XFEM)[21-25]等等.這類方法通過在基函數(shù)中引入先驗的加強函數(shù),可有效求解各類連續(xù)、非連續(xù)問題.

特別針對裂紋問題,1999 年美國西北大學(xué)Ted Belytschko 研究組提出XFEM[21-22],在學(xué)術(shù)界獲得廣泛關(guān)注和迅速發(fā)展.大型商業(yè)有限元軟件紛紛添加XFEM 模塊,標(biāo)志該方法在工業(yè)界獲得認(rèn)可.

XFEM 通過引入C?1型局部加強函數(shù),在函數(shù)空間逼近裂紋的非連續(xù)性.由此,裂紋幾何獨立于模擬網(wǎng)格,裂紋擴展時亦無需網(wǎng)格重分重映.通過引入解析性質(zhì)的加強函數(shù),其計算精度亦獲得大幅提高.

XFEM 在模擬裂紋擴展方面取得巨大成功,但同時也面臨兩大困擾:總體剛度矩陣高度病態(tài)和動力學(xué)計算時額外自由度上能量無法正確傳遞[23-24].前者導(dǎo)致迭代法求解收斂緩慢甚至不收斂;后者導(dǎo)致動力學(xué)求解實施困難.有鑒于此,作者團隊基于無額外自由度的裂尖插值格式,提出一種改進型XFEM(improved XFEM)[26-31],成功克服上述困難.

XFEM 類方法在求解非連續(xù)問題時,面臨非連續(xù)單元積分問題——此類單元受到非連續(xù)面的分割,若使用常規(guī)高斯積分,精度將嚴(yán)重?fù)p失,故需進行特殊處理.目前有以下幾類處理方法[23,32-33].

(1)高階高斯積分法[34-35]:在單元內(nèi)布置大量高斯積分點(如6×7[34],8×8[35]),判斷積分點與裂紋面的位置關(guān)系,從而對積分點所在積分域進行積分.該方法實現(xiàn)簡單,但精度較低[11].

(2)自適應(yīng)積分法[11,36-37]:將初始非連續(xù)單元進行網(wǎng)格細(xì)分,在細(xì)分后的格子上自適應(yīng)布置高斯積分點,然后進行非連續(xù)單元積分.該方法實現(xiàn)較簡單,精度獲得提高,但需要布置大量高斯積分點.

(3)子網(wǎng)格積分法[38-39]:通過子網(wǎng)格剖分,將非連續(xù)單元與裂紋進行幾何分割,非連續(xù)單元積分在子網(wǎng)格上實現(xiàn).由于該方法一般引入精確幾何操作,精度較高,但實現(xiàn)復(fù)雜、魯棒性差.

(4)矩量擬合(moment fitting)法[40-42]:通過一族函數(shù)(如冪函數(shù))在非連續(xù)單元上的精確積分,構(gòu)造一組非線性方程組,求解該方程組獲得積分點位置和權(quán)重,再使用所求得的積分點和權(quán)重進行非連續(xù)單元積分.該方法不需要對單元進行分割,僅需引入少量積分點就可取得較高精度,但由于每個單元都需要求解非線性方程組,導(dǎo)致求解效率降低.

(5)特殊積分法:這類方法的思想是將原積分進行一定特殊處理后再進行積分,具有實現(xiàn)簡便或精度高的優(yōu)點,但通用性一般.這類方法包括:①等效法.例如:Ventura[43]將Heaviside 函數(shù)等效為二次多項式,Abedian 和Düster[44]將非連續(xù)函數(shù)等效為Legendre 多項式;由此把非連續(xù)積分轉(zhuǎn)換為連續(xù)積分進行計算.②降維法.例如:Sudhakar 和Wall[45]通過高斯定理將體(面)積分轉(zhuǎn)換為面(線)積分.③變換法.例如:Mousavi 和Sukumar[46]使用廣義Duffy 變換,將在三角形(金字塔)上積分變換為在正方形(立方體)上積分.此外,還有其他特殊積分方法[47-48],不再一一贅述.

上述幾類方法在精度、效率、實現(xiàn)或通用性方面存在各自的優(yōu)點和缺點.本文從工程實用化角度出發(fā),提出一種基于水平集的非連續(xù)單元模板分割(單元積分)方法.該方法采取模板查詢的方式進行非連續(xù)單元子剖分,在三角形子網(wǎng)格上布置高斯積分點.相比于精確子剖分方法,該方法避免復(fù)雜幾何操作,可提高計算效率和魯棒性.將該方法與常規(guī)XFEM、改進型XFEM 進行結(jié)合,應(yīng)用于孔洞、夾雜、裂紋等非連續(xù)問題分析中,從而為XFEM 類方法在實際工程應(yīng)用中提供有效支撐.

1 控制方程

1.1 強形式

如圖1 所示非連續(xù)靜力問題,應(yīng)滿足平衡方程

式中,?為梯度算子,σ為應(yīng)力張量,b為體力.

圖1 二維非連續(xù)靜力問題描述Fig.1 Illustration of statics with discontinuities in 2D

位移邊界條件和應(yīng)力邊界條件分別如下所示

式中,u為位移,為在位移邊界Γu處給定的位移值,為在應(yīng)力邊界Γt處給定的應(yīng)力值,Γ0為無應(yīng)力邊界,Γc為裂紋面,Γh為孔洞面,n為各個面上單位外法向量.

對于線彈性問題,其本構(gòu)關(guān)系如下

式中,D為彈性張量,ε為應(yīng)變張量.

小變形下,應(yīng)變ε和位移u滿足如下幾何關(guān)系

1.2 弱形式

根據(jù)虛功原理,平衡方程(1)與應(yīng)力邊界條件(3)可轉(zhuǎn)化為如下弱積分形式

式中,δu為虛位移,δε為虛應(yīng)變.

將幾何關(guān)系(5)代入平衡方程弱積分式(6),可得

2 非連續(xù)界面水平集描述

本文使用水平集法[49]隱式描述非連續(xù)界面,如孔洞面、夾雜面、裂紋面、多材料界面等.非連續(xù)界面可分為強、弱兩種:裂紋面屬于強非連續(xù)界面,夾雜面和多材料界面屬于弱非連續(xù)界面.根據(jù)不同類型加強形式,孔洞面可以是強非連續(xù)界面,也可以是弱非連續(xù)界面.本文不對其進行展開討論.

2.1 非連續(xù)面水平集函數(shù)φ(x)

非連續(xù)面水平集函數(shù)φ(x)定義為求解域? 內(nèi)任意給定一點x∈? 到非連續(xù)面Γd的含符號距離,其數(shù)學(xué)表達式為

式中,||·||為歐式范數(shù),n為非連續(xù)面單位外法向量,sign(·)為符號函數(shù).上式中的非連續(xù)面Γd可以是裂紋面Γc、夾雜面Γi、孔洞面Γh或多材料界面Γm.

2.2 裂紋尖端水平集函數(shù)ψ(x)

裂紋尖端水平集函數(shù)ψ(x)定義為求解域? 內(nèi)任意給定一點x∈? 到裂紋尖端xt的含符號距離,其數(shù)學(xué)表達式為

式中,t為沿裂紋擴展方向的單位向量.針對裂紋所構(gòu)造的水平集函數(shù),如圖2 所示.

圖2 裂紋水平集函數(shù)構(gòu)造(根據(jù)文獻[33]重畫)Fig.2 Construction of level set functions for a crack(redrawn after Ref.[33])

2.3 裂紋尖端極坐標(biāo)系(r,θ)

裂紋可通過水平集函數(shù)φ(x)和ψ(x)進行描述.除此之外,還需要定義裂紋尖端極坐標(biāo)系(r,θ),其數(shù)學(xué)表達式可通過水平集函數(shù)給出

3 擴展有限元

3.1 常規(guī)擴展有限元:弱非連續(xù)問題

對于弱非連續(xù)問題,基于常規(guī)XFEM 的位移場逼近可表示為

3.2 常規(guī)擴展有限元:強非連續(xù)問題

對于裂紋問題,除了考慮非連續(xù)面(裂紋面)階躍加強之外,還需考慮裂尖奇異加強.此時,基于常規(guī)XFEM 的位移場逼近為

式中,I 為全部節(jié)點集合,J 為裂紋面加強節(jié)點集合,K 為裂尖加強節(jié)點集合;i,j,k,α 為序數(shù);為常規(guī)自由度,為裂紋面加強額外自由度,為裂尖加強額外自由度;N(x)為標(biāo)準(zhǔn)有限元形函數(shù),Fjump(x)為階躍加強函數(shù),Ftip(x)為裂尖奇異加強函數(shù).Fjump(x)常使用如下Heaviside 函數(shù)

對于各向同性線彈性材料,Ftip(x)可取為

3.3 改進型擴展有限元:裂紋問題

在使用改進型擴展有限元(IXFEM)求解時,僅考慮裂紋問題,其位移場逼近形式為[26-31]

式中,I 為所有節(jié)點集合,J 為裂紋面加強節(jié)點集合,K 為裂尖加強節(jié)點集合,Kk?K 為節(jié)點k處裂尖加強影響域內(nèi)節(jié)點集合;i,j,k,m為序數(shù);為常規(guī)自由度,為裂紋面加強額外自由度;N(x)為標(biāo)準(zhǔn)有限元形函數(shù),Fjump(x)為階躍加強函數(shù),與式(15)定義相同,(x)為裂尖局部加強函數(shù)

其中,p(x)為基向量,pk為p(xk)的縮寫;A為矩量矩陣,為A?1的第一列,為的第一個元素;δ為Kronecker 符號.p(x)和A的表達式如下所示

式中,c為常數(shù),一般取為2 ～3;h為網(wǎng)格尺寸;Ftip(x)為裂尖奇異加強函數(shù),其定義與式(16)相同.

觀察式(17)與式(14)可看出,改進型XFEM 與常規(guī)XFEM 的區(qū)別僅在于裂尖加強項.改進型XFEM裂尖加強項中不再包含額外自由度,可有效解決常規(guī)XFEM 線性相關(guān)問題.

3.4 混合單元修正

混合單元為同時包含有限元節(jié)點與加強節(jié)點的單元(如圖3 所示).為保證XFEM 算法精度和收斂性,混合單元形函數(shù)需加以修正[51].

圖3 混合單元處理:處理之前(左);處理之后(右)Fig.3 Blending element treatment:Before(left);after(right)

引入如下斜坡函數(shù)[52]

則在混合單元?bld上,位移逼近具有如下加權(quán)形式

4 控制方程離散

通過Galerkin 法,并應(yīng)用位移逼近式(12)或式(14)或式(17),可將平衡方程弱積分式(7)離散為如下線性方程組形式

式中,K為剛度矩陣,U為自由度向量,F為載荷向量.對于不同位移逼近,K,U,F形式略有不同.

若采用式(12)進行離散,式(23)具有如下形式

若采用式(14)進行離散,式(23)具有如下形式

若采用式(17)進行離散,式(23)具有如下形式

式(24)～式(26)中K和F可統(tǒng)一表示為

其中,α,β=u,a或u,b,c或u,b,表示所采用的位移逼近方式;B為應(yīng)變矩陣,N為形函數(shù)矩陣.

5 非連續(xù)單元分割方法

本方法的核心思想是,根據(jù)水平集在單元內(nèi)的近似屬性,假定裂紋在單元內(nèi)平直分布,由此根據(jù)水平集值特定形態(tài)對非連續(xù)單元進行分割.該方法由于采用模板查詢的方式,相比于精確分割,可避免復(fù)雜幾何操作,從而提高計算效率和魯棒性.

所分割的非連續(xù)單元分為兩類:一類是被非連續(xù)面完全貫穿的單元,稱為“切割單元”;一類是被非連續(xù)面(即裂紋面)部分貫穿的單元,稱為“裂尖單元”.下面分別給出這兩類單元的分割子剖分.

5.1 切割單元

5.1.1 三角形單元

首先,遍歷三角形切割單元水平集值所有形態(tài)并建立標(biāo)準(zhǔn)單元分割模板庫,如表1 和圖4 所示.

表1 三角形切割單元水平集值符號形態(tài)Table 1 LSV sign patterns of triangular cut elements

圖4 三角形切割單元分割形態(tài)Fig.4 The partitioning patterns of triangular cut elements

表1 中,case 表示所有水平集值符號的可能情況,共33=27 種;status 表示單元水平集值符號個數(shù)所組成的三元組,每個三元組表示一種狀態(tài),共10種;type 表示切割單元類型,共3 種.每種類型可分為(a)、(b)兩種亞型,由sign 進行標(biāo)識.三角形切割單元的3 種切割類型及其亞型,如圖4 所示.

然后,根據(jù)三角形切割單元水平集值,對非標(biāo)準(zhǔn)單元進行形態(tài)查詢(查詢type 和sign)和模板插值.對棱ij上交點p的位置進行線性插值

其中,λ 可由水平集值插值得到

交點p的坐標(biāo)可由i和j的坐標(biāo)和水平集值給出

最后,套用標(biāo)準(zhǔn)單元分割模板實現(xiàn)單元分割和子三角剖分,如算法1 所述.序號ijk的具體排列根據(jù)不同切割單元類型預(yù)先排定,可通過查表得到.

5.1.2 四邊形單元

首先,遍歷四邊形切割單元水平集值所有形態(tài)并建立標(biāo)準(zhǔn)單元分割模板庫.四邊形切割單元中,共34=81 種情況、15 種狀態(tài).

在不考慮一些特殊情況(如裂紋在單元內(nèi)發(fā)生較大轉(zhuǎn)折、兩條裂紋同時切割一個單元等)的假設(shè)下,四邊形切割單元可分為5 種類型,每個類型分為(a)、(b)兩種亞型,如圖5 所示.

圖5 四邊形切割單元分割形態(tài)Fig.5 Partitioning patterns of quadrilateral cut elements

然后,根據(jù)四邊形切割單元水平集值,對非標(biāo)準(zhǔn)單元進行形態(tài)查詢(查詢type 和sign)和模板插值,其中模板插值與三角形切割單元中的實現(xiàn)一致.

最后,套用標(biāo)準(zhǔn)單元分割模板實現(xiàn)單元分割和子三角剖分,如算法2 所述.與三角形單元類似,四邊形單元中的序號ijkl具體排列也已根據(jù)不同切割單元類型預(yù)先排定,可通過查表得到.

5.2 裂尖單元

僅考慮裂尖在單元內(nèi)的情況.裂尖單元分割形態(tài)如圖6 所示,分別對其進行子剖分,詳見算法3.

值得注意的是,本文方法在實現(xiàn)過程中:(1)預(yù)先將若干標(biāo)準(zhǔn)單元進行子三角化,形成標(biāo)準(zhǔn)單元分割模板庫,并在計算機程序中進行存儲;這一存儲過程發(fā)生在程序編譯期,不占用計算時間.(2)在處理任意單元分割時,僅需查詢存儲在計算機中的標(biāo)準(zhǔn)單元分割模板庫進行匹配,就可對任意單元進行分割;這樣可避免復(fù)雜的幾何操作,從而大大提高程序計算效率和魯棒性.

圖6 裂尖單元分割形態(tài)Fig.6 The partitioning patterns of tip elements

5.3 非連續(xù)單元積分

為便于區(qū)分,以下各式中,以x表示全局坐標(biāo),ξ表示有限元母單元上的自然坐標(biāo),表示子三角形母單元上的局部坐標(biāo).

切割單元、裂尖單元非連續(xù)積分可分別表示為

式(32)和式(33)中,Nsub表示子剖分三角形個數(shù),Ngp表示子三角形上高斯點個數(shù),二者上標(biāo)“+”、“–”、“t”分別表示裂紋面上、下和裂尖材料域;|J|為全局坐標(biāo)與自然坐標(biāo)之間雅可比矩陣行列式;ξkl為高斯點自然坐標(biāo),Wkl為高斯點權(quán)重,分別可表示為

6 裂紋相關(guān)計算

6.1 應(yīng)力強度因子

相互作用積分法計算應(yīng)力強度因子(SIF),具有精度高、適應(yīng)范圍廣的優(yōu)點.故本文采用該方法進行計算,其定義式如下[53]

相互作用積分與真實場和附加場應(yīng)力強度因子之間存在以下關(guān)系

式中,E′為等效楊氏模量.對于平面應(yīng)變問題:E′=E/(1 ?ν2);對于平面應(yīng)力問題:E′=E.

在實際計算中,常在圍線Γ0外再增設(shè)一回路C(如圖7),將積分項乘以一光滑權(quán)函數(shù)q(x),并將線積分轉(zhuǎn)化為等效面積分[54]

其中,q(x)在Γ 內(nèi)取1,在C0外取0.A為由Γ,C0,C+和C?圍成的閉合區(qū)域;mj為A的單位外法向量.

圖7 相互作用積分域定義Fig.7 Interaction integral domain

6.2 裂紋擴展

本文采用最大周向拉應(yīng)力強度因子理論[55]計算裂紋擴展方向.令以下臨界斷裂韌度對θ 偏導(dǎo)為

可得裂紋擴展角度

其中,θ0∈(?π,π).當(dāng)KII=0 時,取θ0=0.

7 算例分析

7.1 孔洞問題

如圖8(a)所示,在一邊長l=2 m 的方板中心有一半徑為a=0.2 m 的圓孔.板受x向均勻拉伸作用,拉應(yīng)力σ0=1 Pa.板的彈性模量E=1000 Pa,泊松比ν=0.3.計算網(wǎng)格如圖8(b)所示,網(wǎng)格數(shù)為40×40.假設(shè)該問題為平面應(yīng)變問題.在極坐標(biāo)系下,該問題應(yīng)力場解析解可表示為[33]

圖8 方板中心帶一圓孔:(a)問題描述;(b)計算網(wǎng)格Fig.8 A square plate with a circular hole at its center:(a)Problem definition;(b)simulation mesh

圖9 σx 應(yīng)力云圖對比:本文(左側(cè));解析解(右側(cè))Fig.9 σx stress contour:Present(left);exact(right)

使用XFEM 聯(lián)合本文非連續(xù)單元分割方法計算得到的應(yīng)力云圖,與解析解對比如圖9 所示.從該圖可看出,本文與解析解計算結(jié)果非常接近.特別是在孔洞周圍,本文方法計算得到的應(yīng)力集中因子為2.992,與理論值3.0 非常接近.進一步地,取y軸上a≤y≤L/2 區(qū)間范圍內(nèi)的應(yīng)力,對比本文方法與解析解計算結(jié)果,如圖10 所示.該圖定量說明本文計算結(jié)果與解析解符合較好.該孔洞算例表明,本文發(fā)展的非連續(xù)單元分割方法在XFEM 求解孔洞問題上具有有效性.

圖10 本文方法計算結(jié)果與解析解對比Fig.10 Comparison between present and exact solutions

為研究本文分割方法計算效率,在本算例中,同時使用本文發(fā)展的非連續(xù)單元分割方法以及德勞內(nèi)三角剖分(Delaunay triangulation)方法,進行非連續(xù)單元分割.以下以Sub-DT 表示德勞內(nèi)三角剖分方法;以Sub-LSV 表示本文非連續(xù)單元分割方法.

由于計算模型的對稱性,本文僅給出1/4 孔洞與固體單元(共7 個)分割計算時間對比,如表2 所示.從該表可知,使用德勞內(nèi)方法(Sub-DT)計算總時間為219μs,而使用本文方法(Sub-LSV)計算總時間僅為57 μs.該算例說明,本文非連續(xù)單元分割方法在計算效率上明顯優(yōu)于傳統(tǒng)德勞內(nèi)方法.

表2 使用不同分割方法所消耗計算時間Table 2 Computational costs with different partitioning methods

7.2 夾雜問題

如圖11(a)所示,在一邊長l=2 m 的方板中心有一半徑為a=0.2 m 的圓孔.板受y向均勻壓縮作用,壓應(yīng)力σ0=1 Pa.基質(zhì)材料彈性模量E1=1000 Pa,泊松比ν1=0.3;夾雜材料彈性模量E2=3000 Pa,泊松比ν2=0.3.假設(shè)該問題為平面應(yīng)變問題.

圖11 方板中心帶一夾雜:(a)問題描述;(b)有限元網(wǎng)格Fig.11 A square plate with a circular inclusion at its center:(a)Problem definition;(b)finite element mesh

XFEM 計算網(wǎng)格仍采用圖8(b)中的40×40 規(guī)則網(wǎng)格.該問題無解析解,為驗證方法可行性,同時使用有限元法(FEM)進行計算,網(wǎng)格如圖11(b).該網(wǎng)格具有1685 個節(jié)點,3288 個三角形單元.

使用XFEM 聯(lián)合本文非連續(xù)單元分割方法計算得到uy位移云圖,與FEM 對比,如圖12 所示.XFEM與FEM 計算的uy最大值均為1.809 mm,二者一致.取直線y=x上?L/2≤x≤L/2 區(qū)間范圍內(nèi)ux和uy位移,對比二者結(jié)果,如圖13 所示,圖中曲線兩兩重合.進一步分析可知,二者總位移平均偏差為0.04%.該夾雜算例表明,本文發(fā)展的非連續(xù)單元分割方法在XFEM 求解夾雜問題上同樣具有有效性.

圖12 uy 位移云圖對比:本文(左側(cè));有限元(右側(cè))Fig.12 uy displacement contour:Present(left);FEM(right)

7.3 穩(wěn)定裂紋問題

如圖14(a)所示,一邊裂紋板受x向剪切作用,剪應(yīng)力τ0=1 Pa.板寬b=7 m、高2h=16 m,裂紋長a=3.5 m.板彈性模量E=1000 Pa,泊松比ν=0.3.計算網(wǎng)格如圖14(b)所示,網(wǎng)格數(shù)為19×39.假設(shè)為平面應(yīng)變問題.該問題中裂尖處應(yīng)力強度因子的參考解為[56]

在本算例中,同時使用本文發(fā)展的非連續(xù)單元分割方法(Sub-LSV)以及德勞內(nèi)三角剖分方法(Sub-DT),進行非連續(xù)單元分割.

圖13 本文方法計算結(jié)果與有限元對比Fig.13 Comparison between present and FEM solutions

圖14 邊裂紋板受單向剪切:(a)問題描述;(b)計算網(wǎng)格Fig.14 Plate with an edge crack under unidirectional shear:(a)Problem definition;(b)simulation mesh

分別聯(lián)合Sub-LSV,Sub-DT 和XFEM,IXFEM 計算裂尖處應(yīng)力強度因子,如表3 所示.一方面,通過對比Sub-LSV 和Sub-DT 兩種分割方法可知,不論是應(yīng)用于XFEM 還是應(yīng)用于IXFEM,Sub-LSV 均可取得接近于Sub-DT 的計算精度.Sub-LSV 精度略低的原因是:其分割裂尖單元為5 個子三角形,少于Sub-DT分割的6 個子三角形,導(dǎo)致精度略有降低.另一方面,通過對比XFEM 和IXFEM 可知,IXFEM 具有明顯更高的計算精度.該穩(wěn)定裂紋算例表明,本文發(fā)展的非連續(xù)單元分割方法在XFEM 和IXFEM 求解穩(wěn)態(tài)裂紋問題上均具有有效性.

表3 使用不同方法計算得到的應(yīng)力強度因子Table 3 Stress intensity factors from different approaches

7.4 擴展裂紋問題

如圖15(a)所示,為一帶邊裂紋雙臂梁模型,長l=6 m,寬w=2 m,裂紋長a=2.05 m.雙臂梁彈性模量E=1000 Pa,泊松比ν=0.3.梁右端固支,左端受非對稱集中力載荷P1=1 N,P2=1.01 N.計算網(wǎng)格如圖15(b)所示,網(wǎng)格數(shù)為25×75.假設(shè)裂紋擴展18 步,每一步擴展長度為0.1 m.

圖15 雙臂梁帶一邊裂紋:(a)問題描述;(b)計算網(wǎng)格Fig.15 Double cantilever beam with an edge crack:(a)Problem definition;(b)simulation mesh

分別使用Sub-LSV,Sub-DT 聯(lián)合IXFEM 計算裂紋擴展路徑,如圖16 所示.從該圖可看出,使用IXFEM+Sub-LSV 計算得到的裂紋擴展路徑,與使用IXFEM+Sub-DT 計算得到的裂紋擴展路徑具有高度重合性.該算例進一步說明,對于擴展裂紋問題,本文發(fā)展的非連續(xù)單元分割方法同樣具有有效性.

圖16 裂紋擴展路徑對比Fig.16 Comparison of crack propagation paths

7.5 三維裂紋問題

本文發(fā)展的非連續(xù)單元模板分割方法,可推廣至三維情形.但由于三維情況下,彎曲裂紋切割單元需要進行等參變換等諸多算法研究,因此另文敘述.此處,預(yù)先給出一三維邊裂紋算例,以供參考.

如圖17(a)所示,為一三維貫穿型邊裂紋板,該板寬w=7 m,高h(yuǎn)=16 m,厚t=4 m,裂紋長a=3.5 m.板彈性模量E=1000 Pa,泊松比ν=0.0.板底端固支,頂端受y軸正向拉伸載荷σ0=1 Pa 作用.計算網(wǎng)格如圖17(b)所示,網(wǎng)格數(shù)為19×39×10.

分別使用高階高斯積分法和本文切割方法,并聯(lián)合XFEM 進行求解.計算所得到的uy位移場對比,如圖18 所示.使用XFEM+高階高斯積分法,位移最大值為73.92 mm,最小值為?4.978 mm.使用XFEM+本文分割方法,位移最大值為73.91 mm,最小值為?4.976 mm.兩種方法的位移最大值偏差為0.01%,位移最小值偏差為0.04%.該算例說明,對于三維裂紋問題,本文發(fā)展的非連續(xù)單元分割方法同樣具有有效性.

圖17 三維貫穿型邊裂紋模型:(a)問題描述;(b)計算網(wǎng)格Fig.17 Three-dimensional cut-through edge crack model:(a)Problem definition;(b)simulation mesh

圖18 uy 位移場對比:(a)高階高斯積分法;(b)本文分割方法Fig.18 Comparison of uy displacement fields:(a)High-order Gauss integration method;(b)presented partitioning method

8 結(jié)論

本文提出一種基于單元水平集的模板分割方法,用于非連續(xù)單元子剖分和數(shù)值積分.首先,遍歷單元水平集值所有形態(tài)并建立標(biāo)準(zhǔn)單元分割模板庫;然后,根據(jù)單元水平集值,對非標(biāo)準(zhǔn)單元進行形態(tài)查詢和模板插值;最后,套用標(biāo)準(zhǔn)單元分割模板實現(xiàn)單元高效分割和子剖分.

將該方法與常規(guī)XFEM、改進型XFEM 進行結(jié)合,從而應(yīng)用于孔洞、夾雜、裂紋等非連續(xù)問題分析中.算例分析表明,本文提出的分割方法具有較高的計算精度.由于不引入復(fù)雜幾何操作,該模板分割方法同時具有較高計算效率和魯棒性,可為XFEM 類方法在非連續(xù)問題中提供有效支撐.

作為方法預(yù)研,本文僅給出二維非連續(xù)單元分割算法,用以探究模板分割可行性.事實上,筆者已初步將該算法推廣至三維情形.在三維情況下,對于彎曲裂紋切割單元需要進行等參變換等諸多算法研究,所以將另文敘述.本文分割方法同時還可用于界面與計算網(wǎng)格不對齊的其他數(shù)值方法,如數(shù)值流形法(NMM)、有限格子法(FCM)等等.