橫觀各向同性板巖三軸力學(xué)特性及其本構(gòu)模型

陳曄磊，李地元，王立川 ，李化云，張鵬飛，吳劍

(1. 西華大學(xué) 建筑與土木工程學(xué)院，四川 成都 610039；2. 中南大學(xué) 資源與安全工程學(xué)院，湖南 長沙 410083；3. 中南大學(xué) 土木工程學(xué)院，湖南 長沙 410075；4. 中鐵十八局集團(tuán)有限公司，天津 300222；5. 中鐵西南科學(xué)研究院有限公司，成都 611731)

巖石的各向異性特征是巖石力學(xué)領(lǐng)域的研究重點(diǎn)之一[1]，橫觀各向同性是各向異性的一種特殊形式，其結(jié)構(gòu)特征、力學(xué)參數(shù)和應(yīng)力-應(yīng)變關(guān)系與垂直于各向同性平面的旋轉(zhuǎn)軸對稱[2]。板巖具有顯著橫觀各向同性和結(jié)構(gòu)面弱化特征，其破壞模式與層理夾角密切相關(guān)[3]。中國西部區(qū)域多山，地質(zhì)條件復(fù)雜，廣泛分布著大量板巖。隨著交通基礎(chǔ)設(shè)施的大量建設(shè)，沿線將不可避免地出現(xiàn)大量穿越板巖地層的隧道等地下工程。因板巖層理發(fā)育，層間膠結(jié)能力弱，導(dǎo)致工程結(jié)構(gòu)常出現(xiàn)隆起、彎折、斷裂等嚴(yán)重病害現(xiàn)象[4]。為此，國內(nèi)外不少學(xué)者對橫觀各向同性巖體的力學(xué)特性[5-6]及其本構(gòu)模型建立[7]等方面開展了大量研究。巖土體的本構(gòu)模型對有限元計(jì)算有著直接的影響[8-9]，彈塑性本構(gòu)模型包含了屈服準(zhǔn)則、流動(dòng)法則、硬化定律等理論。Mohr-Coulomb 準(zhǔn)則和Drucker-Prager準(zhǔn)則在巖土材料的塑性分析中得到了廣泛的應(yīng)用[10-12]，彌補(bǔ)了經(jīng)典塑性力學(xué)僅適用于金屬等材料的不足，但這2 種屈服準(zhǔn)則主要適用于各向同性材料的分析。在對橫觀各向同性巖體屈服準(zhǔn)則的改進(jìn)中，層理夾角和軟弱層的力學(xué)特性是需要被考慮的因素。層理夾角對本構(gòu)模型的影響不容忽視，與實(shí)際地下工程的情況相符。以隧道工程為例，如圖1 所示，層理與洞壁切線的夾角隨巖石位置的變化而變化(A位置夾角約為45°，B位置夾角約為0°)。在一些本構(gòu)模型的研究中[12-13]，通常會(huì)從微觀角度去考慮軟弱層的影響，這從一定程度上提高了計(jì)算精度，但計(jì)算繁瑣。當(dāng)軟弱層很小或計(jì)算大型開挖模擬時(shí)，這種方法是不經(jīng)濟(jì)的[14]。橫觀各向同性材料本構(gòu)模型參數(shù)較多，需多組試驗(yàn)才能確定模型參數(shù)[15]，所以參數(shù)的確定是值得研究的。

圖1 層理夾角與隧道輪廓線的關(guān)系Fig. 1 Relationship between bedding angle and tunnel profile

綜上所述，本文以板巖為研究對象，通過常規(guī)三軸壓縮試驗(yàn)的結(jié)果，分析其力學(xué)特性與層理夾角間的關(guān)系。在實(shí)驗(yàn)基礎(chǔ)上，推導(dǎo)出在整體坐標(biāo)中的彈性剛度矩陣；采用完全等效模型將層理面效應(yīng)納入巖體的連續(xù)體描述中，建立適用于板巖的彈塑性本構(gòu)模型，并求解模型參數(shù)。最后通過UMAT 子程序二次開發(fā)將其應(yīng)用于ABAQUS 數(shù)值模擬軟件中，結(jié)合常規(guī)三軸試驗(yàn)的結(jié)果驗(yàn)證其合理性。

1 板巖三軸壓縮試驗(yàn)及力學(xué)特性分析

巖石三軸壓縮試驗(yàn)主要包括常規(guī)三軸壓縮試驗(yàn)(σ1>σ2=σ3)和真三軸壓縮試驗(yàn)(σ1>σ2>σ3)，本文開展的是圓柱體試樣的常規(guī)三軸壓縮試驗(yàn)。試驗(yàn)所采用的圍壓分別為5 MPa 和10 MPa(根據(jù)板巖取樣地所在隧道工程的地應(yīng)力量級(jí)大小確定)，試樣層理夾角分別為0°，30°，45°，60°和90°，以此分析圍壓和層理夾角對板巖力學(xué)特性的影響，并為后續(xù)本構(gòu)模型的建立及參數(shù)的求解提供可靠的試驗(yàn)數(shù)據(jù)。

1.1 試驗(yàn)儀器及試樣加載

本文采用GCTS RTX-1500 巖石三軸儀，對不同層理夾角的板巖試樣(自然風(fēng)干狀態(tài))開展常規(guī)三軸壓縮力學(xué)試驗(yàn)。該試驗(yàn)機(jī)軸向最大試驗(yàn)力為1 500 kN，最大圍壓可達(dá)140 MPa。試樣取自峨漢高速豹貍崗隧道中的板巖，加工制成不同層理傾角的標(biāo)準(zhǔn)圓柱試樣。試樣直徑為(50±3) mm，高度為100 mm，兩端不平整度不超過0.05 mm，端面垂直于試樣軸線，最大偏差不超過0.25°。

針對不同層理夾角的板巖試樣進(jìn)行圍壓分別為5 MPa 和10 MPa 的常規(guī)三軸壓縮試驗(yàn)。本次試驗(yàn)圍壓以0.05 MPa/s的速度加載至預(yù)設(shè)值，并在試驗(yàn)過程中保持不變。軸壓以0.01 MPa/s的速度預(yù)加載至2 MPa，待試樣上方柱狀鋼座接觸壓力機(jī)上壓頭后，以0.05 MPa/s的速度加載至試樣破壞。試樣破壞后如圖2所示。

圖2 常規(guī)三軸壓縮試驗(yàn)破壞試樣Fig. 2 Failure specimens for conventional triaxial compression tests

1.2 試驗(yàn)結(jié)果

圖3為不同層理夾角試樣在不同圍壓下的差應(yīng)力-應(yīng)變關(guān)系。由圖3 可知，圍壓為5 MPa 和10 MPa時(shí)，在試樣的差應(yīng)力-應(yīng)變關(guān)系圖中均可觀察到明顯的峰值跌落現(xiàn)象。但在水平層理試樣中，當(dāng)圍壓為5 MPa時(shí)，由于加載過程中試樣壓碎破壞的響聲明顯，表現(xiàn)出較強(qiáng)的脆性特征，為了保護(hù)試驗(yàn)儀器，防止試樣碎屑掉入壓力倉，在未觀察到明顯的應(yīng)力跌落現(xiàn)象時(shí)便停止了加載，但仍可以看出試樣達(dá)到應(yīng)力峰值后下降的趨勢。本文主要研究巖石試樣破壞前的彈塑性變形階段，所以并不影響最終分析結(jié)果。

如圖3 所示，板巖的差應(yīng)力-應(yīng)變關(guān)系受巖石層理夾角和圍壓的影響。巖石試樣的破壞強(qiáng)度隨著圍壓的增加而增大，初始階段沒有觀察到差應(yīng)力-應(yīng)變曲線上凹，說明無明顯微裂隙閉合階段，主要發(fā)生彈性變形。在達(dá)到破壞強(qiáng)度之前，差應(yīng)力-應(yīng)變曲線下凹，此時(shí)巖石試樣發(fā)生應(yīng)變硬化，彈性模量減小，產(chǎn)生部分塑性變形。達(dá)到峰值強(qiáng)度之后，巖石試樣的微裂隙擴(kuò)張、貫通，呈現(xiàn)應(yīng)變軟化現(xiàn)象。最后，巖石試樣進(jìn)入殘余強(qiáng)度階段，可認(rèn)為巖石試樣已失去承載能力。

圖3 不同層理夾角試樣在2 種圍壓下的差應(yīng)力-應(yīng)變關(guān)系Fig. 3 Differential stress-strain curves of specimens with different bedding angles under two confining pressures

1.3 力學(xué)特性分析

基于板巖常規(guī)三軸試驗(yàn)的結(jié)果，可以獲得其在不同層理夾角下的三軸力學(xué)特性，為后續(xù)本構(gòu)模型的構(gòu)建及彈性參數(shù)的求解提供依據(jù)。

如圖4所示，試樣峰值強(qiáng)度與層理夾角呈非線性關(guān)系。當(dāng)圍壓為5 MPa 時(shí)，層理夾角為45°的巖石試樣的峰值強(qiáng)度最小，其值為112.4 MPa，層理夾角為0°的巖石試樣的峰值強(qiáng)度最大，其值為164.6 MPa。當(dāng)圍壓為10 MPa 時(shí)與圍壓為5 MPa 時(shí)的規(guī)律一致，峰值強(qiáng)度最小為45°層理夾角試樣，峰值強(qiáng)度最大為0°層理夾角試樣。圍壓從5 MPa到10 MPa 對0°～90°層理夾角的試樣峰值強(qiáng)度的提升幅度分別為59%，67%，70%，67%和69%。試樣殘余強(qiáng)度與圍壓呈正相關(guān)，隨著圍壓的增加，巖石試樣的殘余強(qiáng)度也隨之增大，巖石的塑性增強(qiáng)。試樣殘余強(qiáng)度與層理夾角呈非線性關(guān)系，60°層理夾角試樣的殘余強(qiáng)度最小。圍壓對巖石試樣徑向應(yīng)變的影響小于對軸向應(yīng)變的影響以及對應(yīng)力的影響。如圖3(d)所示，當(dāng)圍壓由5 MPa 增加到10 MPa 后，巖石試樣彈性階段的軸向應(yīng)變增大了17%，徑向應(yīng)變增大了6%。這是因?yàn)閲鷫旱脑黾樱箮r石內(nèi)部的晶體凝聚力增大，巖石的破壞表現(xiàn)為滑移破壞，圍壓又從一定程度上限制了其橫向滑移，從而增大了巖石的峰值強(qiáng)度，增強(qiáng)了巖石的塑性。

圖4 2種圍壓下板巖試樣峰值強(qiáng)度與層理夾角的關(guān)系Fig. 4 Relationship between the peak strength and bedding angle of slate specimens under two confining pressures

2 彈塑性本構(gòu)模型

2.1 彈性本構(gòu)模型

目前，巖土材料的彈性階段主要采用在局部坐標(biāo)系中的廣義胡克定律來描述。該定律認(rèn)為在三維空間中，應(yīng)力與應(yīng)變呈線性關(guān)系，其表達(dá)式為[16]：

式中：E為彈性模量，ν為泊松比，G=E/2(1+ν)，εi(i=x,y, z)為主應(yīng)變分量，σi(i=x,y, z)為主應(yīng)力分量。

2.1.1 局部坐標(biāo)彈性本構(gòu)模型

如圖5所示，當(dāng)巖體位于局部坐標(biāo)系中，即不考慮坐標(biāo)軸旋轉(zhuǎn)時(shí)，板巖的彈性本構(gòu)模型為：

圖5 局部坐標(biāo)系中的板巖試樣Fig. 5 Slate samples in local coordinate system

局部坐標(biāo)中的板巖彈性柔度矩陣[Ce]為：

式中：Ex=Ey，為xy層面內(nèi)的彈性模量；νxy為xy層面內(nèi)的泊松比；Ez為垂直于層面的彈性模量；νyz=νxz為垂直于層面的泊松比；Gyz=Gzx為垂直于層面的剪切模量。

2.1.2 整體坐標(biāo)彈性本構(gòu)模型

在實(shí)際工程中，最小主應(yīng)力與巖體層理面的夾角會(huì)因開挖后的應(yīng)力重分布而發(fā)生改變。如圖6所示，假定層理面始終平行于局部坐標(biāo)系中的xy平面，最小主應(yīng)力始終平行于整體坐標(biāo)系中的x＇y＇平面。當(dāng)局部坐標(biāo)系繞x軸旋轉(zhuǎn)時(shí)，最小主應(yīng)力與巖體層理面的夾角θ發(fā)生改變，該夾角簡稱為層理夾角。

圖6 局部坐標(biāo)與整體坐標(biāo)的轉(zhuǎn)換Fig. 6 Conversion of local coordinates to global coordinates

韓昌瑞等[17]基于正交各向異性材料的彈性理論，通過坐標(biāo)軸旋轉(zhuǎn)，實(shí)現(xiàn)了廣義胡克定律從局部坐標(biāo)系到整體坐標(biāo)系的轉(zhuǎn)換，從而體現(xiàn)了層理夾角因開挖后應(yīng)力重分布而發(fā)生的改變。將廣義胡克定律在整體坐標(biāo)系中表示為：

整體坐標(biāo)中的板巖彈性柔度矩陣[Ce]'為：

式中：Sij(i,j=1~6)為每個(gè)單位應(yīng)力分量所產(chǎn)生的應(yīng)變分量。

當(dāng)層理夾角θ=0 或90°時(shí)，[Ce]'可以退化為局部坐標(biāo)系中的彈性柔度矩陣[Ce]。通過對整體坐標(biāo)系中的彈性柔度矩陣[Ce]'求逆，可得整體坐標(biāo)系中的彈性剛度矩陣[De]'：

式中：Kij(i,j=1~6)為每個(gè)單位應(yīng)變分量所對應(yīng)的應(yīng)力分量。

2.2 塑性本構(gòu)模型

巖土材料應(yīng)力-應(yīng)變關(guān)系的塑性部分通常需要屈服準(zhǔn)則、流動(dòng)法則、硬化定律等理論來描述。但目前常用的屈服準(zhǔn)則只適用于各向同性材料，關(guān)聯(lián)流動(dòng)法則一般適用于金屬類材料，不適用于巖土材料。本文考慮層理夾角的影響，對Mohr-Coulomb準(zhǔn)則進(jìn)行了推廣，并采用非關(guān)聯(lián)流動(dòng)法則描述塑性流動(dòng)，從而推導(dǎo)出合理的板巖塑性本構(gòu)模型。

2.2.1 屈服準(zhǔn)則

以Mohr-Coulomb 定律為基礎(chǔ)的屈服準(zhǔn)則被廣泛采用，但是該準(zhǔn)則的屈服面存在尖角，計(jì)算時(shí)會(huì)出現(xiàn)奇點(diǎn)問題，且常用于各向同性材料中，不能反映材料方向性。HILL[16]基于米塞斯準(zhǔn)則提出了能反映材料方向性，適用于巖土材料的屈服準(zhǔn)則：

式中：F，G，H，L，M，N為模型材料參數(shù)，用以描述材料的方向性，其數(shù)值通過常規(guī)三軸壓縮試驗(yàn)的結(jié)果求解。

對于Mohr-Coulomb 這類在π平面上為非圓形的屈服函數(shù)，可用廣義Von·mises 圓曲線去逼近，從而解決了計(jì)算時(shí)的奇點(diǎn)問題，表達(dá)式如下所示[18]：

式中：φ為內(nèi)摩擦角，c為黏聚力。

在式(8)中的應(yīng)力張量第二不變量J2引入模型材料參數(shù)，使其能反映材料方向性。并且因?yàn)榘鍘r的各向異性對稱于z軸，模型材料參數(shù)F=G=M=L，H=N，則J2可被推廣為J＇2：

將式(13)代入式(8)可得：

式中：θσ為應(yīng)力洛德角μσ為洛德參數(shù)，I1為應(yīng)力張量第一不變量，I1=σx+σy+σz。

2.2.2 流動(dòng)法則

流動(dòng)法則是描述塑性流動(dòng)的規(guī)律。Von·mises假設(shè)存在某種塑性勢函數(shù)Q(σ)，其塑性流動(dòng)方向與塑性勢函數(shù)Q(σ)的梯度或外法線方向保持一致。基于塑性位勢理論可得：

在傳統(tǒng)塑性力學(xué)中，常采用關(guān)聯(lián)流動(dòng)法則。該法則將塑性勢函數(shù)看作屈服函數(shù)，即Q(σ)=f(σ)。大量的土工試驗(yàn)表明[19]：相關(guān)聯(lián)的流動(dòng)法則不考慮應(yīng)力主軸的旋轉(zhuǎn)，若Mohr-Coulomb 模型采用關(guān)聯(lián)流動(dòng)法則，將會(huì)導(dǎo)致出現(xiàn)遠(yuǎn)大于實(shí)際的剪脹變形。綜上所述，關(guān)聯(lián)流動(dòng)法則不適用于巖土類材料，根據(jù)廣義塑性理論的規(guī)定選取非關(guān)聯(lián)流動(dòng)法則，Q(σ)≠f(σ)。依據(jù)文獻(xiàn)[19]中的方法，將Q(σ)假設(shè)為：

式中：β為剪脹角，根據(jù)文獻(xiàn)[20]得β=φ/2，R為塑性擬合參數(shù)。

2.2.3 硬化定律

硬化定律是描述應(yīng)力增量對塑性應(yīng)變影響的準(zhǔn)則，它決定了屈服面的變化情況。本文采用塑性主應(yīng)變?yōu)橛不瘏⒘縃[19]：

2.3 參數(shù)求解

2.3.1 彈性參數(shù)

基于圍壓為5 MPa時(shí)的常規(guī)三軸試驗(yàn)，可得到Ez，Ex，Ey，νyz，νxz，νxy等6 個(gè)基本彈性參數(shù)，為本構(gòu)模型的數(shù)值實(shí)現(xiàn)提供彈性參數(shù)依據(jù)。

WANG 等[15]通過水平層理試樣求解Ez，νyz和νxz：

式中：dσi為主應(yīng)力增量，dεi為主應(yīng)變增量。

對于豎直層理試樣，圍壓不變，dσx=dσy=0，可得如下各彈性參數(shù)：

式中：dεr為徑向應(yīng)變，dεr=(dεx+dεy)/2。

基于5 MPa圍壓作用下水平層理和豎直層理試樣的常規(guī)三軸試驗(yàn)結(jié)果，板巖彈性參數(shù)的計(jì)算結(jié)果如表1所示。

表1 板巖彈性參數(shù)Table 1 Elastic parameters of slate

2.3.2 塑性參數(shù)

板巖彈塑性本構(gòu)模型的屈服函數(shù)中的F和H為模型材料擬合參數(shù)，通過分離差應(yīng)力-應(yīng)變關(guān)系中的塑性應(yīng)變增量求得其具體數(shù)值，為本構(gòu)模型的數(shù)值實(shí)現(xiàn)提供塑性參數(shù)依據(jù)。

式(15)中的塑性因子dλ可以表示為：

水平層理試樣的常規(guī)三軸試驗(yàn)中圍壓保持不變，σx=σy，由勢函數(shù)得：

徑向塑性應(yīng)變增量：

軸向塑性應(yīng)變增量：

其中塑性應(yīng)變增量可以通過試驗(yàn)分離。結(jié)合常規(guī)三軸試驗(yàn)，聯(lián)立式(24)和式(25)得到模型材料擬合參數(shù)F=4.67，H=1。

3 本構(gòu)模型的數(shù)值實(shí)現(xiàn)

本節(jié)基于ABAQUS 6.14，對UMAT 子程序進(jìn)行二次開發(fā)。用戶材料子程序UMAT 是ABAQUS用于定義特殊材料屬性的二次開發(fā)接口[21]。

3.1 計(jì)算模型與參數(shù)

建立標(biāo)準(zhǔn)圓柱試樣模型，高100 mm，直徑為50 mm。模型再現(xiàn)常規(guī)三軸試驗(yàn)全過程。在模擬過程中，采用應(yīng)力加載的方式。模型單元類型為C3D8I，邊界條件采取固定底端豎直方向(Z方向)位移，限制模型頂端和底端的轉(zhuǎn)動(dòng)，模型的計(jì)算參數(shù)如表2所示。

表2 ABAQUS模型計(jì)算參數(shù)Table 2 Calculation parameters of ABAQUS model

3.2 計(jì)算結(jié)果及模型驗(yàn)證

在圍壓為5 MPa的條件下，分別對不同層理夾角的試樣進(jìn)行了數(shù)值模擬計(jì)算，其計(jì)算結(jié)果與板巖常規(guī)三軸試驗(yàn)對比，以驗(yàn)證板巖本構(gòu)模型的合理性。本文以不同層理夾角試樣的軸向應(yīng)變?yōu)槔瑪?shù)值模擬與常規(guī)三軸試驗(yàn)的對比結(jié)果如圖7所示。

圖7 模型計(jì)算結(jié)果驗(yàn)證Fig. 7 Validation of model calculation results

圖7表明：板巖彈塑性本構(gòu)模型的計(jì)算結(jié)果與常規(guī)三軸試驗(yàn)的結(jié)果基本吻合，差應(yīng)力-應(yīng)變的變化趨勢一致，當(dāng)達(dá)到板巖的峰值強(qiáng)度時(shí)，應(yīng)變略微偏小。除90°層理夾角試樣以外，其余層理夾角試樣均在20～30 MPa 范圍內(nèi)開始出現(xiàn)拐點(diǎn)，這是因?yàn)樵嚇娱_始出現(xiàn)塑性變形逐漸屈服，與試驗(yàn)結(jié)果一致。數(shù)值模擬結(jié)果再次體現(xiàn)了板巖的材料方向性對板巖的變形和破壞強(qiáng)度都產(chǎn)生了一定的影響。除此之外，在數(shù)值模擬的過程中發(fā)現(xiàn)層理夾角為30°，45°和60°的試樣模型計(jì)算迭代步數(shù)明顯多于水平層理和豎直層理，導(dǎo)致數(shù)值計(jì)算的速度較慢，其原因是當(dāng)層理夾角為30°，45°和60°時(shí)，UMAT 子程序計(jì)算循環(huán)次數(shù)增加，后期研究可繼續(xù)優(yōu)化。

4 結(jié)論

1) 開展了不同層理夾角、不同圍壓條件下的板巖常規(guī)三軸壓縮試驗(yàn)，基于試驗(yàn)結(jié)果分析了板巖的力學(xué)特性。板巖試樣的峰值強(qiáng)度、殘余強(qiáng)度都隨著圍壓的增加而增大，但試驗(yàn)范圍內(nèi)的圍壓對其徑向變形的影響有限，板巖峰值強(qiáng)度和殘余強(qiáng)度均與層理夾角呈非線性關(guān)系，45°層理夾角試樣強(qiáng)度最小，0°層理夾角試樣強(qiáng)度最大。

2) 建立了板巖彈塑性本構(gòu)模型，模型彈性部分采用了整體坐標(biāo)系下的廣義胡克定律，并推導(dǎo)出了整體坐標(biāo)系中的彈性剛度矩陣；采用改進(jìn)的Mohr-Coulomb 屈服準(zhǔn)則和勢函數(shù)、非關(guān)聯(lián)流動(dòng)法則和應(yīng)變硬化準(zhǔn)則描述其塑性部分。基于常規(guī)三軸試驗(yàn)結(jié)果，求解得到了試樣的彈、塑性參數(shù)。

3) 通過對UMAT 子程序的二次開發(fā)，完成了板巖彈塑性本構(gòu)模型在ABAQUS 中的數(shù)值實(shí)現(xiàn)，數(shù)值計(jì)算結(jié)果與常規(guī)三軸試驗(yàn)結(jié)果的對比表明板巖彈塑性本構(gòu)模型是合理的。