考慮初始滲透壓力的裂隙巖體本構模型二次開發及其驗證

張衛東， 陳士海

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

考慮初始滲透壓力的裂隙巖體本構模型二次開發及其驗證

張衛東， 陳士海

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

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

裂隙巖體; 初始滲透壓力； 柔度張量； 本構關系； 二次開發

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

1 含水裂隙巖體的本構關系

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

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

(2)

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

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

(3)

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

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

τ′=(1-Cv)τ.

(5)

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

[7]中的方法，求得單位體積巖體內裂隙產生的附加彈性應變能φd.

且有

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

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

(8)

1.2 裂隙擴展產生的損傷演化柔度張量

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

大量的現場勘測結果和理論計算表明，巖體中的裂隙在壓剪應力的作用下，將在最大和最小主壓應力組成的平面內沿著最大主壓應力σ1的方向穩定擴展[11]，其翼形分支裂隙發生擴展.

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

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

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

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

(11)

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

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

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

1.2.2 拉剪應力狀態下的損傷演化柔度張量 在拉剪應力狀態下，巖體中的裂隙受拉后,出現張開現象，其裂隙表面的黏結力將變為零,而不能傳遞拉應力.在最大周向應力準則的基礎之上，當分支裂隙尖端瞬時應力強度因子KI0達到巖體的斷裂韌度KlC時，裂隙便開始擴展.

拉剪應力狀態下，分支裂隙尖端瞬時應力強度因子為

(12)

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

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

可得拉剪應力狀態下考慮初始滲透壓力的裂隙巖體其裂隙擴展產生的損傷演化柔度張量[13].

2 裂隙巖體自定義本構的二次開發

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

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

3 模型的驗證及結論

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

模型中，節點施加的初始應力條件為

(14)

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

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

圖2 數值計算模型Fig.2 Numerical calculation model

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

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

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

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

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

(a) SZ1節點

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

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

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

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

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

4 結論

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

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

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

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

參考文獻：

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

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

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

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

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

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

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

[8] 朱維申,李術才,陳衛忠.節理巖體破壞機理和錨固效應及工程應用[M].北京：科學出版社,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] 易順民,朱珍德.裂隙巖體損傷力學導論[M].北京:科學出版社,2005:24-36.

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

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

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

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

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

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

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

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-)，男，教授，博士，主要從事巖土與地下工程防災減災，以及工程結構抗震與抗爆的研究.E-mail:cshblast@163.com.

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

TU 91

A

1000-5013(2017)03-0319-06