基于FLAC3D二次開發(fā)的含預(yù)制缺陷巖體數(shù)值模擬

周廷強(qiáng)，李江華

（1．重慶三峽學(xué)院 土木工程學(xué)院，重慶404100；2．河南工業(yè)職業(yè)技術(shù)學(xué)院 城市建設(shè)學(xué)院，河南 南陽473000）

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

1 數(shù)值計(jì)算模型

1.1 庫倫破壞準(zhǔn)則與彈脆性本構(gòu)關(guān)系

假定巖石的破壞主要是剪切破壞，巖石的強(qiáng)度，即抗摩擦強(qiáng)度等于巖石本身抗剪切摩擦的黏聚力和剪切面上法向應(yīng)力產(chǎn)生的摩擦力。即平面上的剪切強(qiáng)度準(zhǔn)則為：

式中：τ 為剪切面上的剪應(yīng)力（剪切強(qiáng)度）；σ 為剪切面上的正應(yīng)力；c 為黏聚力；φ 為內(nèi)摩擦角。

在此描述的基礎(chǔ)上，可以討論巖石的破裂條件及其加載應(yīng)力場的關(guān)系。在平面應(yīng)力狀態(tài)下，庫倫準(zhǔn)則可以用莫爾極限圓直觀的圖解描述，Mohr-Coulomb 破壞準(zhǔn)則如圖1。

確切的準(zhǔn)則由強(qiáng)度曲線AL 表示，其斜率為f=tanφ，在τ 軸上截距為c。平面應(yīng)力狀態(tài)下應(yīng)力σ 與τ由最主應(yīng)力σ1與最小主應(yīng)力σ3確定的應(yīng)力圓所確定。

圖1 Mohr-Coulomb 破壞準(zhǔn)則

如果應(yīng)力圓上的點(diǎn)落在強(qiáng)度曲線AL 之下，則說明該點(diǎn)表示的應(yīng)力還沒有達(dá)到材料的強(qiáng)度值，材料不發(fā)生破壞；如果應(yīng)力圓上的點(diǎn)落在強(qiáng)度曲線AL之上，則說明該點(diǎn)的應(yīng)力已超過材料的強(qiáng)度并發(fā)生破壞；如果應(yīng)力圓上的點(diǎn)正好與強(qiáng)度曲線AL 相切（圖1 中D 點(diǎn)），則說明材料處于極限平衡狀態(tài)，材料所產(chǎn)生的剪切破壞將可能在該點(diǎn)所對(duì)應(yīng)的剪切面上發(fā)生。

事實(shí)上，細(xì)觀單元的破裂不僅僅只有剪切破壞，將細(xì)觀單元破裂分為2 種，即拉伸破裂和剪切破裂。當(dāng)應(yīng)變狀態(tài)達(dá)到了最大拉伸應(yīng)變準(zhǔn)則引起拉伸破裂；壓應(yīng)力或者剪應(yīng)力導(dǎo)致應(yīng)力狀態(tài)滿足剪切破壞準(zhǔn)則，則發(fā)生剪切破裂。帶拉伸破壞準(zhǔn)則的庫倫準(zhǔn)則的表達(dá)式如下：

式中：σc為巖石的單軸抗壓強(qiáng)度；σt為巖石的單軸抗拉強(qiáng)度

假定巖石單元在細(xì)觀上的本構(gòu)關(guān)系服從彈-脆性損傷模型（圖2），在階段I，巖石單元具有線彈性的簡單應(yīng)力-應(yīng)變特性，該階段為彈性段。在階段II，巖石單元達(dá)到一定的強(qiáng)度準(zhǔn)則后發(fā)生破裂轉(zhuǎn)換為弱介質(zhì)，則時(shí)候巖石單元具有殘余強(qiáng)度特性。圖中的εc和εr分別代表巖石單元達(dá)到強(qiáng)度準(zhǔn)則時(shí)對(duì)應(yīng)的應(yīng)變，σcr與σtr分別代表巖石單元破裂后殘余抗壓與抗拉強(qiáng)度。

1.2 FLAC3D自定義本構(gòu)開發(fā)

FLAC3D提供了1 種相對(duì)便利的開發(fā)環(huán)境，允許用戶按照自己的意愿開發(fā)自定義本構(gòu)。該本構(gòu)利用C++開發(fā)后編譯成DLL（動(dòng)態(tài)鏈接庫）文件，用戶在需要使用該本構(gòu)時(shí)只需調(diào)用即可。自定義本構(gòu)的主要功能是返回產(chǎn)生相應(yīng)應(yīng)變?cè)隽亢蟮男聭?yīng)力值，同時(shí)該模型必須提供其他的信息，如模型的名字及該模型中用到的材料參數(shù)，同時(shí)需要制定模型同軟件交互的方式。文獻(xiàn)中[12-13]已經(jīng)提及基于C++的FLAC3D自定義本構(gòu)開發(fā)的方法與流程，這里只對(duì)其2 個(gè)關(guān)鍵進(jìn)行介紹。

圖2 彈-脆性損傷本構(gòu)關(guān)系

1）本構(gòu)模型的注冊(cè)。每1 個(gè)用戶自定義本構(gòu)模型被編譯進(jìn)1 個(gè)DLL 文件，該DLL 文件在FLAC3D中必須被實(shí)例化。按照慣例，在被用作FLAC3D插件的DLL 文件中包含4 個(gè)輸出函數(shù)：getName、getMajorVersion、getMinorVersion 和creatInstance。此外還需提供1 個(gè)被稱為DLLMain 的存根函數(shù)，當(dāng)函數(shù)庫從系統(tǒng)中加載或?qū)懺谑菚?huì)被調(diào)用。其中DLLMain 文件始終是相同的；輸出的getName 始終返回1 個(gè)以“model”開頭的字符串，且該字符串必須是獨(dú)立的，以同其他本構(gòu)模型區(qū)分開來。getMajorVersion 函數(shù)不應(yīng)當(dāng)被改變，表示本構(gòu)模型的主要版本；getMinorVersion 表示的本構(gòu)模型的小版本；creatInstance 函數(shù)實(shí)際上創(chuàng)建并返，1 個(gè)類別實(shí)體，它被存儲(chǔ)于注冊(cè)表中，用以創(chuàng)建其他的實(shí)體。

2）循環(huán)過程中本構(gòu)與FLAC3D間信息的傳遞。FLAC3D與用戶自定義本構(gòu)間最重要的聯(lián)系是被稱為run 的成員函數(shù)，它計(jì)算出循環(huán)過程中模型的力學(xué)響應(yīng)。1 個(gè)被稱為State 的構(gòu)架被用于與模型進(jìn)行信息交流。成員函數(shù)run 的主要任務(wù)是計(jì)算新應(yīng)變?cè)隽繉?duì)應(yīng)的新應(yīng)力。在非線性模型中，與模型的內(nèi)部狀態(tài)交流信息也是非常重要的，以方便狀態(tài)的輸出。開發(fā)的彈脆性損傷模型需要輸入的參數(shù)有彈性模量E，泊松比μ，密度ρ，黏聚力c，內(nèi)摩擦角φ，抗壓強(qiáng)度σc，抗拉強(qiáng)度σt，殘余抗壓強(qiáng)度σcr，殘余抗拉強(qiáng)度σtr共9 個(gè)參數(shù)。

2 含預(yù)制缺陷脆性材料數(shù)值模擬

2.1 數(shù)值計(jì)算模型

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

圖3 數(shù)值計(jì)算模型

利用文獻(xiàn)[14]中提及的網(wǎng)格劃分方法，采用TCL腳本語言對(duì)SURPAC 進(jìn)行二次開發(fā)，利用FLAC3D模型構(gòu)建的數(shù)據(jù)轉(zhuǎn)換方法，能夠有效避免產(chǎn)生重復(fù)節(jié)點(diǎn)，從而實(shí)現(xiàn)FLAC3D模型的快速構(gòu)建。數(shù)值計(jì)算網(wǎng)格如圖4，劃分網(wǎng)格均為1 mm×1 mm×1 mm 的正方體網(wǎng)格。

圖4 數(shù)值計(jì)算網(wǎng)格

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

2.2 數(shù)值計(jì)算結(jié)果

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

圖5 力學(xué)參數(shù)隨缺陷角度的變化關(guān)系

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

圖6 最大主應(yīng)力隨缺陷角度的變化關(guān)系圖

對(duì)試件破壞過程中的微破裂類型進(jìn)行分析，塑性區(qū)隨缺陷角度的變化關(guān)系如圖7。

圖7 塑性區(qū)隨缺陷角度的變化關(guān)系

當(dāng)缺陷角度為0°和15°時(shí)，拉伸破壞在整個(gè)破壞模式中相對(duì)其他角度有更大的占比，約40%，主要集中在缺陷的上下部分，該部分只有少量的剪切破壞摻雜。2 條缺陷之間區(qū)域?yàn)榧羟胸炌ǎ茐妮^少。當(dāng)缺陷角度為30°和45°時(shí)，拉伸破壞相對(duì)于0°和15°時(shí)較為集中，主要分布于2 條缺陷上下側(cè)，缺陷中間區(qū)域以剪切破壞為主，但相對(duì)于0°和15°時(shí)，缺陷內(nèi)側(cè)尖端區(qū)域以拉伸破壞為主，缺陷外側(cè)尖端以混合型破壞和剪切破壞為主。整體上看，拉伸破壞占比相對(duì)0°和15°時(shí)有些微下降，約20%，和剪切破壞占比均有所增加。當(dāng)缺陷角度為60°時(shí)，拉伸破壞區(qū)急劇縮減到5%，主要集中在缺陷內(nèi)測靠尖端區(qū)域，整個(gè)試件明顯開始以剪切破壞為主。當(dāng)缺陷角度為75°和90°，拉破壞占比下降到1%以內(nèi)，試件的破壞模式轉(zhuǎn)變?yōu)榧羟衅茐摹?/p>

3 結(jié) 語

1）通過C++平臺(tái)FLAC3D進(jìn)行二次開發(fā)，開發(fā)了1 種彈脆性損傷本構(gòu)模型，該模型可以很好的反應(yīng)巖石單元細(xì)觀破壞。對(duì)其彈脆性本構(gòu)進(jìn)行詳細(xì)描述，同時(shí)也展示FLAC3D二次開發(fā)的基本過程。

2）通過自主開發(fā)的FLAC3D彈脆性本構(gòu)模型，對(duì)含不同傾角的預(yù)制缺陷巖體進(jìn)行數(shù)值模擬計(jì)算，得到其強(qiáng)度隨缺陷角度的增加先下降后上升的規(guī)律，同時(shí)其細(xì)觀破裂由低缺陷角度時(shí)的拉伸-剪切型的混合破壞轉(zhuǎn)變?yōu)楦呷毕萁嵌葧r(shí)的剪切破壞。