包含摩爾-庫侖準則的單元劈裂法模擬圍壓下節理擴展

楊 帆， 張振南

(1.上海大學土木工程系，上海200072;2.上海交通大學船舶海洋與建筑工程學院，上海200240)

巖體內部含有大量節理，一般在計算時將節理 簡化成裂紋.對于節理擴展匯合模式，國內外學者［1-6］已進行了大量的理論與試驗研究.對于主要節理，一般采用節理單元法來再現節理面之間的相互作用，其中Goodman類型的節理單元［5］能夠很好地模擬節理的力學特性.但是在應用該方法時，有限元網格劃分需要事先考慮節理分布和幾何形狀，以便設置節理單元.另外，在節理擴展時還要重新修正網格，以便設置節理單元，使模擬過程得以繼續，這些都給節理擴展過程的數值模擬帶來很大不便.為了避免前處理網格劃分問題和節理擴展后網格修正問題，張振南等［7-8］提出了單元劈裂法(element partition method，EPM)，該方法利用常應變三角單元的幾何特點構造了一種特殊的三節點接觸單元，用以再現節理面之間的接觸和摩擦效應，還可以在原有網格劃分方案的基礎上直接對裂紋擴展進行模擬，使得節理擴展數值模擬更為簡單、高效.目前國際上流行的擴展有限元法(extended finite element method，XFEM)［9-10］也可以在不用重新劃分網格的基礎上對裂紋擴展進行數值模擬.與EPM相比，XFEM更為精確，但計算過程復雜;而EPM的精確度取決于網格尺寸，計算過程簡單，因而為實際工程的預先粗略計算提供了很好的計算手段.

摩爾-庫侖(Mohr-Coulomb)準則在巖土工程界是一個廣泛使用的破壞準則.當材料微元滿足該準則時，微元發生破壞，相當于在微元層面上形成一個裂紋微段.本研究試圖以摩爾-庫侖準則作為一種單元劈裂準則，應用單元劈裂法對巖體節理擴展行為進行模擬，以探索這種方法在節理擴展模擬過程中的有效性.

1 單元劈裂法

單元劈裂法［7］利用三角單元的幾何特點構造特殊的三節點接觸單元(見圖1)，即當有裂紋貫穿時，總有一個節點(節點3)在裂紋的一側，另外2個節點(節點1和2)在裂紋的另一側.在裂紋一側的一個節點可與另外2個節點分別構成2個接觸點對.根據這2個接觸點對可以推導出三節點接觸單元的剛度矩陣［7］為

圖1 受裂紋切割的單元網格與三節點接觸單元Fig.1 Mesh intersected by crack and the tri-node contact element

式中，Kn，Ks分別為法向剛度系數和切向剛度系數，L為節理與三角單元相交的長度，h為節理厚度，Hij為表示節理張開與否的參量，

在整體坐標系內，劈裂單元的剛度矩陣可表示為

式中，Q為轉換矩陣，

其中θ為節理與坐標軸的夾角.

2 節理力學模型

由于原始節理(即巖石等脆性材料的已有裂隙)間含有填充物，使得節理面之間處于一種弱膠結狀態，從而具有一定的黏聚力，如圖2所示.節理面何時產生相對滑動取決于節理面之間的黏結強度和作用于節理面的法向和切向應力.在同一節理面內，膠結強度越大，節理面越不容易產生相對滑動;同時，法向應力越大，節理面之間的摩擦應力也越大，就越不容易產生相對滑動.而一旦產生相對滑動， 節理面之間的相對摩擦阻力就會突然降低.為了描述這一過程，采用如下節理面相對滑動準則，即

式中，Fs，Fn分別為節理面切向應力和法向應力，Cj為原始節理面內黏聚物黏聚力，μ為節理接觸面的摩擦系數.式(4)實質上是摩爾-庫侖準則的簡單描述形式.

圖2 原始節理示意圖Fig.2 Illustration of intact joint

在節理面滿足滑動準則之前，節理面的切向剛度與法向剛度取相同量級;當節理面滿足滑動準則時，節理面法向剛度保持不變，而切向剛度降為0 (或很小量級).因而，原始節理的三節點節理單元剛度矩陣可表示為

式中，

3 節理擴展過程的模擬

節理擴展過程就是節理尖端材料劈裂過程(見圖3)，在節理尖端取一微元，微元破壞機理可由摩爾-庫侖準則來描述.當微元應力狀態滿足摩爾-庫侖準則時，微元發生破壞，形成一個裂紋微段，表現在數值模型上，則為裂尖單元發生劈裂.因而，可以將摩爾-庫侖準則作為單元劈裂準則，即

式中，τf為巖體破壞面的抗剪強度，c為巖體黏聚力，σn為垂直破壞面的法向應力，φ為內摩擦角.

本研究假設單元劈裂面垂直于最大主應變ε1方向.同時，為了確定裂紋面的位置，假設滑移面經過三角單元的中心，如圖3所示.

4 不同圍壓下節理擴展數值模擬

4.1 數值模擬參數

為檢驗本方法的有效性，本研究應用自編的有限單元程序對圍壓條件下雙裂紋擴展進行數值模擬.試件尺寸、裂紋分布模式及試驗結果均取自文獻［6］.試件數值模擬參數如表1所示，尺寸及裂紋分布如圖4所示，其中假設巖石破壞前為線彈性.在模擬過程中，采用三節點常應變三角形單元，單元總數為28 282，節點總數為14 400;采用位移控制加載方案，每一荷載步為 0.000 6εt×635=0.160 02× 10－3mm.針對不同的巖橋傾角，每種情況采用0.35，0.70，1.50 MPa 3種不同的圍壓進行數值模擬.

4.2 0.35 MPa圍壓下的裂紋擴展分析

圖5為傾角β/α=45°/0°時不同荷載步下的裂紋擴展情況.圖中可以看出，隨著荷載步的增加，試件的外翼裂紋和內翼裂紋從初始裂紋的尖端開始豎直向上擴展;緊接著是次裂紋從初始裂紋尖端順著初始裂紋方向擴展，最終在巖橋中間位置匯合，貫穿整個巖橋.所模擬的整個裂紋擴展過程與文獻［6］中的試驗結果基本吻合.

圖4 含有預置裂紋的試件［6］Fig.4 Specimen containing pre-existed cracks［6］

圖5 β/α=45°/0°時不同荷載步下的裂紋擴展Fig.5 Fracture propagation at different loading steps when β/α=45°/0°

圖6為傾角β/α=45°/30°時不同荷載步下的裂紋擴展情況.圖中可以看出，隨著荷載步的增加，初始裂紋的2個尖端均出現了向上擴展的翼裂紋.當巖橋區出現次裂紋擴展時，張拉型的內翼裂紋與剪切型的次裂紋發生匯合，連接了初始裂紋的2個內部尖端.由此可見，巖橋區的破壞是由剪切和張拉應力共同作用的結果.

圖6 β/α=45°/30°時不同荷載步下的裂紋擴展Fig.6 Fracture propagation at different loading steps when β/α=45°/30°

圖7為傾角β/α=45°/45°時不同荷載步下的裂紋擴展情況.同其他傾角情況相似，首先是外翼裂紋向豎直方向擴展，達到一定程度以后，初始裂紋內部尖端的巖橋區出現了內翼裂紋擴展.當內翼匯合后，外翼裂紋繼續沿著豎直方向擴展，直至試件最終發生破壞.圖8～圖10分別為傾角β/α=45°/60°，45°/75°，45°/90°時不同荷載步下的裂紋擴展情況.由圖可見，隨著荷載步逐漸增加，翼裂紋仍然是最先沿荷載方向擴展，隨后在巖橋區產生剪切型裂紋與張拉型裂紋的匯合，從而導致試件最終發生破壞.

圖7 β/α=45°/45°時不同荷載步下的裂紋擴展Fig.7 Fracture propagation at different loading steps when β/α=45°/45°

圖8 β/α=45°/60°時不同荷載步下的裂紋擴展Fig.8 Fracture propagation at different loading steps when β/α=45°/60°

圖9 β/α=45°/75°時不同荷載步下的裂紋擴展Fig.9 Fracture propagation at different loading steps when β/α=45°/75°

為了便于將試驗結果與模擬結果作對比分析，現將不同巖橋傾角的模擬結果列于圖11.從圖中可以看出，當巖橋傾角大于0°時，巖橋出現了剪切裂紋匯合(即原有裂紋的內側裂尖順著原有裂紋方向產生剪切裂紋匯合)和張拉裂紋匯合(即在巖橋區內產生豎直方向的裂紋匯合).總體上看，試驗所得的裂紋擴展模式和數值模擬所得的裂紋擴展模式基本一致.

圖10 β/α=45°/90°時不同荷載步下的裂紋擴展Fig.10 Fracture propagation at different loading steps when β/α=45°/90°

但從圖11中也可以看出，數值模擬結果與試驗結果存在一些差異.如在試驗中試件會同時出現幾條翼裂紋，而且有些是在裂紋中部開始擴展，而數值模擬結果中沒有同時出現幾條翼裂紋的情況，這主要是由于本模型中沒有考慮巖石的非均質性問題.眾所周知，巖石具有非均質特性，因而在受載過程中，裂紋會在薄弱點成核、擴展等，宏觀上會出現隨機裂紋現象.而在本模型中，巖石為均質材料，因而沒有出現試驗中的隨機次生裂紋.但數值模擬結果能夠基本再現裂紋擴展模式，與試驗結果在整個趨勢上是一致的.因而，應用包含摩爾-庫侖準則的單元劈裂法可以模擬圍壓狀態下的裂紋擴展和匯合行為.

4.3 不同圍壓下裂紋擴展對比分析

為分析不同圍壓下的裂紋擴展特點，本研究針對不同的巖橋傾角，分別施加0.35，0.70，1.50 MPa 3種不同圍壓對裂紋擴展進行模擬.圖12為當β/α=45°/45°時，3種不同圍壓下的裂紋擴展模式.從圖中可以看出，在圍壓由 0.35 MPa增大到0.70 MPa的過程中，裂紋擴展模式基本沒有變化，只是在0.70 MPa圍壓下，張拉型翼裂紋擴展變“慢”了.但當圍壓增大到1.50 MPa時，張拉型翼裂紋不存在了，取而代之的是剪切型裂紋從內側裂尖開始擴展.這是因為在低圍壓狀態下，圍壓還不足以約束張拉型翼裂紋擴展;但隨著圍壓逐漸增大，這種約束力越來越強，以至于張拉型翼裂紋擴展變“慢”;最后，隨著圍壓的進一步增大，圍壓完全約束住了張拉型翼裂紋擴展，從而出現剪切裂紋擴展，這一趨勢符合圍壓條件下的裂紋擴展特征.圖13為當β/α= 45°/45°時，不同圍壓下的軸向應力-應變曲線.從圖中可以看出，隨著圍壓的增大，應力曲線的峰值也在增大，這與實際情況相符.

圖11 裂紋匯合模式的試驗結果［6］與模擬結果對比Fig.11 Comparison of the coalescence pattern between the experimental［6］and the simulated results

圖12 不同圍壓下裂紋擴展對比(β/α=45°/45°)Fig.12 Comparison offracturepropagation under different confining stresses(β/α=45°/45°)

5 結束語

通過與現有雙裂紋圍壓試驗結果進行對比可以發現，包含摩爾-庫侖準則的單元劈裂法可以較好地模擬出巖體節理的擴展和匯合過程，能基本再現巖體試件的破環特征.這是由于單元劈裂法能夠在網格不變的情況下模擬節理的擴展和匯合行為，因此避免了重新劃分網格，減少了計算量，從而在數值模擬節理巖體破壞的過程中具有一定的優勢.但單元劈裂法沒有考慮劈裂后所形成的2個塊體自身的彈塑性變形，所以該方法只是一種近似方法.

圖13 不同圍壓下的軸向應力-應變曲線(β/α=45°/45°)Fig.13 Uniaxial stress-strain curves under different confining stresses(β/α=45°/45°)

［1］ YANGS Q，DAIY H，HANL J，et al.Experimental study on mechanical behavior of brittle marble samples containing different flaws under uniaxial compression［J］.Engineering Fracture Mechanics，2009，76：1833-1845.

［2］ 欒茂田，張大林，楊慶，等.有限覆蓋無單元法在裂紋擴展數值分析問題中的應用［J］.巖土工程學報，2003，25(5)：527-531.

［3］ 車法星，黎立云，劉大安.類巖材料多裂紋體斷裂破壞試驗及有限元分析［J］.巖石力學與工程學報，2000，19(3)：295-298.

［4］ 黃明利，馮夏庭，王水林.多裂紋在不同巖石介質中的擴展貫通機制分析［J］.巖土力學，2002，23(2)：142-146.

［5］ GOODMANR E，TAYLORR L，BREKKET L.A model for the mechanics of joint rock［J］.Journal of the Soil Mechanics and Foundations Division，1968，94：637-659.

［6］ MUGHIEDAO，KARASNEHI.Coalescence of offset rock joints under biaxial loading［J］.Geotechnical and Geological Engineering，2006，24：985-999.

［7］ 張振南，陳永泉.一種模擬節理巖體破壞的新方法：單元劈裂法［J］.巖土工程學報，2009，31(12)：1858-1864.

［8］ ZHANGZ N，CHEN Y Q.Simulation of fracture propagation subjected to compressive and shear stress field using virtual multidimensional internal bonds［J］.International Journal of Rock Mechanics and Mining Sciences，2009，46(6)：1010-1022.

［9］ MOESN，DOLBOWJ，BELYTSCHKOT.A finite element method for crack growth without remeshing［J］.Int J Numer Methods Engrg，1999，46(1)：131-150.

［10］ LARSSONR，FAGERSTROMM.A framework for fracture modelling based on the material forces concept with XFEM kinematics［J］.Int J Numer Methods Engrg，2005，62(13)：1763-1788.