基于FLAC3D二次開發的含預制缺陷巖體數值模擬

周廷強，李江華

（1．重慶三峽學院 土木工程學院，重慶404100；2．河南工業職業技術學院 城市建設學院，河南 南陽473000）

經過漫長的地質作用以及各種工程擾動，巖體內部存在著各種節理、裂隙、孔隙、空洞等缺陷，在外部載荷的作用下，巖體的變形與破壞明顯受控于巖體中的缺陷，呈現出明顯的結構特征。同時缺陷的存在使得巖體的強度特征和變形破壞規律變得更加復雜，研究含缺陷巖體的失穩斷裂規律對于保障缺陷結構工程的穩定性具有重要的理論價值和實踐意義[1]。對此，國內外學者從室內試驗到數值模擬做了大量的研究工作，其中通過室內試驗還原缺陷巖體的結構特征、地質條件及工程力學環境模擬巖體的變形破壞過程和機制已經得到廣泛應用。無論是單條[2]、雙條[3]還是多條[4]預制缺陷的單軸抗壓[5]、直剪[6]以及不同圍壓[7]條件下的室內試驗研究均已得到開展。但室內試驗建立的含缺陷試件比較復雜，且不可控因素較多，取樣技術難，成本高，重復性差，這些都成為研究的瓶頸，且工程尺度巖體往往存在尺寸效應，室內試驗尺度不能直接應用于工程尺度巖體，利用數值方法研究節理巖體力學性質成為目前很普遍的手段。Reyes 等[8]運用有限元方法研究了含2 條預制缺陷巖體試樣的起裂、擴展和貫通模式，并利用物理試驗結果對數值試驗結果正確性進行了初步驗證。指出適當的選用數值模型可以得到與模型試驗相同的裂紋的起裂、擴展以及貫通模式。Bobet[9]等利用FROCK 程序建立了1 個適用于拉、剪復合破壞模式的準則的模型，并且研究了預制缺陷與外荷載形成不同角度時裂紋尖端的應力場。唐春安[10]基于細觀損傷理論自行開發的RFPA數值計算軟件，在細觀層次上對巖石、混凝土等非均勻脆性材料的裂紋擴展規律進行了系統研究。梁正召等[11]采用細觀損傷數值模擬方法模擬了單軸壓縮條件下巖石試樣的破壞過程，數值模擬得到了表面裂隙內部擴展、貫通過程，動態再現翼型裂紋、殼體裂紋的形態，探討三維裂紋內部的受力機制，推測可能發生的斷裂類型。上述研究從實驗和數值模擬角度完善了缺陷巖體中強度變形規律以及裂紋擴展貫通原理，取得了積極有益的進展。但上述的試驗與模擬研究中很少涉及到工程尺度上缺陷巖體，限制了研究成果的應用。為此提出了1 種彈脆性損傷本構模型，并利用C++語言進行二次開發將該模型嵌入到FLAC3D軟件中，該方法不僅可以表征裂紋擴展規律，也可以得到裂紋擴展時微觀上的拉剪破裂類型，而且基于FLAC3D軟件的二次開發使得該方法適用于大變形下的巖石力學問題。作為案例分析，利用該模型對含不同傾角的雙缺陷巖體進行數值模擬得到其強度及微裂紋擴展演化規律。

1 數值計算模型

1.1 庫倫破壞準則與彈脆性本構關系

假定巖石的破壞主要是剪切破壞，巖石的強度，即抗摩擦強度等于巖石本身抗剪切摩擦的黏聚力和剪切面上法向應力產生的摩擦力。即平面上的剪切強度準則為：

式中：τ 為剪切面上的剪應力（剪切強度）；σ 為剪切面上的正應力；c 為黏聚力；φ 為內摩擦角。

在此描述的基礎上，可以討論巖石的破裂條件及其加載應力場的關系。在平面應力狀態下，庫倫準則可以用莫爾極限圓直觀的圖解描述，Mohr-Coulomb 破壞準則如圖1。

確切的準則由強度曲線AL 表示，其斜率為f=tanφ，在τ 軸上截距為c。平面應力狀態下應力σ 與τ由最主應力σ1與最小主應力σ3確定的應力圓所確定。

圖1 Mohr-Coulomb 破壞準則

如果應力圓上的點落在強度曲線AL 之下，則說明該點表示的應力還沒有達到材料的強度值，材料不發生破壞；如果應力圓上的點落在強度曲線AL之上，則說明該點的應力已超過材料的強度并發生破壞；如果應力圓上的點正好與強度曲線AL 相切（圖1 中D 點），則說明材料處于極限平衡狀態，材料所產生的剪切破壞將可能在該點所對應的剪切面上發生。

事實上，細觀單元的破裂不僅僅只有剪切破壞，將細觀單元破裂分為2 種，即拉伸破裂和剪切破裂。當應變狀態達到了最大拉伸應變準則引起拉伸破裂；壓應力或者剪應力導致應力狀態滿足剪切破壞準則，則發生剪切破裂。帶拉伸破壞準則的庫倫準則的表達式如下：

式中：σc為巖石的單軸抗壓強度；σt為巖石的單軸抗拉強度

假定巖石單元在細觀上的本構關系服從彈-脆性損傷模型（圖2），在階段I，巖石單元具有線彈性的簡單應力-應變特性，該階段為彈性段。在階段II，巖石單元達到一定的強度準則后發生破裂轉換為弱介質，則時候巖石單元具有殘余強度特性。圖中的εc和εr分別代表巖石單元達到強度準則時對應的應變，σcr與σtr分別代表巖石單元破裂后殘余抗壓與抗拉強度。

1.2 FLAC3D自定義本構開發

FLAC3D提供了1 種相對便利的開發環境，允許用戶按照自己的意愿開發自定義本構。該本構利用C++開發后編譯成DLL（動態鏈接庫）文件，用戶在需要使用該本構時只需調用即可。自定義本構的主要功能是返回產生相應應變增量后的新應力值，同時該模型必須提供其他的信息，如模型的名字及該模型中用到的材料參數，同時需要制定模型同軟件交互的方式。文獻中[12-13]已經提及基于C++的FLAC3D自定義本構開發的方法與流程，這里只對其2 個關鍵進行介紹。

圖2 彈-脆性損傷本構關系

1）本構模型的注冊。每1 個用戶自定義本構模型被編譯進1 個DLL 文件，該DLL 文件在FLAC3D中必須被實例化。按照慣例，在被用作FLAC3D插件的DLL 文件中包含4 個輸出函數：getName、getMajorVersion、getMinorVersion 和creatInstance。此外還需提供1 個被稱為DLLMain 的存根函數，當函數庫從系統中加載或寫在是會被調用。其中DLLMain 文件始終是相同的；輸出的getName 始終返回1 個以“model”開頭的字符串，且該字符串必須是獨立的，以同其他本構模型區分開來。getMajorVersion 函數不應當被改變，表示本構模型的主要版本；getMinorVersion 表示的本構模型的小版本；creatInstance 函數實際上創建并返，1 個類別實體，它被存儲于注冊表中，用以創建其他的實體。

2）循環過程中本構與FLAC3D間信息的傳遞。FLAC3D與用戶自定義本構間最重要的聯系是被稱為run 的成員函數，它計算出循環過程中模型的力學響應。1 個被稱為State 的構架被用于與模型進行信息交流。成員函數run 的主要任務是計算新應變增量對應的新應力。在非線性模型中，與模型的內部狀態交流信息也是非常重要的，以方便狀態的輸出。開發的彈脆性損傷模型需要輸入的參數有彈性模量E，泊松比μ，密度ρ，黏聚力c，內摩擦角φ，抗壓強度σc，抗拉強度σt，殘余抗壓強度σcr，殘余抗拉強度σtr共9 個參數。

2 含預制缺陷脆性材料數值模擬

2.1 數值計算模型

利用彈脆性損傷本構模型（圖2），研究不同缺陷角度對巖石的強度特性、應力變化與裂紋擴展規律及破壞模式的影響。數值模擬的模型示意圖如圖3，其中模型為50 mm×50 mm×100 mm 的長方體試件；缺陷角度為缺陷面與水平面的夾角，分別為0°、15°、30°、45°、60°、75°、90°共7 中工況；缺陷長度為預制缺陷與試樣自由表面間的交線長度，為20 mm；裂紋厚度為預制缺陷與試樣自由表面間的交線厚度，為2 mm；巖橋角度為兩缺陷間最短連線與水平面的夾角，為45°；巖橋長度為兩缺陷間最短連線的長度，為14 mm。

圖3 數值計算模型

利用文獻[14]中提及的網格劃分方法，采用TCL腳本語言對SURPAC 進行二次開發，利用FLAC3D模型構建的數據轉換方法，能夠有效避免產生重復節點，從而實現FLAC3D模型的快速構建。數值計算網格如圖4，劃分網格均為1 mm×1 mm×1 mm 的正方體網格。

圖4 數值計算網格

模擬砂巖選取力學參數[15]：彈性模量為2.0 GPa，泊松比2.5，密度2 500 kg/m3，黏聚力2.5 MPa，內摩擦角46°，單軸抗壓強度為20 MPa，單軸抗拉強度為1.0 MPa，殘余抗壓強度為2.0 MPa，殘余抗拉強度為0.1 MPa。采用兩端位移加載方式，上下兩端加載速率均為5×10-5mm/時步[16]。

2.2 數值計算結果

力學參數隨缺陷角度的變化關系如圖5，缺陷角度0°、15°、30°、45°、60°、75°、90°試件的峰值強度分別為13.75、13.45、12.74、11.23、9.45、9.73、12.27 MPa，很明顯的可以看出隨著傾角的增大，峰值強度成非線性變化趨勢，從0°到60°,隨著缺陷傾角的增大，峰值強度逐漸降低，從60°到90°,隨著缺陷傾角的增大，峰值強度逐漸增加。于之相對應的彈性模量表現出與峰值強度類似的變化規律，隨著缺陷傾角的增大，彈性模量先減小后增大，彈性模量分別為1.54、1.48、1.43、1.24、1.16、1.23、1.47 GPa，傾角 為60°的彈性模量最小。

圖5 力學參數隨缺陷角度的變化關系

最大主應力隨缺陷角度的變化關系圖如圖6。當缺陷傾角為0°和15°時，首先在預制缺陷的兩側與巖橋區域產生一定的應力集中，巖橋迅速貫通，同時在缺陷尖端萌生拉應力引起的翼型裂紋，翼型裂紋沿著加載方向延伸一定程度不再擴展，在預制缺陷尖端產生新的次生裂紋，最終壓剪裂紋貫通使試件發生破壞。缺陷傾角為30°和45°時破壞模式相似，2 條預制缺陷均以翼型裂紋形式擴展，應力集中出現在巖橋區域，使巖橋區域發生剪切線貫通破壞，兩外端的裂紋繼續沿加載方向擴展，最終由翼型裂紋貫通試件引起整體失穩破壞。缺陷傾角為60°和75°時，在預制缺陷尖端同時產生翼型和反翼型裂紋，巖橋區域出現1 條剪切裂紋，與兩預制裂隙尖端的裂紋貫通使得試件發生宏觀破壞。缺陷傾角為90°時，巖橋貫通不明顯，最終由預制缺陷與試件上下表面貫通引發宏觀破壞。

圖6 最大主應力隨缺陷角度的變化關系圖

對試件破壞過程中的微破裂類型進行分析，塑性區隨缺陷角度的變化關系如圖7。

圖7 塑性區隨缺陷角度的變化關系

當缺陷角度為0°和15°時，拉伸破壞在整個破壞模式中相對其他角度有更大的占比，約40%，主要集中在缺陷的上下部分，該部分只有少量的剪切破壞摻雜。2 條缺陷之間區域為剪切貫通，拉破壞較少。當缺陷角度為30°和45°時，拉伸破壞相對于0°和15°時較為集中，主要分布于2 條缺陷上下側，缺陷中間區域以剪切破壞為主，但相對于0°和15°時，缺陷內側尖端區域以拉伸破壞為主，缺陷外側尖端以混合型破壞和剪切破壞為主。整體上看，拉伸破壞占比相對0°和15°時有些微下降，約20%，和剪切破壞占比均有所增加。當缺陷角度為60°時，拉伸破壞區急劇縮減到5%，主要集中在缺陷內測靠尖端區域，整個試件明顯開始以剪切破壞為主。當缺陷角度為75°和90°，拉破壞占比下降到1%以內，試件的破壞模式轉變為剪切破壞。

3 結 語

1）通過C++平臺FLAC3D進行二次開發，開發了1 種彈脆性損傷本構模型，該模型可以很好的反應巖石單元細觀破壞。對其彈脆性本構進行詳細描述，同時也展示FLAC3D二次開發的基本過程。

2）通過自主開發的FLAC3D彈脆性本構模型，對含不同傾角的預制缺陷巖體進行數值模擬計算，得到其強度隨缺陷角度的增加先下降后上升的規律，同時其細觀破裂由低缺陷角度時的拉伸-剪切型的混合破壞轉變為高缺陷角度時的剪切破壞。