Hansbo類有限單元法的非連續(xù)分片試驗(yàn)

凌道盛，石吉森，張如如，王云崗

（1.浙江大學(xué) 軟弱土與環(huán)境土工教育部重點(diǎn)實(shí)驗(yàn)室，浙江 杭州310058；2.浙江大學(xué) 巖土工程研究所，浙江 杭州310058；3.浙江大學(xué) 城市學(xué)院，浙江 杭州310015）

裂紋（剪切帶）等非連續(xù)變形在土木工程、機(jī)械工程、航空航天等工程領(lǐng)域結(jié)構(gòu)破壞過(guò)程中普遍存在，非連續(xù)變形分析已成為相關(guān)工程領(lǐng)域的研究熱點(diǎn)之一.由于非連續(xù)變形的產(chǎn)生與擴(kuò)展是一個(gè)路徑事先未知的動(dòng)態(tài)過(guò)程，將裂紋限制在單元邊界的常規(guī)有限單元法存在明顯的局限性［1-2］，能夠表征單元內(nèi)部非連續(xù)變形的新型分析方法得到越來(lái)越多的關(guān)注.

為克服常規(guī)有限元法的困難，石根華［3］最早（1991年）基于單位分解（Partition of unity）和數(shù)值流形的思想，提出了適用于非連續(xù)變形分析的數(shù)值流形法［4-6］，雖然單位分解法的概念直到1997年才被Babu?ka等［7］正式提出.同樣基于單位分解思想，Belytschko等［8-9］通過(guò)引入Heaviside函數(shù)和裂尖函數(shù)描述裂紋變形的非連續(xù)性和裂尖應(yīng)力的奇異性，提 出 了 擴(kuò) 展 有 限 元 法（extended finite element method）［10-11］.Hansbo等［12-13］提出了基于常規(guī)有限元位移插值描述單元內(nèi)非連續(xù)變形的方法，被廣泛稱為Hansbo和Hansbo方法（Hansbo ＆ Hansbo's Method）.已經(jīng)證明［14］，Hansbo和Hansbo方法與未引入裂尖奇異函數(shù)的擴(kuò)展有限元法是完全等價(jià)的.Song 等［15］提 出 的 虛 節(jié) 點(diǎn) 法（phantom node method）實(shí)質(zhì)上是Hansbo和Hansbo方法的另一種形式.Ling等［16-18］提出了基于網(wǎng)格分離和位移映射的有限元位移插值技術(shù)及增強(qiáng)有限單元法（augmented finite element method）.與擴(kuò)展有限元法和Hansbo等提出的方法不同，增強(qiáng)有限單元法將開(kāi)裂后的單元視為2個(gè)完全獨(dú)立的單元，使得復(fù)雜裂紋擴(kuò)展模擬更加容易實(shí)現(xiàn).需要指出的是，Hansbo和Hansbo方法是增強(qiáng)有限單元法的一個(gè)特例，也是裂紋擴(kuò)展問(wèn)題中最常用的一種形式.目前，所有與Hansbo和Hansbo方法本質(zhì)相同的方法被統(tǒng)稱為Hansbo和Hansbo 類方法（Hansbo ＆ Hansbo's type of methods）.

為避免體積不可壓縮條件下模擬裂紋擴(kuò)展時(shí)出現(xiàn)網(wǎng)格自鎖，Dolbow 等［19］在擴(kuò)展有限元法位移插值中引入加強(qiáng)假定應(yīng)變（enhanced assumed strain），并設(shè)計(jì)了一個(gè)非連續(xù)分片試驗(yàn)（discontinuous path test）來(lái)測(cè)試其構(gòu)造的單元的性能.

本文基于作者提出的增強(qiáng)有限單元法，利用Dolbow 等設(shè)計(jì)的非連續(xù)分片試驗(yàn)對(duì)Hansbo 和Hansbo類有限單元法進(jìn)行測(cè)試，發(fā)現(xiàn)該類方法普遍不能以機(jī)器精度通過(guò)非連續(xù)分片試驗(yàn)，從理論上揭示了分片試驗(yàn)失敗的原因，并提出了相應(yīng)的解決方法.

1 增強(qiáng)有限單元法［16-18］

增強(qiáng)有限元法將常規(guī)有限單元分離為幾何上相互獨(dú)立的數(shù)學(xué)單元和物理單元，利用數(shù)學(xué)單元構(gòu)造離散位移函數(shù)，利用物理單元定義真實(shí)物理區(qū)域.與數(shù)值流形法和廣義有限單元法不同，增強(qiáng)有限單元法的關(guān)鍵在于引入映射法則將數(shù)學(xué)單元和物理單元的離散位移場(chǎng)有機(jī)關(guān)聯(lián)在一起，從而極大提高了有限單元法的適用性和靈活性，如圖1所示.最簡(jiǎn)單的一種映射法則是物理單元采用與數(shù)學(xué)單元完全相同的離散位移函數(shù)，即簡(jiǎn)單覆蓋法則.簡(jiǎn)單覆蓋法則也是數(shù)值流形法和廣義有限單元法在構(gòu)造物理單元（或積分網(wǎng)格）離散位移函數(shù)時(shí)采用的法則.為簡(jiǎn)單起見(jiàn)，本文也限于簡(jiǎn)單覆蓋法則.

圖1 增強(qiáng)有限單元Fig.1 Augmented finite element

利用分離的數(shù)學(xué)單元和物理單元可以實(shí)現(xiàn)任意強(qiáng)非連續(xù)變形模擬.考慮如圖2（a）所示被裂紋分割的單元，為簡(jiǎn)便起見(jiàn)，假定裂紋穿過(guò)前數(shù)學(xué)單元和物理單元完全重合.裂紋穿過(guò)后，物理單元被分割為2個(gè)不同的物理區(qū)域，如圖2（b）和（c）中陰影部分所示.為描述2個(gè)物理區(qū)域間的非連續(xù)變形，分別為每個(gè)物理區(qū)域構(gòu)造數(shù)學(xué)單元，并令2個(gè)新的數(shù)學(xué)單元（圖2中的1I2I3I4I和1II2II3II4II）具有和原數(shù)學(xué)單元（圖2中的1234）相同的幾何形狀，占據(jù)相同的空間區(qū)域.為保證與鄰單元間的位移連續(xù)性，數(shù)學(xué)節(jié)點(diǎn)1II、2II、3I、4I為原數(shù)學(xué)節(jié)點(diǎn)1、2、3、4，其余數(shù)學(xué)節(jié)點(diǎn)為新富集的數(shù)學(xué)節(jié)點(diǎn).不難看出，描述非連續(xù)變形時(shí)，基于簡(jiǎn)單覆蓋法則的增強(qiáng)有限單元法與Hans-bo和Hansbo方法本質(zhì)上是相同，不同的是，Hansbo和Hansbo方法將開(kāi)裂后的單元視為一個(gè)含裂紋的單元，而增強(qiáng)有限單元法將其當(dāng)作2個(gè)完全獨(dú)立的單元.

圖2 非連續(xù)變形的增強(qiáng)有限元描述Fig.2 Simulation of discontinuous deformation in AFEM

數(shù)學(xué)單元和物理單元的分離使得物理單元可以具有任意的幾何形狀，為實(shí)現(xiàn)剛度矩陣等在單元上的數(shù)值積分，增強(qiáng)有限單元法［18］采用與擴(kuò)展有限元法相同的求積技術(shù)，將物理區(qū)域細(xì)化為若干三角形子域，對(duì)每個(gè)子域采用Hammer積分.對(duì)于四邊形或可細(xì)化為四邊形的物理單元，也可采用Gauss積分.為簡(jiǎn)潔起見(jiàn)，本文用Q4和T3分別表示4節(jié)點(diǎn)四邊形數(shù)學(xué)單元和3 節(jié)點(diǎn)三角形數(shù)學(xué)單元；用H3和H1 分別表示三角形物理子域3 點(diǎn)和1 點(diǎn)Hammer積分；用G2和G1分別表示四邊形物理子域2×2和1×1點(diǎn)Gauss積分.于是，Q4H3代表含裂紋單元為4節(jié)點(diǎn)四邊形數(shù)學(xué)單元，積分采用三角形物 理 子 域3 點(diǎn)Hammer 積 分，如 圖3（a）所 示；T3G2H1則代表含裂紋單元為3節(jié)點(diǎn)三角形數(shù)學(xué)單元，積分混合采用G2和H1這2種方法，且G2優(yōu)先，如圖3（c）所示.圖3（b）還給出了Q4G2示意圖.

圖3 網(wǎng)格類型和積分方法Fig.3 Different meshes and integration methods

2 非連續(xù)分片試驗(yàn)

Dolbow 等［19］設(shè)計(jì)的非連續(xù)分片試驗(yàn)如圖4示，分析域Ω 為3×3的正方形區(qū)域，界面無(wú)作用力的裂紋Γ 將分析域分成2個(gè)相同的區(qū)域Ω＋和Ω-，兩側(cè)分別作用大小為2t和t MPa·m 的均布拉力.與文獻(xiàn)［19］不同，本文不考慮材料不可壓縮這種特殊情況，取材料彈性模量E＝1 000MPa，泊松比ν＝0.25.該問(wèn)題滿足平面應(yīng)變假定，其應(yīng)力精確解為

式中：σx為x 方向正應(yīng)力，σy為y 向正應(yīng)力，τxy為剪應(yīng)力，t為常規(guī)變量，可以為任意數(shù)值.

Dolbow 等［19］利用沒(méi)有富集裂尖奇異函數(shù)的擴(kuò)展有限元法對(duì)該問(wèn)題進(jìn)行分析，采用圖4 所示的3×3的網(wǎng)格，其中中間一行單元為含裂紋單元.本文以增強(qiáng)有限單元法為例，利用該算例對(duì)Hansbo和Hansbo類方法進(jìn)行分片試驗(yàn)測(cè)試.需要指出的是，對(duì)于圖4示網(wǎng)格，增強(qiáng)有限單元法采用12個(gè)單元，其中裂紋附近的6個(gè)單元的物理單元不再與對(duì)應(yīng)的數(shù)學(xué)單元重合.不難看出，Ω＋和Ω-對(duì)于增強(qiáng)有限單元法而言是2個(gè)相互獨(dú)立的分析域，可以分別進(jìn)行計(jì)算.不失一般性，本文僅對(duì)Ω-進(jìn)行分析.

圖4 非連續(xù)分片試驗(yàn)Fig.4 Discontinuous patch test

為對(duì)比研究數(shù)學(xué)單元類型、剛度矩陣積分方法對(duì)分片試驗(yàn)結(jié)果的影響，本文以圖4網(wǎng)格為基礎(chǔ)，構(gòu)造如圖5所示的10種網(wǎng)格.圖中粗實(shí)線代表數(shù)學(xué)網(wǎng)格，陰影部分代表物理區(qū)域，虛線代表物理單元的子域劃分.其中，（a）和（b）采用常規(guī)四邊形有限單元；（c）和（d）采用與圖4相同的數(shù)學(xué)網(wǎng)格；（g）和（h）采用梯形數(shù)學(xué)單元，且梯形兩平行邊與裂紋面平行；（i）和（j）中，裂紋面所在的數(shù)學(xué)單元為平行四邊形單元.

為便于評(píng)價(jià)計(jì)算精度，利用積分點(diǎn)應(yīng)力定義相對(duì)誤差指標(biāo)為

式中：σ為3個(gè)應(yīng)力分量中的任意一個(gè)；σexact和σh分別為應(yīng)力精確值和計(jì)算值；為特征應(yīng)力值，這里取為所有積分點(diǎn)中求最大值.10種網(wǎng)格計(jì)算結(jié)果如表1所示，其中試驗(yàn)結(jié)果以計(jì)算誤差是否達(dá)到機(jī)器誤差（10-15）為評(píng)價(jià)標(biāo)準(zhǔn).

圖5 分片試驗(yàn)網(wǎng)格Fig.5 Different meshes for discontinuous patch tests

表1 分片試驗(yàn)應(yīng)力相對(duì)誤差指標(biāo)統(tǒng)計(jì)表Tab.1 Relative error indexes of stresses for different meshes

由表1可以看出，采用三角形數(shù)學(xué)單元的網(wǎng)格能以機(jī)器精度通過(guò)分片試驗(yàn)；除Q4G2-S、Q4G2-M、Q4G2-P和Q4H3-P 這4 種情況外，采用四邊形數(shù)學(xué)單元的網(wǎng)格都不能以機(jī)器精度通過(guò)分片測(cè)試.

為進(jìn)一步考察數(shù)學(xué)單元和物理單元形狀對(duì)計(jì)算精度的影響，采用如圖6所示的Q4G2和Q4H3網(wǎng)格進(jìn)行測(cè)試.圖6中f 為一無(wú)量綱的量，4個(gè)單元的形狀隨著f 的變化而改變，當(dāng)f＝1時(shí)，即為圖5（c）和（d）示網(wǎng)格.

圖6 網(wǎng)格變化示意圖Fig.6 Mesh varying with f

如表2所給出了不同f 值時(shí)應(yīng)力相對(duì)誤差指標(biāo)計(jì)算結(jié)果.由表可見(jiàn)，單元形狀對(duì)計(jì)算誤差有顯著的影響.當(dāng)f＝0.5時(shí)，數(shù)學(xué)和物理單元的形狀相對(duì)較好，計(jì)算誤差相對(duì)較??；當(dāng)f＝3.0 時(shí)，應(yīng)力相對(duì)誤差指標(biāo)分別達(dá)到7.29×10-3和2.32×10-3.

表2 網(wǎng)格變化時(shí)應(yīng)力相對(duì)誤差指標(biāo)Tab.2 Relative error indexes of stresses for different fs

由表可見(jiàn)，應(yīng)力相對(duì)誤差指標(biāo)接近10-2.無(wú)論從實(shí)用角度還是理論角度，這樣的誤差都不是可以接受的，尤其是對(duì)于裂尖附近的含裂紋單元.這為Hansbo和Hansbo類方法應(yīng)用裂紋擴(kuò)展分析帶來(lái)不確定性.

3 分片試驗(yàn)失敗原因

無(wú)論是三角形還是四邊形數(shù)學(xué)單元，其離散位移函數(shù)能夠描述區(qū)域Ω＋和Ω-內(nèi)的簡(jiǎn)單拉伸變形.包括增強(qiáng)有限單元法在內(nèi)的Hansbo和Hansbo類方法完全具備描述非連續(xù)變形的能力.由于材料非不可壓縮，網(wǎng)格也沒(méi)有出現(xiàn)明顯畸形，可以排除整體剛度矩陣病態(tài)引起分片試驗(yàn)失敗.因此，數(shù)學(xué)單元和物理單元的幾何分離及剛度矩陣積分技術(shù)引起的剛度矩陣數(shù)值誤差可能是導(dǎo)致分片失敗的原因.

3.1 增強(qiáng)有限單元?jiǎng)偠染仃?/h3>

根據(jù)等參插值技術(shù)，4節(jié)點(diǎn)四邊形數(shù)學(xué)單元內(nèi)任意點(diǎn)的坐標(biāo)和位移可表示為

式中：a1，a2，…，d2為與數(shù)學(xué)單元節(jié)點(diǎn)（數(shù)學(xué)節(jié)點(diǎn)）坐標(biāo)有關(guān)的常數(shù)；e1，e2，…，h2為與數(shù)學(xué)單元節(jié)點(diǎn)位移線性相關(guān)的待求量；ξ和η 為數(shù)學(xué)單元的母單元坐標(biāo).

將ξη 項(xiàng)視為常數(shù)，聯(lián)立式（3）和（4）可得

利用式（3）和（4），數(shù)學(xué)單元的Jacobi矩陣J 及其逆矩陣J-1可表示為

將式（9）和（10）代入幾何方程，可得單元應(yīng)變?yōu)?/p>

式中：

其中，

利用虛功方程導(dǎo)出單元的剛度矩陣為

式中：ΩPE為物理單元所占據(jù)的空間區(qū)域，D 為材料的彈性系數(shù)矩陣.

將式（17）代入式（20），可得

其中

對(duì)于3節(jié)點(diǎn)三角形數(shù)學(xué)單元，不難證明，只需令ξ和η 為數(shù)學(xué)單元的面積坐標(biāo)，令d1，d2，h1，h2，，，和為零，上述推導(dǎo)仍然成立.此時(shí)，數(shù)學(xué)單元的Jacobi矩陣和Jacobi行列式與坐標(biāo)無(wú)關(guān)，且有

當(dāng)物理單元為4節(jié)點(diǎn)四邊形單元時(shí)，單元內(nèi)任意點(diǎn)的坐標(biāo)可插值表示為

式（31）～（34）一般采用Gauss積分進(jìn)行數(shù)值求積.

比較式（31）～（34）和式（35）～（38）可以看出，數(shù)學(xué)單元和物理單元分離使得剛度矩陣積分表達(dá)式更為復(fù)雜，表現(xiàn)在2個(gè)方面：

當(dāng)物理單元為3節(jié)點(diǎn)三角形時(shí)，式（22）～（25）改寫為如下形式，

3.2 剛度矩陣積分誤差

為分析Hansbo和Hansbo類方法不能以機(jī)器精度通過(guò)非連續(xù)分片試驗(yàn)的原因，本文根據(jù)數(shù)學(xué)單元和物理單元的類型，分4種情況分析剛度矩陣的積分誤差.不失一般性，假設(shè)一個(gè)單元只含一個(gè)物理子域，同時(shí)不考慮計(jì)算機(jī)的舍入誤差.

1）四邊形數(shù)學(xué)單元和四邊形物理單元

與Q4G2-M 型網(wǎng)格相比，Q4G2-S 型網(wǎng)格更為特殊，其數(shù)學(xué)單元和物理單元完全重合，此時(shí)ξ ＝、η＝，增強(qiáng)有限元單元退化為常規(guī)有限元單元.

2）四邊形數(shù)學(xué)單元和三角形物理單元

3）三角形數(shù)學(xué)單元和三角形物理單元

4）三角形數(shù)學(xué)單元和四邊形物理單元

4 數(shù)值算例

第3節(jié)從理論上證明了數(shù)學(xué)單元和物理單元分離及剛度矩陣積分技術(shù)引起單元?jiǎng)偠确e分精度降低，是導(dǎo)致非連續(xù)分片試驗(yàn)失敗的原因.本節(jié)通過(guò)三個(gè)算例從數(shù)值方面進(jìn)一步驗(yàn)證上述推導(dǎo)的正確性.

算例1 如圖7所示為直角梯形分析模型，包含1個(gè)常規(guī)有限單元，即物理單元與數(shù)學(xué)單元重合，幾何尺寸、邊界條件如圖所示.在模型的右豎直邊和斜邊上分別施加荷載q1＝1 MPa·m，q2＝1／MPa·m，設(shè)材料的E＝1 MPa，ν＝0.1.該問(wèn)題應(yīng)力的理論解為σx＝1，σy＝τxy＝0.

圖7 常應(yīng)力直角梯形Fig.7 Right trapezoid with constant stress

表3 應(yīng)力相對(duì)誤差指標(biāo)Tab.3 Relative error indexes of stresses

式中：σ 和ε 分別為應(yīng)力矢量和應(yīng)變矢量；上標(biāo)exact和h分別代表精確解和數(shù)值解.

圖8 常應(yīng)力矩形Fig.8 Rectangle with constant stresses

如圖9所示為相對(duì)誤差能量范數(shù)隨h的變化關(guān)系曲線.由圖可以看出，當(dāng)h＝1時(shí)，即數(shù)學(xué)單元為正方形時(shí)，相對(duì)誤差能量范數(shù)達(dá)到了10-16的機(jī)器精度；隨著h的增加或減少，相對(duì)誤差能量范數(shù)逐漸增加，當(dāng)h＝1.5 時(shí) 達(dá) 到3.5‰，當(dāng)h＝0.5 時(shí) 達(dá) 到2.2%.

圖9 相對(duì)誤差能量范數(shù)隨h的變化曲線Fig.9 Variation of relative error in energy norm with respect to h

保持h＝0.8不變，將原物理單元細(xì)劃成M×M個(gè)相同的矩形積分子域，每個(gè)子域仍采用2×2 點(diǎn)Gauss積分求單元?jiǎng)偠染仃?如圖10所示給出了相對(duì)誤差能量范數(shù)隨物理單元單邊細(xì)化次數(shù)M 的變化曲線.由圖可以看出，隨著M 的增加，相對(duì)誤差能量范數(shù)以幾乎恒定的速率迅速降低，當(dāng)M＝5時(shí)的相對(duì)誤差能量范數(shù)僅為M＝1 時(shí)的1.7‰.計(jì)算結(jié)果進(jìn)一步證明了剛度矩陣的數(shù)值計(jì)算誤差是導(dǎo)致分片試驗(yàn)失敗的原因.

圖10 相對(duì)誤差能量范數(shù)隨M 的變化曲線Fig.10 Variation of relative error in energy norm with respect to M

算例3 為驗(yàn)證細(xì)化積分子域?qū)Ψ沁B續(xù)分片試驗(yàn)的影響，將圖5（c）示模型中的含裂紋物理單元細(xì)化為M×M 個(gè)積分子域，當(dāng)M＝2時(shí)的積分子域劃分如圖11所示.

如表4所示給出了不同M 值時(shí)應(yīng)力相對(duì)誤差指標(biāo)計(jì)算結(jié)果.由表可見(jiàn)，隨著細(xì)化次數(shù)的增加，應(yīng)力計(jì)算精度明顯提高.

圖11 當(dāng)M＝2時(shí)物理單元積分子域劃分Fig.11 Subdomain partition of physical elements when M＝2

表4 不同細(xì)化次數(shù)時(shí)應(yīng)力相對(duì)誤差指標(biāo)Tab.4 Relative errors indexes of stresses with different Ms

5 結(jié) 論

本文以增強(qiáng)有限單元法為例，對(duì)Hansbo 和Hansbo類有限單元法進(jìn)行了非連續(xù)分片試驗(yàn)，理論和數(shù)值分析了不能以機(jī)器精度通過(guò)分片試驗(yàn)的原因，得到如下結(jié)論：

（1）采用3節(jié)點(diǎn)三角形數(shù)學(xué)單元的網(wǎng)格能以機(jī)器精度通過(guò)分片試驗(yàn)；除少數(shù)特殊情況外，采用4節(jié)點(diǎn)四邊形數(shù)學(xué)單元的網(wǎng)格不能以機(jī)器精度通過(guò)分片測(cè)試.

（2）數(shù)學(xué)單元與物理單元分離使得Hansbo和Hansbo類有限單元法具有表征單元內(nèi)部非連續(xù)變形的能力.然而，對(duì)于四邊形數(shù)學(xué)單元，數(shù)學(xué)單元與物理單元分離同時(shí)也使得剛度矩陣積分式中被積函數(shù)更為復(fù)雜，一般的三角形子域積分技術(shù)可能進(jìn)一步降低剛度矩陣的積分精度，導(dǎo)致單元不能以機(jī)器精度通過(guò)非連續(xù)分片試驗(yàn).

（3）為保證高精度通過(guò)非連續(xù)分片試驗(yàn)，合理的方法是在非連續(xù)變形區(qū)域采用三角形或平行四邊形數(shù)學(xué)單元，或?qū)⑴c數(shù)學(xué)單元不完全重合的物理單元精細(xì)劃分成積分子域以提高單元?jiǎng)偠染仃嚨挠?jì)算精度.需要指出，這些方法僅是從技術(shù)上提高了計(jì)算精度，若要從根本上解決問(wèn)題，需要一種無(wú)網(wǎng)格依賴性的積分方法，這有待進(jìn)一步研究.

（

）：

［1］AKSOYLU B，BOND S，HOLST M.An odyssey into local refinement and multilevel preconditioning III：Implementation and numerical experiments ［J］.SIAM Journal on Scientific Computing，2003，25（2）：478-498.

［2］ZHAO X，MAO S，SHI Z.Adaptive Finite Element Methods on quadrilateral meshes without hanging nodes［J］.SIAM Journal on Scientific Computing，2010，32（4）：2099-2120.

［3］SHI G H.Manifold method of material analysis［C］∥Transactions of the 9th army conference on applied mathematics and computing.Minneapolis，Minnesota：［s.n］，1991：57-76.

［4］位偉，姜清輝，周創(chuàng)兵.基于有限變形理論的數(shù)值流形方法研究［J］.力學(xué)學(xué)報(bào)，2014，46（1）：78-86.WEI Wei，JIANG Qing-hui，ZHOU Chuang-bing.Study on numerical manifold method based on finite deformation theory［J］.Chinese Journal of Theoretical and Applied Mechanics，2014，46（1）：78-86.

［5］徐棟棟，鄭宏，夏開(kāi)文，等.高階擴(kuò)展數(shù)值流形法在裂紋擴(kuò)展中的應(yīng)用.巖石力學(xué)與工程學(xué)報(bào)［J］，2014，33（7）：1375-1387.XU Dong-dong，ZHENG Hong，XIA Kai-wen，et al.Application of higher-order enriched numerical manifold method to crack propagation［J］.Chinese Journal of Rock Mechanics and Engineering，2014，33（7）：1375-1387.

［6］蘇海東，祁勇峰，龔亞琦，等.任意形狀覆蓋的數(shù)值流形方法初步研究.長(zhǎng)江科學(xué)院院報(bào)［J］，2013，30（12）：91-96.SU Hai-dong，QI Yong-feng，GONG Ya-qi，et al.Preliminary research of numerical manifold method based on covers of arbitrary shape［J］.Journal of Yangtze River Scientific Research Institute，2013，30（12）：91-96.

［7］BABU? KA I，MELENK J M.The partition of unity method［J］.International Journal for Numerical methods in Engineering，1997，40（4）：727-758.

［8］MO?S N，DOLBOW J，BELYTSCHKO T.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46：131-150.

［9］BELYTSCHKO T，BLACK T.Elastic crack growth in finite elements with minimal remeshing［J］.International Journal for Numerical Methods in Engineering，1999，45：601-620.

［10］楊志鋒，周昌玉，代巧.基于擴(kuò)展有限元法的彈塑性裂紋擴(kuò)展研究.南京工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版［J］，2014，36（4）：50-57.YANG Zhi-feng，ZHOU Chang-yu，DAI Qiao.Elasticplastic crack propagation based on extended finite element method［J］.Journal of Nanjing University of Technology：Natural Science Edition，2014，36（4）：50-57.

［11］師訪，高峰，楊玉貴.正交各向異性巖體裂紋擴(kuò)展的擴(kuò)展有限元方法研究［J］.巖土力學(xué)，2014，35（4）：1203-1210.SHI Fang，GAO Feng，YANG Yu-gui.Application of extended finite element method to study crack propagation problems of orthotr opic rock mass［J］.Rock and Soil Mechanics，2014，35（4）：1203-1210.

［12］HANSBO A，HANSBO P.A finite element method for the simulation of strong and weak discontinuities in solid mechanics［J］.Computer Methods in Applied Mechanics and Engineering，2004，193：3523-3540.

［13］HANSBO A，HANSBO P.An unfitted finite element method，based on Nitsche's method for elliptic interface problems［J］.Computer Methods in Applied Mechanics and Engineering，2002，191：5537-5552.

［14］AREIAS P，BELYTSCHKO T.A comment on the article“A finite element method for simulation of strong and weak discontinuities in solid mechanics”by A.Hansbo and P.Hansbo ［Comput.Methods Appl.Mech.Engrg.193（2004）3523–3540］［J］.Computer Methods in Applied Mechanics and Engineering，2006，195（9）：1275-1276.

［15］SONG J H，AREIAS P，BELYTSCHKO T.A method for dynamic crack and shear band propagation with phantom nodes［J］.International Journal for Numerical Methods in Engineering，2006，67（6）：868-893.

［16］LING D S，YANG Q D，COX B N.An augmented finite element method for modeling arbitrary discontinuities in composite materials［J］.International Journal of Fracture，2009，156（1）：53-73.

［17］凌道盛，卜令方，涂福彬.粘聚裂紋擴(kuò)展的強(qiáng)化有限元h型網(wǎng)格自適應(yīng)模擬［J］.計(jì)算力學(xué)學(xué)報(bào)，2014，31（2）：241-247.LING Dao-sheng，BU Ling-fang，TU Fu-bing.Modelling of cohesive crack propagation using enhanced finite element method via h-adaptive technique［J］.Chinese Journal of Computational Mechanics，2014，31（2）：241-247.

［18］LING D S，BU L F，TU F B，et al.A finite element method with mesh-separation based approximation technique and its application in modeling crack propagation with adaptive mesh refinement［J］.International Journal for Numerical Methods in Engineering，2014，99（7）：487-521.

［19］DOLBOW J E，DEVAN A.Enrichment of enhanced assumed strain approximations for representing strong discontinuities：addressing volumetric incompressibility and the discontinuous patch test［J］.International Journal for Numerical Methods in Engineering，2004，59：47-67.