橫觀各向同性板巖三軸力學特性及其本構模型

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

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

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

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

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

1 板巖三軸壓縮試驗及力學特性分析

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

1.1 試驗儀器及試樣加載

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

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

圖2 常規三軸壓縮試驗破壞試樣Fig. 2 Failure specimens for conventional triaxial compression tests

1.2 試驗結果

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

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

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

1.3 力學特性分析

基于板巖常規三軸試驗的結果，可以獲得其在不同層理夾角下的三軸力學特性，為后續本構模型的構建及彈性參數的求解提供依據。

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

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

2 彈塑性本構模型

2.1 彈性本構模型

目前，巖土材料的彈性階段主要采用在局部坐標系中的廣義胡克定律來描述。該定律認為在三維空間中，應力與應變呈線性關系，其表達式為[16]：

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

2.1.1 局部坐標彈性本構模型

如圖5所示，當巖體位于局部坐標系中，即不考慮坐標軸旋轉時，板巖的彈性本構模型為：

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

局部坐標中的板巖彈性柔度矩陣[Ce]為：

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

2.1.2 整體坐標彈性本構模型

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

圖6 局部坐標與整體坐標的轉換Fig. 6 Conversion of local coordinates to global coordinates

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

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

式中：Sij(i,j=1~6)為每個單位應力分量所產生的應變分量。

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

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

2.2 塑性本構模型

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

2.2.1 屈服準則

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

式中：F，G，H，L，M，N為模型材料參數，用以描述材料的方向性，其數值通過常規三軸壓縮試驗的結果求解。

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

式中：φ為內摩擦角，c為黏聚力。

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

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

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

2.2.2 流動法則

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

在傳統塑性力學中，常采用關聯流動法則。該法則將塑性勢函數看作屈服函數，即Q(σ)=f(σ)。大量的土工試驗表明[19]：相關聯的流動法則不考慮應力主軸的旋轉，若Mohr-Coulomb 模型采用關聯流動法則，將會導致出現遠大于實際的剪脹變形。綜上所述，關聯流動法則不適用于巖土類材料，根據廣義塑性理論的規定選取非關聯流動法則，Q(σ)≠f(σ)。依據文獻[19]中的方法，將Q(σ)假設為：

式中：β為剪脹角，根據文獻[20]得β=φ/2，R為塑性擬合參數。

2.2.3 硬化定律

硬化定律是描述應力增量對塑性應變影響的準則，它決定了屈服面的變化情況。本文采用塑性主應變為硬化參量H[19]：

2.3 參數求解

2.3.1 彈性參數

基于圍壓為5 MPa時的常規三軸試驗，可得到Ez，Ex，Ey，νyz，νxz，νxy等6 個基本彈性參數，為本構模型的數值實現提供彈性參數依據。

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

式中：dσi為主應力增量，dεi為主應變增量。

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

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

基于5 MPa圍壓作用下水平層理和豎直層理試樣的常規三軸試驗結果，板巖彈性參數的計算結果如表1所示。

表1 板巖彈性參數Table 1 Elastic parameters of slate

2.3.2 塑性參數

板巖彈塑性本構模型的屈服函數中的F和H為模型材料擬合參數，通過分離差應力-應變關系中的塑性應變增量求得其具體數值，為本構模型的數值實現提供塑性參數依據。

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

水平層理試樣的常規三軸試驗中圍壓保持不變，σx=σy，由勢函數得：

徑向塑性應變增量：

軸向塑性應變增量：

其中塑性應變增量可以通過試驗分離。結合常規三軸試驗，聯立式(24)和式(25)得到模型材料擬合參數F=4.67，H=1。

3 本構模型的數值實現

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

3.1 計算模型與參數

建立標準圓柱試樣模型，高100 mm，直徑為50 mm。模型再現常規三軸試驗全過程。在模擬過程中，采用應力加載的方式。模型單元類型為C3D8I，邊界條件采取固定底端豎直方向(Z方向)位移，限制模型頂端和底端的轉動，模型的計算參數如表2所示。

表2 ABAQUS模型計算參數Table 2 Calculation parameters of ABAQUS model

3.2 計算結果及模型驗證

在圍壓為5 MPa的條件下，分別對不同層理夾角的試樣進行了數值模擬計算，其計算結果與板巖常規三軸試驗對比，以驗證板巖本構模型的合理性。本文以不同層理夾角試樣的軸向應變為例，數值模擬與常規三軸試驗的對比結果如圖7所示。

圖7 模型計算結果驗證Fig. 7 Validation of model calculation results

圖7表明：板巖彈塑性本構模型的計算結果與常規三軸試驗的結果基本吻合，差應力-應變的變化趨勢一致，當達到板巖的峰值強度時，應變略微偏小。除90°層理夾角試樣以外，其余層理夾角試樣均在20～30 MPa 范圍內開始出現拐點，這是因為試樣開始出現塑性變形逐漸屈服，與試驗結果一致。數值模擬結果再次體現了板巖的材料方向性對板巖的變形和破壞強度都產生了一定的影響。除此之外，在數值模擬的過程中發現層理夾角為30°，45°和60°的試樣模型計算迭代步數明顯多于水平層理和豎直層理，導致數值計算的速度較慢，其原因是當層理夾角為30°，45°和60°時，UMAT 子程序計算循環次數增加，后期研究可繼續優化。

4 結論

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

2) 建立了板巖彈塑性本構模型，模型彈性部分采用了整體坐標系下的廣義胡克定律，并推導出了整體坐標系中的彈性剛度矩陣；采用改進的Mohr-Coulomb 屈服準則和勢函數、非關聯流動法則和應變硬化準則描述其塑性部分。基于常規三軸試驗結果，求解得到了試樣的彈、塑性參數。

3) 通過對UMAT 子程序的二次開發，完成了板巖彈塑性本構模型在ABAQUS 中的數值實現，數值計算結果與常規三軸試驗結果的對比表明板巖彈塑性本構模型是合理的。