豫南燕山期花崗巖蠕變特性及非線性蠕變損傷模型

唐志強， 吉鋒*， 許漢華， 馮文凱， 何蕭

(1.成都理工大學地質災害防治與地質環境保護國家重點實驗室， 成都 610059； 2.中國有色金屬工業昆明勘察設計研究院有限公司， 昆明 650051； 3.云南省巖土工程與地質災害重點實驗室， 昆明 650051)

秦嶺-大別造山帶作為世界上規模最大的超高壓變質帶，其后造山階段的殼幔作用以及造山帶的折返、隆升和剝蝕是個復雜的動力學過程；具體的演化主要分為俯沖、同碰撞和后碰撞三個階段，后碰撞階段主要出現走滑變形和伸展斷裂活動，出現大規模的花崗巖漿活動[1]。信陽雞公山一帶花崗巖屬于燕山期花崗巖體的一部分，主要出露在桐商斷裂和龜梅斷裂以南，北大別構造帶上，以高鉀鈣堿性I型花崗巖為主。伴隨著豫南地區眾多抽水蓄能電站的建設施工，該地區深埋地下洞室的圍巖長期穩定性問題亟待解決。許多研究表明，深埋地下洞室的圍巖變形失穩破壞往往與時間密切相關[2-3]。因此，研究該地區特有高鉀鈣堿性 I型花崗巖的蠕變力學性質，對于豫南燕山期花崗巖地區深埋地下洞室圍巖長期穩定性評價及支護設計具有極其重要的意義。

目前針對蠕變模型的研究方法主要是通過對試樣的蠕變試驗結果進行分析，根據試驗蠕變變形表現出的特征，對模型中可能存在的基本元件進行辨識，并將各基本元件進行串、并聯等組合，預測在不同應力作用史及變化的物理環境下巖體的變形及破壞與時間的關系[4]。元件模型通常是線性模型，而許多工程實踐和試驗研究表明，巖石的蠕變往往是一個非線性的過程。鑒于元件模型存在缺陷，國內外學者針對如何更吻合的描述巖石蠕變全過程展開了大量研究。Zhao等[5]通過連接胡克體、胡克體與塑性滑塊并聯、Kelvin體和廣義Bingham體，建立了包含多種變形分量的非線性彈塑性流變模型；Liu等[6]采用分數階黏壺替代西原正夫模型中的Newton黏壺，從而得到巖石的非線性蠕變本構模型；韓陽等[7]提出了一種非線性黏壺元件替換Burgers模型中的兩個線性黏壺元件，建立了一種非定常參數Burgers模型；李任杰等[8]用分數階微積分軟體元件替換Kelvin模型中整數階黏壺元件,同時引入非定常化非線性黏性元件,采用元件組合方式,構建了巖石硬性結構面蠕變模型。這些研究主要采用非線性元件與傳統元件串并聯或替代傳統元件來建立非線性模型。但非線性模型組合似乎只是為了是理論曲線與試驗曲線相吻合，非線性元件的參數卻缺乏具體的物理意義[9-10]。

巖石的蠕變破壞本質是巖石在長期荷載作用下，內部裂紋不斷產生與發展，巖石材料力學性質逐漸劣化損傷直至破壞的過程[11]。在蠕變模型中考慮材料力學參數隨應力狀態以及時間的損傷變化，建立巖石的蠕變損傷模型更能反映蠕變的真實狀態。李杰林等[12]通過定義基于核磁共振孔隙度的巖石凍融損傷因子與基于體應變的蠕變損傷因子，并將損傷因子引入基本流變元件中，建立了凍融巖石的非線性蠕變損傷本構模型及方程；Lu等[13]通過對缺陷砂巖進行凍融循環試驗和三軸壓縮試驗，建立了凍融砂巖在長期荷載作用下的損傷演化方程；鐘祖良等[14]基于試驗數據將傳統Burgers模型改進為全損傷Burgers蠕變模型，該模型能較好地描述花崗質片麻巖的蠕變特性。綜上，中外學者針對特定地質環境下的不同巖石構建蠕變損傷本構模型，并應用于計算模擬，顯著提升了巖石變形預測精度。然而，目前針對豫南地區燕山期高鉀鈣堿性 I型花崗巖的蠕變力學性質的研究仍較薄弱。隨著豫南某水電站大型地下廠房洞室開挖，對該地區燕山期花崗巖的蠕變力學特性進行研究顯得極為重要。

為此，以該電站地下廠房內的燕山晚期中粒二長花崗巖為研究對象，采用YSJ-01-00型巖石三軸蠕變試驗機開展了三軸壓縮蠕變實驗，分析了花崗巖在不同圍壓下長期強度的變化規律。構建了能夠描述加速蠕變階段的蠕變損傷模型，并對改進后的西原模型進行了參數辨識及擬合。研究成果可為豫南燕山期花崗巖地區深埋地下洞室圍巖長期穩定性評價提供試驗依據和理論指導。

1 花崗巖三軸蠕變試驗

1.1 研究區花崗巖礦物成分及分布規律

豫南某水電站工程區位于秦嶺～大別造山帶上，所在區域劃分為3個一級大地構造單元，包括中朝準地臺、秦嶺地槽褶皺系和揚子準地臺。進一步可劃分為9個二級大地構造單元，即中朝準地臺的華熊臺緣坳陷、華北坳陷，秦嶺褶皺系的北秦嶺褶皺帶、南秦嶺褶皺帶、桐柏—大別褶皺帶、南陽—襄樊坳陷、潢川坳陷和揚子準地臺的上揚子臺隆、兩湖斷坳，如圖1所示。

圖1 區域大地構造單元劃分

該電站的上水庫和地下廠房區均位于出露在龜山-梅山斷裂(產狀為N85°W/SW60°～80°)和桐柏-商城斷裂帶(產狀為N60°W/SW50°～80°)之間，北大別構造帶上的燕山期花崗巖株內，為高鉀鈣堿性I型花崗巖。為揭示此種花崗巖的礦物成分及分布規律，沿燕山期花崗巖株橫剖面方向取6個樣品(圖2)進行磨片測試。

圖2 研究區巖性分布及取樣位置

根據巖石薄片鑒定結果，研究區花崗巖的礦物成分主要為堿性長石、斜長石、石英，次要礦物見黑云母、白云母，副礦物為鋯石和磁鐵礦。其中堿性長石含量為33%～39%，呈半自形-他形短柱狀；斜長石含量為32%～40%，多為自形-半自形的厚板狀；石英含量為21%～25%，無色，他形粒狀；黑云母含量為2%～3%，呈褐黃色-褐色多色性，鱗片狀，具一組極完全解理；白云母含量為1%～3%，無色，具閃突起，干涉色鮮艷。各類巖石礦物顆粒粒徑沿花崗巖株橫剖面方向呈先減小后增大的分布規律(圖3)。

圖3 顆粒粒徑分布規律

1.2 蠕變試驗設備

本次試驗采用的試驗設備是由成都理工大學地質災害防治與地質環境保護國家重點試驗室自主研發的YSJ-01-00型巖石三軸蠕變試驗機(圖4)，該試驗機的圍壓范圍為0～30 MPa，軸向荷載范圍為0～600 kN，軸向荷載測試精度和圍壓測試精度均為0.5%F.S，軸向荷載和圍壓恒定時間可以保持6個月以上，滿足本次力學試驗的要求。

圖4 試驗加載系統

1.3 試樣制備

本次試驗的巖石試樣取自河南某水電站大型地下廠房內，巖性為花崗巖；此種花崗巖呈灰白色，中～細粒結構，塊狀構造，礦物粒徑在1～4 mm。參照《水利水電工程巖石試驗規程》的試件加工精度要求，制作成Φ50 mm×100 mm的標準圓柱樣。試樣的高度、直徑允許誤差為±0.3 mm；試件兩端面的不平行度允許偏差為±0.05 mm；端面應垂直于試件軸線，允許偏差為±0.25°。

該水電站地下廠房區巖體應力主要為自重應力場，最大主應力σ1為8.59 MPa，巖石單軸抗壓強度平均值為101.17 MPa。為模擬研究區地應力條件下花崗巖的力學特征及蠕變特性，分別在圍壓σ3=5 MPa和σ3=7.5 MPa條件下進行及常規三軸壓縮試驗，物理力學特性指標如表1所示。

表1 花崗巖物理力學參數

1.4 蠕變試驗方案

基于單軸抗壓試驗和常規三軸壓縮試驗結果，分別設置圍壓σ3=5 MPa和σ3=7.5 MPa兩個試驗組，加載方式采用單試件逐級加載試驗方式“陳氏加載法”[15-16]，根據單軸抗壓強度以及不同圍壓下三軸最大偏應力進行應力水平分級，具體加載應力路徑如圖5所示。采用荷載控制加載速率方式，加載速率為2.5 MPa/min，考慮到蠕變試驗周期較長，每一級荷載的持續時間為24 h，然后重復操作直至試樣破壞，停止加載。整個實驗過程中需保持圍壓恒定不變。

圖5 不同圍壓下加載應力路徑

2 試驗結果分析

2.1 蠕變特征分析

根據花崗巖三軸蠕變力學試驗曲線，可分析不同應力狀態下花崗巖的蠕變變形隨時間變化的規律。圖6為花崗巖試樣在圍壓σ3=5、7.5 MPa時的三軸壓縮蠕變軸向應變-時間曲線。表2為花崗巖試樣在不同圍壓下的蠕變試驗數據。

圖6 試樣軸向應變-時間

分析表2中數據可知：同一圍壓條件下，每一級應力增量相等時，隨著軸向偏應力的增加，瞬時軸向應變增量逐漸減小，表明試樣的彈性模量逐漸增大，試樣呈現應變硬化特性；兩種圍壓條件下，軸向蠕變總量與瞬時彈性應變的比值εc/ε0均不超過10%(其中，ε0為瞬時彈性應變，εc為軸向蠕變總量)，表明花崗巖試樣在蠕變全過程中軸向總應變以瞬時彈性應變為主；圍壓7.5 MPa，偏應力42.5 MPa時，軸向蠕變總量εc為0%，表明此時偏應力水平并未達到花崗巖試樣的蠕變起始應力，驗證了堅硬巖也存在蠕變變形閾值現象[17]；當圍壓升高時，在相當的偏應力水平作用下，軸向蠕變總量會減小，如圍壓5 MPa，偏應力150 MPa時，軸向蠕變總量為0.084%，圍壓7.5 MPa，偏應力170 MPa時，軸向蠕變總量為0.04%；后者的偏應力水平比前者高20 MPa，但軸向蠕變總量卻小于前者，表明圍壓提高，可以制約花崗巖的蠕變變形。

表2 不同圍壓作用下花崗巖試樣三軸壓縮蠕變實驗數據

2.2 長期強度分析

巖石蠕變長期強度對評價工程長期穩定性而言是極其重要的。通過三軸壓縮蠕變試驗可知，附加應力水平較小時，巖石的蠕變變形速率會逐漸減小，最終趨于穩定，蠕變變形也會趨于收斂，巖石不會發生破壞。當附加應力水平大于或等于某一定值時，試樣的蠕變變形會隨時間逐漸增大，且巖石的蠕變速率會在某一階段迅速增加，試樣最終發生破壞。這個應力水平的臨界值就認為是巖石的長期強度。

常用的巖石蠕變長期強度的確定方法有3種：一是等時應力-應變曲線法；二是過渡蠕變法(主要根據廣義Kelvin模型和Burgers模型對數據的辨識度進行判斷)；三是蠕變曲線第一拐點法。

采用等時應力-應變曲線法來確定長期強度。試驗時每一級荷載作用持續時間為24 h，分別采取0、2、6、12、18、24 h對應的應變量，繪制六簇等時應力-應變關系曲線，如圖7所示。

圖7 不同圍壓下等時應力-應變圖

觀察六簇等時應力-應變曲線(圖7)可知，曲線主要由兩部分組成，一部分由左側斜率較大線段組成，另一部分由斜率較小線段組成，中間存在一個轉折點，其對應的軸向應變突然增大，故該轉折點所對應的軸向應力即為所求的長期強度值。

當圍壓為5 MPa時，巖石長期強度為180.23 MPa；當圍壓為7.5 MPa時，長期強度為168.9 MPa；與常規三軸試驗相比，在5 MPa和7.5 MPa條件下，長期強度分別下降了6.37%和21.0%。

2.3 破壞特征分析

花崗巖試件在圍壓5 MPa、7.5 MPa時，三軸壓縮蠕變試驗的破壞形態如圖8所示。從破壞時間來看，由圖6(b)的最后一級蠕變破壞曲線可知，圍壓5 MPa下試樣蠕變總時長約為1.69 h，發生加速蠕變至試樣破壞的總歷時為0.11 h，其加速蠕變持續時間較為短暫，僅占總時長的6.5%。花崗巖試樣的蠕變破壞呈現出脆性破壞特征。

圖8 不同圍壓下試樣蠕變破壞形態及素描圖

從破壞形式來看，兩種不同圍壓條件下，試樣均表現為雙剪切面“Y”形剪切破壞，且都發育有一條主貫通裂紋，主裂紋附近存在破壞較徹底的情況，產生大量粉末狀固體顆粒，分析其原因是花崗巖在破壞過程中，由于附加應力克服了主裂紋面上的法向應力所產生的摩擦力以及巖石本身的內聚力，導致試樣沿主裂紋處發生了強烈的剪切滑移，破壞了剪切面上巖石的晶體結構，所以產生了大量粉末狀的固體顆粒。圍壓7.5 MPa的試樣表現出整體濕潤，這是由于試樣受力破壞時將包裹試樣的熱縮管炸裂，圍壓缸內的液壓油浸入所致。

由圖8可知，圍壓5 MPa下，試樣破壞時2條主裂紋的宏觀破壞面夾角約36°，主裂紋交匯處位于試樣底部；圍壓7.5 MPa下，試樣破壞時2條主裂紋的宏觀破壞面夾角約26°，主裂紋交匯處位于試樣中部。對比可知，圍壓7.5 MPa時，主裂紋沿軸線方向的擴展長度及沿徑向的擴展角度均小于圍壓為5 MPa時，但破壞時的軸向偏應力卻比前者大22.5 MPa。因此，相比于軸向偏應力，圍壓對剪切裂紋的擴展具有更為明顯的控制作用。

綜上，花崗巖試樣在三軸壓縮蠕變試驗下的蠕變破壞形式為脆性剪切破壞。

3 非線性損傷蠕變模型研究

由圖6分析可知，當偏應力施加過程中，巖石試樣會產生瞬時彈性變形，則蠕變模型中應存在彈性元件；當附加應力恒定后，巖石試樣蠕變變形速率會逐漸減小趨于穩定，蠕變變形量會隨時間逐漸減小直至收斂，則蠕變模型中應存在與時間相關的黏性元件；當附加偏應力高于長期強度值且蠕變應變達到某一值時，巖石試樣會進入加速蠕變階段，產生不可恢復的塑性變形并導致最終斷裂破壞，則蠕變模型中應存在變形不可恢復的塑性元件。傳統西原模型由胡克體元件、黏彈性體及黏塑性體(Newton體)通過串聯組合而成，能夠較為全面的反應巖石的黏彈塑性變形性質，但對于加速蠕變階段的描述還不夠準確。分析其原因在于傳統西原模型中的Newton的黏性系數為定值，然而在實際蠕變試驗過程中，巖石會發生損傷，巖石內部本身存在的一些隨機分布的微裂紋、孔洞等缺陷在外界影響因素作用下會逐步擴大導致巖石性能劣化，巖石的黏性系數就會發生改變，尤其是加速蠕變階段的黏性系數與蠕變損傷密切相關[18-19]。

因此，考慮蠕變過程中巖石發生的損傷累積，來描述黏性系數隨時間的變化規律，得到一個考慮損傷的變黏性系數黏壺元件，如圖9所示。

圖9 考慮損傷的變黏性系數黏壺元件

圖9中的變黏性系數黏壺元件實質上就是基于損傷累積的Newton體，根據Newton體的本構關系可得出該變黏性系數黏壺元件的本構關系為

(1)

考慮荷載長期作用對黏性系數的影響，可得

η1(D)=η1(1-DC)

(2)

式(2)中：η1為變黏性系數黏壺元件的初始黏性系數；DC為損傷變量，0≤DC<1。

根據文獻[20-22]可知，巖石蠕變過程中的損傷變量與時間呈負指數函數關系，即

DC=1-e-αt

(3)

式(3)中：α為與花崗巖材料性質相關的系數；t為時間。

將式(3)代入式(2)中，可得到變黏性系數黏壺元件的黏性系數與時間關系為

η1(D)=η1e-αt

(4)

再將式(4)代入式(1)中，并對其積分，可得變黏性系數黏壺元件的一維本構方程為

(5)

在三維應力作用下，變黏性系數黏壺元件的三維蠕變本構方程為

(6)

式(6)中：ε為變黏性系數黏壺元件的應變值；σS為長期強度值。

通過將上述得到的變黏性系數黏壺元件替換傳統西原模型中的Newton體，對于三軸壓縮蠕變試驗過程中發現的蠕變變形閾值現象，考慮引入一個開關元件來表征，該開關元件以蠕變起始應力σcr為界限；當σ<σcr時，模型表現出彈性元件性質，只能發生彈性變形，無蠕變變形；當σ≥σcr時，該開關元件失去作用，模型發生彈性變形后，在長期荷載的作用下可發生蠕變變形。得到了基于損傷累計的改進西原模型，用于描述巖石蠕變的全過程尤其是加速蠕變階段的特性。改進的西原模型如圖10所示。

E0、E1及η0分別為瞬時彈性模量、黏彈性模量和黏滯系數

基于傳統西原模型的本構方程，并結合式(6)，可以推導得出改進西原模型的三維蠕變本構方程為

(7)

式(7)中：K0、G0及G1分別為體積模量、剪切模量和黏彈性剪切模量。

4 蠕變損傷模型驗證及參數辨識

為了驗證基于巖石損傷累積得到的改進西原模型的適用性和合理性。根據花崗巖三軸壓縮蠕變試驗結果，對模型進行參數辨識，并得到相應的流變力學參數。借助于七維高科有限公司獨立開發的軟件1-stOpt，采用Levenberg-Marquardt優化算法(即最小二乘算法)+通用全局優化算法對模型參數進行逐次線性迭代求解，并根據相關性系數反映模型參數與蠕變試驗數據的擬合精度。

在常規三軸蠕變試驗中，滿足σ2=σ3，根據胡克定律有

(8)

限于篇幅，僅對圍壓5 MPa的三軸壓縮蠕變試驗數據進行參數辨識，得到各應力水平下的模型參數如表3所示。

表3 圍壓5 MPa時，改進西原模型參數辨識結果

將計算出的模型參數帶入改進的西原模型蠕變方程，得到改進西原模型理論計算曲線和蠕變試驗曲線的對比，如圖11所示。

圖11 圍壓5 MPa時蠕變試驗曲線與改進西原模型理論曲線對比

根據圖11可知，改進的西原模型能夠對花崗巖的蠕變全過程特征進行較好的描述，準確反映了最后一級加載過程中出現的蠕變三階段特性，尤其是對加速蠕變階段的擬合效果也較好，擬合相關系數為0.973 2。因此，改進的西原模型能夠較好地描述河南某水電站大型地下廠房洞室花崗巖的蠕變力學特征，可為豫南燕山期花崗巖地區地下洞室圍巖長期穩定性評價及支護設計提供試驗依據和理論指導。

5 結論

以豫南信陽地區燕山期花崗巖為研究對象，在查明巖石分布情況及礦物成分基礎上，通過設計不同圍壓下的花崗巖三軸壓縮蠕變試驗，對蠕變速率、長期強度及蠕變破壞特征三方面進行綜合分析，構建了考慮損傷累計的非線性損傷變模型，得出以下結論。

(1)工程區燕山期花崗巖的礦物成分主要為堿性長石、斜長石、石英，次要礦物見黑云母、白云母；各類巖石礦物顆粒粒徑沿花崗巖株橫剖面方向呈先減小后增大的分布規律。

(2)花崗巖在蠕變全過程中軸向總應變以瞬時彈性應變為主，且存在蠕變變形閾值；隨著偏應力水平的增加，巖石試樣的彈性模量不斷增大，呈現出應變硬化特性。

(3)花崗巖在三軸壓縮蠕變試驗下的破壞形式為脆性剪切破壞，且圍壓比軸向偏應力對剪切裂紋的擴展具有更為明顯的控制作用；相比于常規三軸壓縮試驗所得到的峰值強度，花崗巖的長期強度分別下降了6.37%和21.0%。

(4)基于試驗結果發現的蠕變變形閾值現象及蠕變過程中巖石發生的損傷累積，構建了改進的西原模型。根據試驗數據對改進西原模型進行了驗證及參數辨識，試驗值與理論值吻合較好，驗證了該非線性蠕變損傷本構模型的適用性及合理性。