考慮初始滲透壓力的裂隙巖體本構(gòu)模型二次開發(fā)及其驗證

張衛(wèi)東， 陳士海

(華僑大學(xué) 土木工程學(xué)院， 福建 廈門 361021)

考慮初始滲透壓力的裂隙巖體本構(gòu)模型二次開發(fā)及其驗證

張衛(wèi)東， 陳士海

(華僑大學(xué) 土木工程學(xué)院， 福建 廈門 361021)

為充分研究含水裂隙巖體的力學(xué)特性，在摩爾庫倫模型的基礎(chǔ)上，充分考慮初始滲透壓力的作用.對本構(gòu)方程中的柔度張量進行適當(dāng)?shù)男拚?，生成考慮初始滲透壓力的裂隙巖體本構(gòu)模型.為了增加對比性，按照同樣的理論推導(dǎo)方法及二次開發(fā)流程，在不考慮初始滲透壓力的影響下，成功地生成未考慮初始滲透壓力的裂隙巖體本構(gòu)模型.在此基礎(chǔ)上，利用C++語言，將以上兩種模型編譯成可運行的動態(tài)鏈接庫文件.為驗證自定義本構(gòu)模型的準(zhǔn)確性，在FLAC3D中建立同一個50 m×50 m×1 m的圓形隧道數(shù)值模型進行對比分析.通過分別求解計算3種本構(gòu)模型,分析各圓形隧道上4個對稱點的位移變化情況.結(jié)果表明:將初始滲透壓力直接引進本構(gòu)方程中，利用C++二次開發(fā)出來的本構(gòu)模型是成功的.

裂隙巖體; 初始滲透壓力； 柔度張量； 本構(gòu)關(guān)系； 二次開發(fā)

自然界中的巖體是一種含有初始損傷的介質(zhì)，它包含著各式各樣的非連續(xù)面，如斷層、節(jié)理、裂隙，而含有初始地下水的裂隙巖體在工程的施工中會經(jīng)常碰到.1999年，朱維申等[1]在Betti能量互易定理和不可逆熱力學(xué)定律的基礎(chǔ)之上，建立了裂隙巖體的三維脆彈性斷裂損傷本構(gòu)模型.2008年，藍航等[2]在幾何損傷力學(xué)理論的基礎(chǔ)上，建立了節(jié)理巖體的采動損傷本構(gòu)模型.2012年，齊銀萍[3]在理論分析本構(gòu)模型的基礎(chǔ)之上，在FLAC3D中編寫出了損傷流變本構(gòu)方程.2014年，文獻[4-5]等利用FLAC3D建立了含水裂隙巖體的數(shù)值模型.但是，以上研究開發(fā)出來的模型大部分是針對無水裂隙巖體，而對于有水裂隙巖體，則是利用模擬軟件中相應(yīng)的滲流模塊進行流固耦合計算，而在本構(gòu)模型當(dāng)中沒有將滲透壓力帶入到本構(gòu)方程中進行二次開發(fā).因此，初始含水裂隙巖體工程的模擬分析存在一定的局限性.本文在摩爾庫倫模型的基礎(chǔ)之上，充分考慮初始滲透壓力的作用，并將其直接引進本構(gòu)方程.

1 含水裂隙巖體的本構(gòu)關(guān)系

從宏觀力學(xué)效果方面考慮，裂隙巖體的初始損傷及其損傷的演化可用柔度張量的改變進行表示，而裂隙的半徑、法向和切向剛度系數(shù)、傳壓和傳剪系數(shù)、初始滲透壓力等都會在不同程度上對柔度張量產(chǎn)生影響.文獻[6]通過坐標(biāo)運算和疊加原理得到了含水裂隙巖體柔度張量的表達式，即

1.1 含水裂隙巖體的初始損傷柔度張量

(2)

式(2)中：KⅠ,KⅡ,KⅢ分別為裂隙尖端的Ⅰ，Ⅱ，Ⅲ型應(yīng)力強度因子；Ω為裂隙的外表面，積分沿整個裂隙表面進行.

現(xiàn)假定裂隙都為圓形，且半徑為a，在查閱應(yīng)力強度因子手冊后，可獲得裂隙尖端的Ⅰ，Ⅱ，Ⅲ型應(yīng)力強度因子[8].

(3)

式(3)中：σ，τ分別為應(yīng)力張量在裂隙面法向和切向上的投影分量；G1，G2為裂隙形狀系數(shù).

在某一應(yīng)力狀態(tài)下，裂隙面上實際的法向和切向應(yīng)力張量是需要進行修正的.為此,引進Cn和Cv,反映裂隙在某一應(yīng)力狀態(tài)下，裂隙面?zhèn)鬟f法向應(yīng)力和切向應(yīng)力的能力[9].對于含水裂隙巖體，還需要考慮初始滲透壓力(p)的作用，因此,綜合以上情況,作用在裂隙面上的有效法向和切向應(yīng)力分別為

τ′=(1-Cv)τ.

(5)

將式(4)，(5)帶入式(3)，化簡可得到單位體積巖體內(nèi)該裂隙產(chǎn)生的附加彈性應(yīng)變能φd，即

[7]中的方法，求得單位體積巖體內(nèi)裂隙產(chǎn)生的附加彈性應(yīng)變能φd.

且有

其中：ni，nj，nk，nl均為裂隙面的法向向量分量.

通過比較式(7)及文獻[7]中的表達式，可得由于初始滲透壓力存在而產(chǎn)生的附加柔度張量，即

(8)

1.2 裂隙擴展產(chǎn)生的損傷演化柔度張量

1.2.1壓剪應(yīng)力狀態(tài)下的損傷演化柔度張量 在壓剪應(yīng)力狀態(tài)下，巖體中的裂隙隨著外加作用載荷的增加，開始閉合摩擦滑動，壓剪起裂，形成分支張型裂隙，新裂面順著最大主應(yīng)力方向不斷延伸發(fā)展，直至位于分支張型裂隙尖端的微裂隙損傷區(qū)相互匯合，導(dǎo)致宏觀裂隙的擊穿貫通，而使巖體局部破壞[10].

大量的現(xiàn)場勘測結(jié)果和理論計算表明，巖體中的裂隙在壓剪應(yīng)力的作用下，將在最大和最小主壓應(yīng)力組成的平面內(nèi)沿著最大主壓應(yīng)力σ1的方向穩(wěn)定擴展[11]，其翼形分支裂隙發(fā)生擴展.

設(shè)遠場Cauchy應(yīng)力張量在裂隙面法向和切向的投影分別為正應(yīng)力σn和剪應(yīng)力τn，即

式(9)中：θ為裂隙與最大主壓應(yīng)力的夾角.

當(dāng)考慮初始滲透壓力的影響時，裂隙面上實際傳遞的法向應(yīng)力σne和切向應(yīng)力τne分別為

壓剪應(yīng)力狀態(tài)下，考慮初始滲透壓力的裂隙巖體其裂隙擴展產(chǎn)生的損傷演化柔度張量可根據(jù)應(yīng)變能等效原理和自洽理論推導(dǎo)得出[12]，壓剪應(yīng)力作用下分支裂隙尖端瞬時應(yīng)力強度因子[13]為

(11)

式(11)中：τe=τne-τs，τs為裂隙面上下表面抵抗外力的能力.

τe=fsσn+Cs表示摩爾庫倫準(zhǔn)則[14]，fs為裂隙面的摩擦系數(shù)(fs=tanφ，φ為巖體介質(zhì)的內(nèi)摩擦角)，Cs為裂隙面的黏結(jié)力，σn為裂隙面的法向應(yīng)力，且有σn=σne.

擴展中的翼形分支裂隙逐漸沿平行于最大主壓應(yīng)力的方向穩(wěn)定擴展，當(dāng)分支裂隙擴展至KI=KlC時，裂隙停止擴展[15].其中，KI為翼形分支裂隙尖端的應(yīng)力強度因子；KIC為裂隙巖體的斷裂韌度.翼形分支裂隙尖端應(yīng)力強度因子計算方法采用Kemeny計算模型[13].

1.2.2 拉剪應(yīng)力狀態(tài)下的損傷演化柔度張量 在拉剪應(yīng)力狀態(tài)下，巖體中的裂隙受拉后,出現(xiàn)張開現(xiàn)象，其裂隙表面的黏結(jié)力將變?yōu)榱?而不能傳遞拉應(yīng)力.在最大周向應(yīng)力準(zhǔn)則的基礎(chǔ)之上，當(dāng)分支裂隙尖端瞬時應(yīng)力強度因子KI0達到巖體的斷裂韌度KlC時，裂隙便開始擴展.

拉剪應(yīng)力狀態(tài)下，分支裂隙尖端瞬時應(yīng)力強度因子為

(12)

式(12)中：開裂角θ0可參考文獻[12]中的表達式.

在文獻[16]的基礎(chǔ)之上，翼形分支裂隙的擴展長度，即

可得拉剪應(yīng)力狀態(tài)下考慮初始滲透壓力的裂隙巖體其裂隙擴展產(chǎn)生的損傷演化柔度張量[13].

2 裂隙巖體自定義本構(gòu)的二次開發(fā)

利用FLAC3D的二次開發(fā)平臺，可以按照自己的意愿進行本構(gòu)模型的自定義.自定義本構(gòu)模型的計算流程圖，如圖1所示.

圖1 自定義本構(gòu)模型的計算流程圖Fig.1 Calculation flow chart of the constitutive model

3 模型的驗證及結(jié)論

在FLAC3D中,建立的統(tǒng)一圓形隧道模型在水平方向上和豎直方向上的長度都為50 m，進深為1 m，隧道半徑為5 m.模型劃分單元數(shù)為2 500，節(jié)點數(shù)為5 200.模型中：所有節(jié)點y方向上的位移被約束，z=-25 m的底面及x=±25 m的左右兩面固定，即模型底面及左右兩面上的節(jié)點約束其位移.

模型中，節(jié)點施加的初始應(yīng)力條件為

(14)

式(14)中：σx，x,σy，y,σz，z分別表示模型節(jié)點受到的x方向、y方向和z方向上的應(yīng)力.

經(jīng)測定，在模型底面滲透壓力為0.275 MPa，模型頂面為0.225 MPa，且沿z軸呈線性分布；隧道圍巖彈性模量為15 MPa，容重為25 kN·m-3，內(nèi)摩擦角為34°，粘聚力為0.85 MPa，抗拉強度為2 MPa，斷裂韌度為1.2 MPa·m-1/2；裂隙特征長度為0.1 m，內(nèi)摩擦角為11.5°，體積密度為0.33，裂隙面軸向夾角為60°，法向剛度為2.5 MPa·m-1，切向剛度為1.15 MPa.

圖2 數(shù)值計算模型Fig.2 Numerical calculation model

為驗證自定義本構(gòu)模型的準(zhǔn)確性，采用“建立同一數(shù)值模型、分別調(diào)用3種不同本構(gòu)模型(即摩爾庫倫模型、未考慮初始滲透壓力p的裂隙巖體本構(gòu)模型，以及考慮初始滲透壓力p的裂隙巖體本構(gòu)模型)”的方法進行數(shù)值模擬，并分別將以上3個數(shù)值模擬記為SZ1，SZ2和SZ3.數(shù)值計算模型,如圖2所示.

圖3 圓形隧道上的對稱節(jié)點布置Fig.3 Layout of symmetric nodes in circular tunnel

為保證數(shù)值模擬中模型所受到的滲透壓力條件的統(tǒng)一性，SZ1和SZ2直接通過“INITIAL pp 2.5e5 grad-0.1e4”設(shè)定初始的孔隙水壓力；SZ3則采用FISH函數(shù)及“PROPERTY XXX”命令定義初始的滲透壓力.計算求解完成后，統(tǒng)一選取各圓形隧道模型上z方向上位移變化相同的對稱點(0,0,5)，(0,0,-5),以及x方向上的位移變化點(5,0,0),(-5,0,0)，如圖3所示.

SZ1，SZ2，SZ3求解后，各圓形隧道對稱位置上節(jié)點的位移(s)變化曲線，如圖4所示.圖4中:n為時步.由圖4可知：各模型節(jié)點的位移最終趨于某一個固定值.

SZ3中模型的最大不平衡力記錄圖，如圖5所示.由圖5可知：最大不平衡力隨著求解時步的增加最終接近于零.因此,SZ1,SZ2,SZ3求解計算到8 000個時步時,系統(tǒng)已經(jīng)達到了平衡狀態(tài).

(a) SZ1節(jié)點

(b) SZ2節(jié)點 (c) SZ3節(jié)點圖4 各節(jié)點位移變化Fig.4 Nodal displacement

圖5 SZ3最大不平衡力記錄Fig.5 Maximal unbalanced force of SZ3

由圖4(a)可知：各圓形隧道對稱位置上的節(jié)點位移變化趨勢是大致相同的，即節(jié)點1 301，3 901，1 353，1；在0～500個時步內(nèi)，4個對稱節(jié)點的位移受到外荷載的作用急劇地增加；在500～8 000個時步內(nèi)，節(jié)點3 901的位移在8～9 mm起伏，基本保持不變；節(jié)點1 353的位移在0.5 mm左右保持不變；節(jié)點1的位移在0.6 mm左右保持不變；在500～4 000個時步內(nèi)，節(jié)點1 301的位移緩慢地增加；在4 000～8 000個時步內(nèi)，節(jié)點1 301的位移在4～5 mm之間呈穩(wěn)定狀態(tài).這說明通過控制數(shù)值模擬條件的統(tǒng)一性驗證自定義本構(gòu)模型的方法是成功的.

圖4(c)中節(jié)點3 901的位移最終趨于穩(wěn)定值9 mm左右，而圖4(a)中節(jié)點3 901的位移最終趨于穩(wěn)定值8 mm左右，這說明采用摩爾庫倫模型進行數(shù)值模擬計算求解后得到的節(jié)點位移偏于保守.

圖4(b)中各對稱節(jié)點最終平衡狀態(tài)的位移均比圖6大，以節(jié)點1 353，1為例，圖4(b)中位移分別為0.6,0.6 mm，而圖4(c)中位移分別為0.4,0.4 mm，這說明當(dāng)裂隙巖體中存在初始狀態(tài)下的地下水時，直接通過命令“INITIAL pp”賦予模型初始的孔隙水壓力，其求解計算得到的對稱節(jié)點上的位移比直接通過修正本構(gòu)模型中剛度矩陣得到對稱節(jié)點上的位移更大.

4 結(jié)論

1) 當(dāng)裂隙巖體中存在初始地下水時，考慮初始滲透壓力的影響是非常有必要的.

2) 在充分考慮初始滲透壓力的作用，并將其直接引進本構(gòu)方程中，利用C++二次開發(fā)出來的本構(gòu)模型是成功的.通過修正摩爾庫倫本構(gòu)模型中的剛度矩陣考慮初始滲透壓力的影響會更結(jié)合實際，其數(shù)值模擬結(jié)果相對于直接設(shè)置初始孔隙水壓力時將更加精確.

3) 當(dāng)初始條件統(tǒng)一時，采用自定義本構(gòu)模型數(shù)值模擬出來的結(jié)果與采用摩爾庫倫模型數(shù)值模擬出來的結(jié)果大同小異，模型中各對稱節(jié)點的位移變化趨勢大致相同，但采用摩爾庫倫模型計算得到的結(jié)果相對于采用自定義本構(gòu)模型計算得到的結(jié)果偏于保守.

因此，綜合以上結(jié)論分析，此次二次開發(fā)出來的本構(gòu)模型是成功的.同時，該模型也給一些包含地下水的巖體工程的數(shù)值模擬工作提供了一定的參考價值.

參考文獻：

[1] 朱維申,張強勇.節(jié)理巖體脆彈性斷裂損傷模型及其工程應(yīng)用[J].巖石力學(xué)與工程學(xué)報,1999,18(3):245-249.

[2] 藍航,姚建國,張華興,等.基于FLAC3D的節(jié)理巖體采動損傷本構(gòu)模型的開發(fā)及應(yīng)用[J].巖石力學(xué)與工程學(xué)報,2008,27(3):572-579.

[3] 齊銀萍.裂隙巖體三維損傷流變模型研究及工程應(yīng)用[D].濟南:山東大學(xué),2012：11-40.

[4] 王昆.含裂隙水巷道變形破壞特征研究[D].淮南:安徽理工大學(xué),2014:9-33.

[5] 王昆,趙光明,孟祥瑞，等.含水裂隙巖體本構(gòu)模型及數(shù)值模擬理論研究[J].地下空間與工程學(xué)報,2015,11(4):901-908.

[6] 平楊,鄭少河,白世偉.考慮滲透壓力的裂隙巖體斷裂損傷本構(gòu)模型研究[C]∥中國巖石力學(xué)與工程學(xué)會第六次學(xué)術(shù)會議.北京:中國科學(xué)技術(shù)出版社,2000:134-136.

[7] 鄭少河.裂隙巖體滲流場： 損傷場禍合理論研究及應(yīng)用[D].武漢：中國科學(xué)院武漢巖土力學(xué)研究所,2000:77-118.

[8] 朱維申,李術(shù)才,陳衛(wèi)忠.節(jié)理巖體破壞機理和錨固效應(yīng)及工程應(yīng)用[M].北京：科學(xué)出版社,2002:53-72.

[9] KEMENY J,COOK N G W.Effective moduli, non-linear deformation and strength of a cracked elastic solid[J].J Rock Mech Min Sci and Geomech Abstr,1986,23(2):107-118.

[10] 易順民,朱珍德.裂隙巖體損傷力學(xué)導(dǎo)論[M].北京:科學(xué)出版社,2005:24-36.

[11] 趙延林,曹平,林杭,等.滲透壓作用下壓剪巖石裂紋流變斷裂貫通機制及破壞準(zhǔn)則探討[J].巖土工程學(xué)報,2008,30(4):511-517.

[12] 張電吉.裂隙巖質(zhì)邊坡滲流場與應(yīng)力場耦合分析及工程應(yīng)用[D]. 武漢:中國科學(xué)院武漢巖土力學(xué)研究所, 2003:11-21.

[13] 柴紅保,曹平,趙延林,等.裂隙巖體損傷演化本構(gòu)模型的實現(xiàn)及應(yīng)用[J].巖土工程學(xué)報,2010,32(7):1047-1053.

[14] 曹平,蒲成志.單壓下有序多裂隙脆性材料破壞機制及其簡化模型[J].中國有色金屬學(xué)報,2011,21(10):2659-2668.

[15] 王濤,韓煊,趙先宇,等.FLAC3D數(shù)值模擬方法及工程應(yīng)用： 深入剖析FLAC3D 5.0[M].北京:中國建筑工業(yè)出版社,2015:168-199.

[16] 鄭穎人,孔亮.巖土塑性力學(xué)[M].北京:中國建筑工業(yè)出版社,2010:155-198.

(責(zé)任編輯： 陳志賢 英文審校： 方德平)

Secondary Development on Fractured Rock Constitutive Model Considering Initial Seepage Pressure

ZHANG Weidong， CHEN Shihai

(College of Civil Engineering, Huaqiao University, Xiamen 361021, China)

In order to study of mechanical properties of water-bearing fractured rock, based on Mohr coulomb model, considering the role of initial seepage pressure, flexibility tensor of the constitutive equation is appropriately amended, the constitutive model of jointed rock considering initial seepage pressure is obtained. By comparison, the constitutive model of jointed rock is also obtained without considering initial seepage pressure, according to the same theoretical derivation method and secondary development process. The C++program of two models is compiled to verify the accuracy of the models, through building the same 50 m×50 m×1 m numerical model for circular tunnel in the FLAC3D, three kinds of constitutive model and 4 symmetric point displacements are analyzed in circular tunnel. The result shows that constitutive model introducing initial seepage pressure is successful accurate. Keywords：fractured rock； initial seepage pressure; flexibility tensor; constitutive relation; secondary development

10.11830/ISSN.1000-5013.201703007

2016-08-02

陳士海(1964-)，男，教授，博士，主要從事巖土與地下工程防災(zāi)減災(zāi)，以及工程結(jié)構(gòu)抗震與抗爆的研究.E-mail:cshblast@163.com.

高等學(xué)校博士學(xué)科點專項基金資助項目(20113718110002); 爆炸沖擊防災(zāi)減災(zāi)國家重點實驗室基金資助項目(DPMEIKF201307); 華僑大學(xué)科研基金資助項目(13BS402); 華僑大學(xué)研究生科研創(chuàng)新能力培育計劃資助項目(1400204010)

TU 91

A

1000-5013(2017)03-0319-06