常 鑫，程遠(yuǎn)方，夏強(qiáng)平，韓修廷

（1.中國石油大學(xué)石油工程學(xué)院，山東青島266580；2.清華大學(xué)航空航天與工程力學(xué)系，北京100084）

一種模擬巖體裂紋擴(kuò)展的三角單元網(wǎng)格開裂技術(shù)

常 鑫1，程遠(yuǎn)方1，夏強(qiáng)平2，韓修廷1

（1.中國石油大學(xué)石油工程學(xué)院，山東青島266580；2.清華大學(xué)航空航天與工程力學(xué)系，北京100084）

基于三角網(wǎng)格的幾何特征，提出一種利用有限元方法模擬巖體裂紋擴(kuò)展的三角單元網(wǎng)格開裂技術(shù)。該方法選取三角網(wǎng)格進(jìn)行單元離散，采用遠(yuǎn)場圍線積分計(jì)算裂尖應(yīng)力強(qiáng)度因子，由最大周向應(yīng)力準(zhǔn)則確定裂紋擴(kuò)展方向，最后通過開裂單元的網(wǎng)格分裂或節(jié)點(diǎn)移動(dòng)，實(shí)現(xiàn)裂紋擴(kuò)展的數(shù)值模擬。以有限寬中心裂紋板、曲線翼型裂紋擴(kuò)展和含孔洞多裂隙巖體的裂紋擴(kuò)展為例進(jìn)行模擬驗(yàn)證。結(jié)果表明：在該方法中，裂紋可以直接劈開一個(gè)單元，或沿單元邊界擴(kuò)展，因此裂紋能夠不受初始網(wǎng)格的限制沿任意路徑擴(kuò)展；與現(xiàn)有的網(wǎng)格重構(gòu)算法相比，該方法只須對裂尖局部單元進(jìn)行網(wǎng)格開裂或節(jié)點(diǎn)移動(dòng)，更加簡便、高效，該方法還具有較好的適用性，能夠準(zhǔn)確模擬拉伸、壓剪等復(fù)雜應(yīng)力狀態(tài)下的裂紋萌生和擴(kuò)展。

有限元法；三角網(wǎng)格開裂；裂紋擴(kuò)展；數(shù)值模擬

巖體既非完全的連續(xù)體，也非完全的離散體，而是富含節(jié)理、裂隙等軟弱結(jié)構(gòu)面的地質(zhì)缺陷體，極易發(fā)生脆性斷裂。大量的巖土工程實(shí)踐表明，巖體的失穩(wěn)破壞與其內(nèi)部節(jié)理、裂隙的擴(kuò)展和貫通密切相關(guān)［1-8］。巖體裂紋擴(kuò)展的實(shí)質(zhì)是連續(xù)面到非連續(xù)面的轉(zhuǎn)換，如何實(shí)現(xiàn)對不連續(xù)面的準(zhǔn)確刻畫是問題的關(guān)鍵。傳統(tǒng)有限元方法（FEM）［9-12］模擬裂紋擴(kuò)展時(shí)通常需要在裂尖構(gòu)建奇異單元，且隨著裂紋擴(kuò)展不斷地進(jìn)行網(wǎng)格重構(gòu)，該方法實(shí)現(xiàn)困難且低效。邊界元法（BEM）［13-14］雖然避免了網(wǎng)格重新劃分的問題，但嚴(yán)重依賴于問題的基本解，對于非線性、非均質(zhì)等問題其優(yōu)勢大大降低。針對以上不足，研究者提出擴(kuò)展有限元方法（XFEM）［15-17］，該方法在模擬裂紋擴(kuò)展時(shí)無需重新劃分網(wǎng)格，而是在原有網(wǎng)格方案基礎(chǔ)上采用特殊的插值技術(shù)處理裂紋問題，雖然較傳統(tǒng)FEM有了較大改進(jìn)，但在處理復(fù)雜裂紋擴(kuò)展時(shí)富集函數(shù)通常難以確定，且引入了附加自由度，計(jì)算量較大。筆者從三角單元幾何特征出發(fā)，提出一種模擬巖體裂紋擴(kuò)展的三角單元網(wǎng)格開裂技術(shù)，該技術(shù)只對裂尖局部單元進(jìn)行網(wǎng)格調(diào)整，不涉及不連續(xù)位移場的處理問題，也無須引入額外自由度，并將數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果進(jìn)行對比，驗(yàn)證三角單元網(wǎng)格開裂方法的有效性。

1　基本原理

網(wǎng)格開裂技術(shù)的基本流程可以概括為：選取三角形非結(jié)構(gòu)化網(wǎng)格進(jìn)行單元離散；采用遠(yuǎn)場圍線積分獲取裂尖應(yīng)力強(qiáng)度因子并確定裂紋擴(kuò)展方向；通過開裂單元的網(wǎng)格分裂或節(jié)點(diǎn)移動(dòng)，實(shí)現(xiàn)裂紋擴(kuò)展的數(shù)值模擬。

1.1裂隙巖體的虛功方程

考慮裂隙巖體Ω，線彈性小變形，其應(yīng)力和位移邊界分別為Γσ和Γu，裂紋體Ω內(nèi)含多條裂紋Γcr，如圖1所示。

圖1　裂隙巖體受力示意圖Fig.1 Stress diagram of rock mass with cracks

根據(jù)虛位移原理，平衡方程的等效積分形式為

分布積分后可得如下弱形式

采用最小二乘有限元法求應(yīng)力，定義泛函：

式中，σ0（u）為由變形算得的應(yīng)力；σ為待求的應(yīng)力。由變分原理知δF（σ）=0，可得

1.2應(yīng)力強(qiáng)度因子計(jì)算及斷裂準(zhǔn)則確定

目前，求解復(fù)雜幾何形狀和復(fù)雜載荷下的混合型應(yīng)力強(qiáng)度因子主要有權(quán)函數(shù)法、虛擬裂紋擴(kuò)展法、邊界配置法和J積分法等［18］，但上述方法一般對裂尖網(wǎng)格質(zhì)量要求較高。本文中采用基于Betti功能互等定理的遠(yuǎn)場圍線積分法計(jì)算裂尖應(yīng)力強(qiáng)度因子，該方法只需要知道遠(yuǎn)離裂尖處某一積分路徑上的應(yīng)力解與位移解，便可計(jì)算出裂尖應(yīng)力強(qiáng)度因子，計(jì)算結(jié)果不受裂紋幾何形狀的影響，而且在遠(yuǎn)場積分，裂紋尖端不需要采用特殊類型的單元，用較粗的有限元網(wǎng)格也可獲得足夠的精度。圍線積分的基本原理如下。

不計(jì)體力時(shí)，均勻各向同性彈性體的Betti功能互等定理為

式中，Γ為彈性體內(nèi)某閉合曲線；ui和ti為某一平衡狀態(tài)下邊界Γ上的位移和作用力；^ui和^ti為另一輔助平衡狀態(tài)下圍線Γ上的位移和作用力。

對于含裂紋體，積分圍線示意圖如圖2所示。

圖2　積分圍線示意圖Fig.2 Principle diagram of contour integral

取閉合曲線Γ=Γδ′+AB+（-Γδ）+CD作為圍線Γ。當(dāng)裂紋面為自由面時(shí)，Betti功能互等定理等價(jià)于如下形式：

若圍線Γδ取半徑為δ的圓的邊界，則令

式中，G為剪切模量；k=（3-ν）/（1+ν）（平面應(yīng)力）或k=3-4ν，ν為材料的泊松比。

由于Γδ′為遠(yuǎn)場繞裂端的圍線，ui和ti取有限元數(shù)值解，^ui和^ti分別取輔助位移場（^u^v）與輔助應(yīng)力場（^σx^σy^τxy），在圖2所示的局部坐標(biāo)系下的表達(dá)式為

采用復(fù)化Simpson計(jì)算遠(yuǎn)場圍線積分，其解為

將式（7）和式（9）代入式（6）可得

式中，KⅠ和KⅡ?yàn)閼?yīng)力強(qiáng)度因子。

比較CⅠ和CⅡ前面的系數(shù)即可得混合型裂紋的應(yīng)力強(qiáng)度因子。

根據(jù)最大周向應(yīng)力準(zhǔn)則，裂紋開裂方向角θC為

當(dāng)最大周向應(yīng)力大于等于臨界值時(shí)，裂紋繼續(xù)擴(kuò)展，此時(shí)有

1.3三角網(wǎng)格開裂算法

網(wǎng)格開裂算法遵循兩大基本原則：首先，保證新舊網(wǎng)格之間具有較高的相容性；其次，保證對初始網(wǎng)格較低的依賴性。該算法主要包含破裂單元確定及單元分裂兩大模塊。

（1）確定破壞單元。根據(jù)最大周向應(yīng)力準(zhǔn)則計(jì)算單元破壞方向，通過矢量叉乘運(yùn)算，搜索破裂方向通過單元，將該單元記為破裂單元。

（2）分裂破壞單元。沿著破裂方向與破裂單元對邊交點(diǎn)連線將破裂單元進(jìn)行單元剖分，如圖3所示。

圖3　三角網(wǎng)格開裂技術(shù)原理示意圖Fig.3 Principle diagram of triangular mesh split method

詳細(xì)的計(jì)算流程如下：

（1）通過遠(yuǎn)場圍線積分，計(jì)算裂尖應(yīng)力強(qiáng)度因子KⅠ、KⅡ以及裂紋開裂方向角θC，根據(jù)復(fù)合裂紋斷裂準(zhǔn)則確定開裂點(diǎn)。

（2）搜索包含開裂點(diǎn)的所有單元及單元信息。

（3）記錄開裂點(diǎn)周圍的邊界信息（裂紋）即上游開裂點(diǎn)信息。

（4）通過裂紋開裂方向矢量與NECP中三角形單元以開裂點(diǎn)為起點(diǎn)的兩條邊的方向矢量做矢量叉乘運(yùn)算確定開裂單元。

（5）增加與開裂點(diǎn)重合的節(jié)點(diǎn)，并求取開裂方向與開裂單元對邊的交點(diǎn)。

（6）沿裂紋開裂方向，虛擬劈開單元，依據(jù)Delaunay準(zhǔn)則判斷網(wǎng)格的畸變性。

（7）若存在畸變網(wǎng)格，調(diào)整網(wǎng)格將網(wǎng)格棱邊移至裂紋開裂方向，如圖4所示。假設(shè)開裂方向?yàn)?，棱邊為，則將B點(diǎn)坐標(biāo)替換為E點(diǎn)的值，同時(shí)將和作為裂紋邊界，其中A′為第（5）步中新增節(jié)點(diǎn)。

（8）若不存在畸變網(wǎng)格，則沿開裂方向?qū)㈤_裂單元分裂為兩個(gè)網(wǎng)格，增加節(jié)點(diǎn)E（圖4），將和記為裂紋邊界，同時(shí)將臨開裂單元ΔBCD，在E點(diǎn)分裂為兩個(gè)網(wǎng)格。

圖4　網(wǎng)格優(yōu)化Fig.4 Mesh optimization

（9）最后修正單元節(jié)點(diǎn)編號，順時(shí)針側(cè)單元用節(jié)點(diǎn)A表示，逆時(shí)針側(cè)單元用節(jié)點(diǎn)A′表示，同時(shí)為新單元映射初始信息。

（10）重復(fù)（1）～（9），進(jìn)行下一開裂點(diǎn)網(wǎng)格開裂操作。

網(wǎng)格調(diào)整方案如圖5所示。

圖5　裂紋擴(kuò)展過程中的網(wǎng)格調(diào)整方案Fig.5 Grid adjustment programs of crack propagation process

2　裂紋擴(kuò)展模擬

2.1有限寬中心裂紋板

對于特殊的含裂紋構(gòu)件，已有學(xué)者求得了閉合形式解（以公式、曲線或圖表的形式表示），并匯編成應(yīng)力強(qiáng)度因子手冊。以有限寬中心裂紋板為例對本程序斷裂參數(shù)的計(jì)算精度進(jìn)行驗(yàn)證?？紤]如圖6所示的有限寬中心裂紋板，板的半高H、半寬w和厚度B分別為200、100和1.0 mm，中心裂紋總長2a為40 mm。在板的遠(yuǎn)處承受均勻的拉應(yīng)力σ=30 MPa。假設(shè)板各向同性、均勻、線彈性，彈性模量E= 200 GMPa，泊松比ν=0.3。應(yīng)力強(qiáng)度［19］表達(dá)式為

其中F為有限寬度修正系數(shù)，典型的計(jì)算公式有以下3種。

Isida公式的近似：

修正的Koiter公式：

修正的Feddersen公式：

考慮模型對稱性，取1/2模型進(jìn)行分析，如圖7所示。采用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行單元離散（共11893個(gè)節(jié)點(diǎn)和23240個(gè)單元）。

圖6　有限寬中心裂紋板Fig.6 Cracked plate of finite width center

圖7　網(wǎng)格劃分方案（1/2模型）Fig.7 Program of grid generation（1/2 model）

圖8　裂紋體變形（×200倍）及裂尖前緣方向的應(yīng)力云圖Fig.8 Crack deformation（×200）and crack tip σystress nephogram

表1　無量綱化的應(yīng)力強(qiáng)度因子Table 1 Dimensionless stress intensity factor

2.2壓剪狀態(tài)下曲線翼型裂紋擴(kuò)展特性模擬

翼型裂紋的起裂、擴(kuò)展及其失穩(wěn)特性一直是研究脆性裂隙巖體漸近失穩(wěn)破壞機(jī)制的基礎(chǔ)，國內(nèi)外眾多學(xué)者對壓剪狀態(tài)下的翼型裂紋擴(kuò)展進(jìn)行了試驗(yàn)研究。為了驗(yàn)證本文中方法的準(zhǔn)確性，將數(shù)值模擬結(jié)果與楊慶等［1］的試驗(yàn)結(jié)果進(jìn)行對比。試驗(yàn)采用熟石膏作為試件材料，其彈性模量E=0.682 GPa，泊松比ν=0.32，抗拉強(qiáng)度為0.384 MPa，斷裂韌度KIC=0.032 MPa·m1/2。試件外形尺寸為250 mm× 180 mm×20 mm，如圖9所示。

圖9　裂紋板幾何模型Fig.9 Geometric model of cracked plate

在試件中預(yù)埋薄金屬片，金屬片的寬度為0.2 mm，當(dāng)石膏凝結(jié)后拔出金屬片，形成含預(yù)制裂紋試件，裂紋厚度為0.2 mm，裂紋的深度為20 mm。模擬不同裂紋角條件下裂紋的擴(kuò)展路徑，基本參數(shù)見表2。

表2　試樣編號及裂紋特征Table 2 Sample number and crack characteristics

數(shù)值模型尺寸與試樣尺寸相同，邊界采用滑動(dòng)邊界，利用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行單元離散。數(shù)值模擬開始后，計(jì)算當(dāng)前載荷作用下的裂紋尖端應(yīng)力強(qiáng)度因子KⅠ和KⅡ，并進(jìn)行裂紋擴(kuò)展判斷，若滿足式（13），則對裂紋尖端相應(yīng)單元進(jìn)行剖分，整個(gè)數(shù)值計(jì)算過程通過Mesh_operate.for程序自動(dòng)實(shí)現(xiàn)。圖10為試樣I-2計(jì)算40步后的網(wǎng)格拓?fù)浣Y(jié)構(gòu)。圖10中紫色曲線為裂紋擴(kuò)展路徑，與原始網(wǎng)格結(jié)構(gòu)相比僅裂紋擴(kuò)展路徑附近單元發(fā)生了網(wǎng)格剖分，新生成的網(wǎng)格依然保持較高的網(wǎng)格質(zhì)量。以試樣I-2為例給出了翼型裂紋擴(kuò)展的有限元實(shí)體模型，如圖11所示。

圖10　網(wǎng)格開裂后新裂紋拓?fù)浣Y(jié)構(gòu)（I-2）Fig.10 New crack topology structure after grid cracking（I-2）

圖12為不同初始裂紋角條件下翼型裂紋形狀及擴(kuò)展角，給出了翼型裂紋擴(kuò)展路徑及擴(kuò)展方向隨裂紋擴(kuò)展長度的變化關(guān)系。由圖12可以看出：在翼型裂紋擴(kuò)展中存在一條過初始裂紋中點(diǎn)平行于最大主應(yīng)力方向的直線，隨著裂紋長度增加，翼型裂紋不斷逼近這條直線；同時(shí)發(fā)現(xiàn)在裂紋擴(kuò)展初期，裂紋擴(kuò)展方向變化較大，隨著裂紋長度增加，裂紋擴(kuò)展方向逐漸逼近90°。

圖11　壓剪作用下翼型裂紋擴(kuò)展有限元模擬結(jié)果（I-2）Fig.11 Crack extension finite element simulation results under compression-shear action（I-2）

圖12　不同初始裂紋角條件下翼型裂紋形狀和擴(kuò)展角Fig.12 Crack path and propagation angle under different initial crack angles

2.3含孔洞多裂紋巖體開裂過程模擬

巖石材料具有明顯的非均質(zhì)性，內(nèi)部含有多種尺度、多種形狀的孔洞、裂隙以及結(jié)構(gòu)面，在工程計(jì)算領(lǐng)域?qū)锥炊嗔鸭y巖體開裂過程的模擬一直是一個(gè)難點(diǎn)問題。已有的研究大都集中在單一孔洞模型，對含孔洞多裂紋巖體的開裂過程和孔洞的存在對裂紋擴(kuò)展路徑的影響研究較少［20］，以雙孔洞裂隙巖體為例對上述問題進(jìn)行數(shù)值模擬，著重分析孔洞對裂紋擴(kuò)展路徑的影響，以證明網(wǎng)格開裂算法在模擬多裂紋擴(kuò)展方面的適用性。模型結(jié)構(gòu)如圖13所示。

模型尺寸為D×L（127 mm×127 mm），厚度忽略不計(jì)，模型簡化為平面應(yīng)力問題。圓孔洞的半徑為2 mm，初始裂紋長度a=2.54 mm，傾角θ=45°。其中：彈性模量E=0.689 GPa，泊松比ν=0.3，受雙向拉伸載荷作用。仍然利用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行單元離散，圖14給出了計(jì)算30步后的右側(cè)孔洞附近的網(wǎng)格拓?fù)浣Y(jié)構(gòu)，圖14中紫色曲線為裂紋擴(kuò)展路徑，其中內(nèi)側(cè)裂紋（I-1、I-2號裂紋）在孔壁附近明顯沿曲線路徑擴(kuò)展。

圖13　含孔洞多裂紋體幾何模型Fig.13 Geometric model of multiple cracks body with holes

圖15為裂紋擴(kuò)展過程中的第一主應(yīng)力云圖。從圖15（a）可以看出，當(dāng)內(nèi)側(cè)裂紋相遇后在孔洞內(nèi)側(cè)產(chǎn)生顯著的應(yīng)力集中，裂紋擴(kuò)展發(fā)生明顯偏折，呈現(xiàn)出“相互吸引”的特征，隨著裂紋進(jìn)一步擴(kuò)展，內(nèi)側(cè)應(yīng)力集中區(qū)域逐漸消失，相反外側(cè)裂紋尖端處應(yīng)力集中明顯，如圖15（b）所示。內(nèi)側(cè)裂紋同時(shí)受孔洞和原始裂隙的雙重影響，擴(kuò)展路徑迂回曲折，而外側(cè)裂紋呈現(xiàn)出近乎直線的擴(kuò)展形態(tài)。圖16為本文數(shù)值模擬結(jié)果與Moes等［21］的研究成果的對比，可以看出兩者具有極高的吻合度。

圖14　網(wǎng)格開裂后新網(wǎng)格拓?fù)浣Y(jié)構(gòu)Fig.14 New mesh topology structure after grid cracking

圖15　裂紋擴(kuò)展過程第一主應(yīng)力云圖（×2）Fig.15 The first principal stress nephogram during crack propagation process（×2）

圖16　含孔洞多裂紋巖體裂紋擴(kuò)展路徑Fig.16 Crack propagation path of rock mass with holes

3　結(jié) 論

（1）采用三角單元網(wǎng)格開裂技術(shù)模擬裂紋擴(kuò)展時(shí)，裂紋可以直接劈開一個(gè)單元或沿單元邊界擴(kuò)展，因此裂紋能夠不受初始網(wǎng)格的限制沿任意路徑擴(kuò)展；與現(xiàn)有的網(wǎng)格重構(gòu)算法相比，該方法只須對裂尖局部單元進(jìn)行網(wǎng)格開裂或節(jié)點(diǎn)移動(dòng)，更加簡便、高效。

（2）三角單元網(wǎng)格開裂技術(shù)不僅具有較高的斷裂參數(shù)計(jì)算精度，而且可以準(zhǔn)確地模擬拉伸、壓剪等復(fù)雜應(yīng)力狀態(tài)下的裂紋萌生和擴(kuò)展，將為巖體基本力學(xué)性質(zhì)以及巖體工程方法的基礎(chǔ)研究提供有力支撐。

［1］楊慶，李強(qiáng)，趙維，等.翼型裂紋擴(kuò)展特性的試驗(yàn)和數(shù)值分析［J］.水利學(xué)報(bào)，2009，40（8）：934-940. YANG Qing，LI Qiang，ZHAO Wei，et al.Experiment study and numerical simulation on propagation characteristrics of wing type cracks under compression［J］.Journal of Hgdraulic Engineering，2009，40（8）：934-940.

［2］陳衛(wèi)忠，李術(shù)才，朱維申，等.巖石裂紋擴(kuò)展的實(shí)驗(yàn)與數(shù)值分析研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2003，22（1）：18-23. CHEN Weizhong，LI Shucai，ZHU Weishen，et al.Experimental and numerical research on crack propagation in rock under compression［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（1）：18-23.

［3］周火明，楊宇，張宜虎，等.多裂紋巖石單軸壓縮漸進(jìn)破壞過程精細(xì)測試［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29（3）：465-470. ZHOU Huoming，YANG Yu，ZHANG Yihu，et al.Fine test on progressive fracturing process of multi-crack rock samples under uniaxial compression［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（3）：465-470.

［4］李世愚，和泰名，尹祥礎(chǔ)，等.巖石斷裂力學(xué)導(dǎo)論［M］.合肥：中國科學(xué)技術(shù)大學(xué)出版社，2010：1-10.

［5］鄭哲敏.連續(xù)介質(zhì)力學(xué)與斷裂［J］.力學(xué)進(jìn)展，1982，12（2）：133-140. ZHENG Zhemin.Continuum mechanics and fracture［J］. Advances in Mechanics，1982，12（2）：133-140.

［6］唐春安，黃明利，張國民，等.巖石介質(zhì)中多裂紋擴(kuò)展相互作用及其貫通機(jī)制的數(shù)值模擬［J］.地震，2001，21（2）：53-58. TANG Chun'an，HUANG Mingli，ZHANG Guomin，et al.Numerical simulation on propagation interaction and coalescence of multi-cracks in rocks［J］.Earthquake，2001，21（2）：53-58.

［7］李強(qiáng)，張力霆，齊清蘭，等.壓縮載荷下閉合裂紋的曲線分支裂紋模型研究［J］.工程力學(xué)，2013，30（4）：294-300. LI Qiang，ZHANG Liting，QI Qinglan，et al.Study of the curve branch crack model for the closed flaw under compressive loading［J］.Engineering Mechanics，2013，30（4）：294-300.

［8］陳勉.頁巖氣儲層水力裂縫轉(zhuǎn)向擴(kuò)展機(jī)制［J］.中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，37（5）：88-94. CHEN Mian.Re-orientation and propagation of hydraulic fractures in shale gas reservoir［J］.Journal of China University of Petroleum（Edition of Natural Science），2013，37（5）：88-94.

［9］楊慶生，楊衛(wèi).斷裂過程的有限元模擬［J］.計(jì)算力學(xué)學(xué)報(bào)，1997，14（4）：33-38. YANG Qingsheng，YANG Wei.Finite element simulation of fracture process［J］.Chinese Journal of Computational Mechanics，1997，14（4）：33-38.

［10］孫翔，劉傳奇，薛世峰.有限元與離散元混合法在裂紋擴(kuò)展中的應(yīng)用［J］.中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，37（3）：126-130，136. SUN Xiang，LIU Chuanqi，XUE Shifeng.Application of combined finite-discrete element method for crack propagation［J］.Journal of China University of Petroleum（E-dition of Natural Science），2013，37（3）：126-130，136.

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

［12］DAVID A，MARCELO E，MARIA C R.Automatic LEFM crack propagation method based on local lepp-delaunay mesh refinement［J］.Advances in Engineering Software，2010，41：111-119.

［13］雷鈞，黃祎豐，楊慶生.裂紋動(dòng)態(tài)擴(kuò)展的邊界元模擬［J］.北京工業(yè)大學(xué)學(xué)報(bào)，2013，39（6）：806-810. LEI Jun，HUANG Weifeng，YANG Qingsheng.Dynamic crack propagation by boundary element method［J］. Journal of Beijing University of Technology，2013，39（6）：806-810.

［14］秦飛，岑章志，杜慶華.多裂紋擴(kuò)展分析的邊界元方法［J］.固體力學(xué)學(xué)報(bào)，2002，23（4）：431-438. QIN Fei，CEN Zhangzhi，DU Qinghua.Boundary element analysis of cracked bodies［J］.Acta Mechanica Solida Sinica，2002，23（4）：431-438.

［15］郭歷倫，陳忠富，羅景潤，等.擴(kuò)展有限元方法及應(yīng)用綜述［J］.力學(xué)季刊，2011，32（4）：612-615. GUO Lilun，CHEN Zhongfu，LUO Jingrun，et al.A review of the extended finite element method and its applications［J］.Chinese Quarterly of Mechanics，2011，32（4）：612-615.

［16］方修君，金峰，王進(jìn)廷.基于擴(kuò)展有限元法的粘聚裂紋模型［J］.清華大學(xué)學(xué)報(bào)：自然科版，2007，47（3）：344-347. FANG Xiujun，JIN Feng，WANG Jinting.Cohesive crack model based on extended finite element method［J］.J Tsing Univ（Sci&Tech），2007，47（3）：344-347.

［17］劉劍，劉曉波，馮發(fā)強(qiáng).基于擴(kuò)展有限元法對含孔洞平板多裂紋擴(kuò)展模擬研究［J］.失效分析與預(yù)防，2013，8（6）：326-330. LIU Jian，LIU Xiaobo，F(xiàn)ENG Faqiang.Simulation study of multiple crack extension of holey plate based on extended finite element method［J］.Failure Analysis and Prevention，2013，8（6）：326-330.

［18］楊曉翔，范家齊，匡震邦.求解混合型裂紋應(yīng)力強(qiáng)度因子的圍線積分法［J］.計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，1996，13（1）：84-89. YANG Xiaoxiang，F(xiàn)AN Jiaqi，KUANG Zhenbang.A contour integral method for stress intensity factors of mixed-mode crack［J］.Computational Structural Mechanics and Applications，1996，13（1）：84-89.

［19］中國航空研究院.應(yīng)力強(qiáng)度因子手冊［M］.北京：科學(xué)出版社，1993：128-129.

［20］楊圣奇，呂朝輝，渠濤.含單個(gè)孔洞大理巖裂紋擴(kuò)展細(xì)觀試驗(yàn)和模擬［J］.中國礦業(yè)大學(xué)學(xué)報(bào)，2009，38（6）：774-781. YANG Shengqi，Lü Chaohui，QU Tao.Investigations of crack expansion in marble having a single pre-existing hole：experiment and simulations［J］.Journal of China University of Mining&Technology，2009，38（6）：774-781.

［21］NICOLAS M，JOHN D，TED B.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46：131-150.

（編輯 李志芬）

A triangular mesh split method for simulating crack propagation in rock matrix

CHANG Xin1，CHENG Yuanfang1，XIA Qiangping2，HAN Xiuting1
（1.School of Petroleum Engineering in China University of Petroleum，Qingdao 266580，China；2.Department of Engineering Mechanics in Tsinghua University，Beijing 100084，China）

Based on the geometrical characteristics of triangular meshes，a triangular mesh split method was proposed to simulate the propagation of cracks in rock matrix.Firstly，the rock matrix with primary cracks was divided by triangular elements. Then the stress intensity factors were calculated using a far-field contour integral，and the orientation of the crack propagation can be determined via the maximum circumferential stress criterion.Finally，the numerical simulation of crack propagation was realized by the mesh splitting or node movement.The applicability of the triangular mesh split algorithm was verified in case studies，including the crack propagation in a centrally cracked panel with a finite width，along curved wings and in a multi-fractured rock matrix with holes.The results show that the crack can split a unit element or can be extended along the mesh boundary，therefore the crack can be extended along arbitrary path without the limitation of the original mesh.Compared with the existing mesh reconstruction algorithm，this method is more convenient and efficient for use，and it can be also used effectively to simulate the initiation and propagation of cracks under complex stresses including tension，compression and shearing.

finite element method；triangular mesh split；crack propagation；numerical simulation

TU 457

A

1673-5005（2015）03-0105-08

10.3969/j.issn.1673-5005.2015.03.014

2014-07-26

“十二五”國家科技重大專項(xiàng)（2011ZX05037-004）；教育部長江學(xué)者和創(chuàng)新團(tuán)隊(duì)發(fā)展計(jì)劃（IRT1086）

常鑫（1987-），男，博士研究生，主要從事非常規(guī)儲層水力壓裂裂紋擴(kuò)展形態(tài)研究。E-mail：changxin7521@163.com。

引用格式：常鑫，程遠(yuǎn)方，夏強(qiáng)平，等.一種模擬巖體裂紋擴(kuò)展的三角單元網(wǎng)格開裂技術(shù)［J］.中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2015，39（3）：105-112.

CHANG Xin，CHENG Yuanfang，XIA Qiangping，et al.A triangular mesh split method for simulating crack propagation in rock matrix［J］.Journal of China University of Petroleum（Edition of Natural Science），2015，39（3）：105-112.