粗粒土應力腐蝕蠕變影響因素的離散元分析

王 瑤，魏進兵，溫仁節，王 琛

(1.四川大學 水利水電學院， 四川 成都 610000；2.四川蜀禹水利水電工程設計有限公司， 四川 成都 610000)

粗粒土是一種工程中廣泛應用的土體材料，大量的國內外試驗結果和監測資料[1-4]均表明其具有明顯的長期變形特性，粗粒土的蠕變性質研究一直是巖土工程領域的研究熱點之一[5]。

土的蠕變性質研究一般包括蠕變機理、蠕變性質與規律、蠕變模型及其應用等方面，其中，蠕變機理研究有助于對蠕變性質與規律的認識，亦有助于建立更符合土體變形規律的蠕變模型。迄今，一般認為粗粒土的蠕變機理包括[6-10]：(1) 土顆粒棱角破碎；(2) 顆粒內部的微裂縫受應力作用隨時間發展為顆粒劈裂；(3) 土顆粒材料本身的流變性質。前兩種蠕變機理都認為粗粒土蠕變的主要細觀機理為顆粒破碎，顆粒位置發生重分布進而導致了蠕變現象。

考慮顆粒破碎的粗粒土蠕變機理研究可基于宏觀試驗和細觀數值試驗展開。宏觀試驗現今主要為土樣級配分析，其方法簡便且可為細觀數值試驗提供細觀參數標定標準；缺點是僅僅能獲取某些特征點的顆粒破碎特性，無法獲取整個蠕變過程中連續的顆粒破碎結果。細觀試驗(如圖像技術)和數值試驗(如離散單元法)可實現連續監測，從而獲取試樣蠕變全過程相關參數。離散單元法[11]等數值試驗方法相對于圖像技術等細觀試驗方法來說，經濟迅捷，并可監測蠕變過程中的顆粒運動、組構和力鏈演變、破碎形式及破碎分布等。Tran等[12]在PFC2D中將堆石顆粒視為由剛性顆粒膠結成的團粒，并讓粒間膠結鍵強度隨時間增長而降低，成功模擬了粗粒土在荷載作用下隨時間增長而產生裂縫直至發生顆粒破碎的現象，并由此來研究堆石壩的蠕變。Potyondy[13]基于亞臨界裂縫擴展理論[14]認為膠結鍵半徑受應力腐蝕而逐漸縮減，在PFC 4.0中編寫了相應的應力腐蝕程序，成功模擬了巖石試樣的蠕變問題。在PFC 5.0版本中，應力腐蝕程序的實現及其模擬粗粒土蠕變的適宜性尚需進一步研究，為進一步模擬粗粒土蠕變細觀機理，也需探討參數對蠕應變的影響并確定適合的細觀模型參數。

本文基于應力腐蝕理論在PFC 5.0中編制了考慮顆粒破碎的粗粒土蠕變離散元程序，分析接觸模型參數中顆粒摩擦系數、膠結鍵粘結強度與顆粒接觸剛度及應力腐蝕模型參數中材料參數、激活應力對蠕變的影響規律及其敏感性；基于室內試驗結果，驗證應力腐蝕模型的應用，標定適用于紅砂巖粗粒土的細觀數值試驗參數，以期為粗粒土蠕變細觀機理的進一步研究奠定基礎。

1 粗粒土應力腐蝕蠕變的數值模型

1.1 應力腐蝕理論與程序

應力腐蝕原理是由巖石斷裂力學中的亞臨界裂縫拓展理論類比而來，將土體視為一種含有微裂縫且不連續的材料，認為應力腐蝕是長期應力作用下導致土體裂縫發展的主要原因[14]。應力腐蝕機制中存在兩個重要的界限值，即發生應力腐蝕與裂縫拓展突變的界限值，分別稱為起裂韌度和斷裂韌度。當應力達到起裂韌度時，微裂縫開始緩慢拓展；當應力接近或超過斷裂韌度時，裂縫發展速度會產生突變，迅速發生顆粒破碎，宏觀上表現為蠕變變形顯著增長；當大量的顆粒破碎后，新的顆粒破碎發展速率趨于緩慢，宏觀上表現為蠕變趨于穩定。

基于這一理論，可在PFC中通過膠結鍵拉應力的函數來模擬膠結鍵的膠結半徑隨時間縮減的過程，從而模擬巖土材料的裂縫隨時間發展的過程[13]。膠結鍵直徑縮減速度公式：

(1)

式中:dc為膠結鍵直徑；σa為激活應力；σt為膠結鍵抗拉強度；σ為膠結鍵上的最大拉應力；β1和β2均為材料參數。作用于膠結鍵上的最大拉應力是判斷是否發生應力腐蝕的基準，當σ≥σa時，才會開始發生應力腐蝕，當σ≥σt時，膠結鍵則會破壞進而發生顆粒破碎。

在進行蠕變的離散元分析時，蠕變時間一般有兩種考慮方式：其一是用黏彈性接觸模型(如：Burgers模型)模擬真實的物理時間。然而因離散元軟件的物理時間步通常小于百萬分之一秒，對幾天乃至幾年的蠕變模擬需借助云計算或超算才能實現；其二是通過考慮蠕變機理(如應力腐蝕模型)，用蠕變機理模型中的時間變量替代真實的物理時間來模擬蠕變，可大幅度節省計算時間，更具有可行性[13]，本文采用該方法進行分析。

在PFC 5.0中編寫的應力腐蝕程序流程如下：(1) 對試樣施加伺服偏應力，判定試樣是否破壞，若試樣破壞，則結束蠕變模擬，若試樣未破壞，則執行應力腐蝕程序；(2) 判定是否達到設置的蠕變模擬時間，若當前蠕變時間已達到設置的模擬時間，則結束模擬，若未達到設置的模擬時長，則執行應力腐蝕函數；(3) 遍歷每一個膠結鍵，判定是否發生應力腐蝕，若不發生應力腐蝕，則結束蠕變模擬，否則計算應力腐蝕速度，縮減膠結鍵半徑，同時更新蠕變時間；(4) 再次施加應力邊界條件，按上述相同步驟反復迭代計算，直至數值模擬滿足設置的蠕變模擬時長、試樣不發生應力腐蝕和試樣破壞三個判定條件之一，即結束蠕變模擬。

1.2 試樣建模

1.2.1 墻體生成

PFC中的墻體一般為剛性墻且無質量，墻體與墻體之間不存在相互作用，故不能運用牛頓第二定律或對墻體直接施加力或彎矩作用，只能通過對墻體賦予速度的方式給顆粒體系施加伺服力作用，墻體與顆粒間的相互作用則通過改變墻體剛度的方式進行調整。數值建模參考室內試驗儀器尺寸，寬100 mm、高200 mm。如圖1所示。試樣包含四面墻體，上下兩面墻分別模擬試樣帽和試樣底座，用于施加伺服偏應力；左右兩面墻為柔性伺服邊界，模擬橡皮膜對試樣施加伺服圍壓。

圖1 試樣墻體示意圖

1.2.2 顆粒生成

顆粒生成包括級配設計和Cluster團粒生成。

在顆粒級配設計中，為更好地分析蠕變過程中顆粒破碎與分維值的關系，同時避免孔隙結構的影響，根據分形理論[15]將土料設計為具有自相似分形結構的級配，有:

(2)

式中:M(d

圖2 級配曲線

在顆粒生成時，每個粒組對應使用一種模板生成顆粒，為使顆粒破碎后體積守恒，模板顆粒設置為兩兩相切。忽略各粒組土顆粒形狀的隨機性，統一用一種模板顆粒代替。其中，20 mm～10 mm模板代表扁長顆粒，10 mm～5 mm模板代表尖銳棱角顆粒，5 mm～2 mm模板代表渾圓度較高顆粒，模板顆粒如圖3所示。

圖3 團粒模板

根據面積等效原理可得每個模板的等效半徑與顆粒數目為：

(3)

式中:R為大顆粒的面積等效半徑，n組成模板的小顆粒個數，ri第i個小顆粒的半徑，通過半徑最終計算等效粒徑與各模板團粒個數見表1。

表1 團粒粒徑和數目

建模時擬使用的Cluster團粒因其能破碎的特性可有效用于模擬顆粒破碎，但在PFC中無法直接生成。為解決這一問題先根據級配利用不可破碎的Clump模板生成試樣，再用自編寫的函數轉換成可破碎的Cluster團粒，如圖4所示。

生成試樣后，由轉換命令給每個接觸兩側球體賦予額外變量，識別不同的接觸，當接觸兩側球體的額外變量相同時，該接觸定義為顆粒內部接觸，否則定義為Cluster之間的接觸。Cluster內部接觸采用平行粘結模型，模擬試樣顆粒的膠結作用；Cluster之間采用線性接觸模型，模擬粗粒土之間的接觸關系。

圖4 數值試樣模型

1.3 基準參數

表2 模型參數值

2 粗粒土應力腐蝕蠕變影響因素分析

在粗粒土蠕變數值試驗分析中，參數敏感性分析是極為重要的一環，通過分析不同參數對蠕變的敏感性，可定性確定對蠕變起主要影響的參數，排除低相關性參數，為今后進一步的蠕變機理研究提供一定的參考[19-21]。本文選取蠕應變作為目標函數值，在基準參數的基礎上，通過改變單一參數，控制其它因素不變來分析變化參數的敏感性。由于各參數的量綱并不統一，且取值相差較大，故引入無量綱參數敏感度函數來度量參數敏感性。針對蠕變問題，參數敏感度函數可定義為蠕應變的相對變化值δε與影響參數相對變化值δx之比為：

(4)

(5)

(6)

式中:Δε為蠕應變的絕對變化值；Δx為參數的絕對變化值；S(x)為參數敏感度函數；ε為蠕變基準值；x為敏感參數基準值。

本文進行圍壓σ3=300 kPa，應力水平S=0.75條件下的參數敏感性分析，分析各參數對最終蠕應變εf的影響，計算敏感度函數值，其值越大參數越敏感。

2.1 接觸模型參數對蠕變的影響

最終蠕應變εf-摩擦系數f關系如圖5(a)所示。隨著摩擦系數的增大，顆粒之間的相對運動愈加困難，所消耗的能量亦越多，蠕應變越小；但這種影響并不是無止境的，因影響粗粒土變形的主要因素還是顆粒破碎，摩擦系數僅能在很小范圍內影響到蠕應變，所以蠕應變減小速率逐漸變緩。顆粒接觸摩擦系數可定性的反映母巖顆粒的表面摩擦特征，顆粒表面愈粗糙的粗粒土蠕變越小。

圖5(b)給出了最終蠕應變εf-膠結鍵抗拉強度σt關系，εf隨σt的增大而減小。原因為膠結鍵抗拉強度的增大使得膠結鍵愈發難以破壞，同時也降低了顆粒應力腐蝕速度，這兩種影響都會導致顆粒破碎量減少，進而使蠕應變降低。表明粗粒土的蠕應變隨母巖抗拉強度的增大而減小。

顆粒剛度kn與顆粒密度ρ對最終蠕應變的影響較小，如圖5(c)和圖5(d)所示。顆粒剛度決定了顆粒的接觸位移大小，故在一定程度上決定了顆粒間的相對位置，但顆粒剛度決定的彈性應變在伺服階段已基本完成，而基于應力腐蝕的蠕變變形主要來自于顆粒裂縫發展導致的顆粒破碎，故顆粒剛度對蠕應變的影響很小；顆粒密度變化僅影響顆粒加速度，而伺服調節始終會使顆粒達到平衡位置，故顆粒密度對粗粒土蠕變的影響不大。

綜上，粗粒土的母巖性質對其蠕變有決定性影響，粗粒土顆粒的表面摩擦系數越大蠕應變越小，抗拉強度越大蠕應變越小。

圖5 接觸腐蝕模型參數-蠕應變關系曲線

2.2 應力腐蝕模型參數對蠕變的影響

應力腐蝕模型參數包括激活應力σa和材料參數β1、β2。蠕應變隨激活應力σa的增大而呈線性減小，如圖6(a)所示。由粗粒土母巖性質決定的激活應力σa作為應力腐蝕的閾值條件，直接影響到參與應力腐蝕的膠結鍵數目，隨著激活應力增大，顆粒破碎急劇減少，蠕應變隨之減小。本文激活應力為0.3 MPa～0.9 MPa，該線性關系在此范圍外或不成立。

蠕應變隨著β1和β2的增大而增大，如圖6(b)和圖6(c)所示。由式(1)可知膠結鍵應力腐蝕速度與β1成正比，β2是應力腐蝕速度函數的指數參數，故當β1、β2增大時，應力腐蝕速度變快，使得相同時間內的顆粒破碎量顯著增大，蠕應變隨之增大。隨著蠕應變的增大，粗粒土孔隙逐漸變小，顆粒位置調整受限，蠕應變增長放緩。

各參數的參數敏感度如表3所示，四個接觸模型參數中f、ρ、kn的敏感度函數值均不大于0.05，僅σt大于0.2，總體表現為對蠕應變影響較不敏感；而應力腐蝕模型參數的敏感度函數值均大于0.2，尤其β2達到了0.3，綜合可知應力腐蝕模型參數比接觸模型參數更加敏感。

2.3 參數標定

在離散元數值試驗中，為了確保計算模型最終能實際模擬土體的宏觀物理力學性質，標定出一套適用于室內試驗結果的細觀參數至關重要，可為進一步分析粗粒土的蠕變機理奠定基礎。

圖6 應力腐蝕模型參數-蠕應變關系曲線

表3 參數敏感度

室內試驗選取粗粒土為人工破碎弱風化紅砂巖，比重由虹吸桶法測得Gs=2.68,制樣干密度通過嘗試法即在承膜桶內將分維值為1.5的試驗土料擊實至人工所能達到的最密實狀態獲得ρd=1.70 g/cm3，孔隙比e=0.57這一干密度較工程中偏小，其原因為試樣剔除了1 mm以下土料，且分維值D=1.5偏小，致試樣具有較大的孔隙結構。試樣尺寸與級配與數模試樣保持一致，采用SR-20型三軸蠕變儀進行三軸飽和固結排水試驗，采用單級加載方式，不同圍壓下的破壞偏應力通過砝碼間隔50 kPa加載，直至軸向應變達到10%或試樣發生持續不穩定變形為止。本文以試驗圍壓σ3=300 kPa，應力水平S=0.25、0.50、0.75和0.90為例進行分析。由圖7室內試驗蠕變結果可知：紅砂巖無黏性粗粒土的蠕應變-時間曲線可分為衰減蠕變和勻速蠕變兩個階段，勻速蠕變階段試樣變形微小，在應力水平較低時迅速趨于穩定。相同圍壓下，蠕變隨應力水平增大而增大，其原因為隨著應力水平的增大，土顆粒接觸應力相應增大后產生了更多的顆粒破碎，進而導致更大的位置調整，致使蠕變隨之增大。這與數值模擬理論分析結果相類似。

圖7 數值模擬與室內試驗蠕變結果對比圖

在接觸模型和應力腐蝕模型參數敏感性分析的基礎上，本文標定了適用于紅砂巖的模型參數標定值，見表2。室內試驗與數值模擬所得的蠕應變εc-時間t關系曲線(σ3=300 kPa)對比如圖7所示，數值曲線的衰減蠕變時長比室內試驗短，后期變形較室內試驗偏小，但大體趨勢一致。結果表明，基于應力腐蝕理論，采用離散元模擬粗粒土的蠕變是可行的，可基于該程序與參數開展粗粒土蠕變細觀機理的進一步研究。

3 結 論

(1) 接觸模型參數中，顆粒摩擦系數和膠結鍵黏結強度對蠕變影響較大，顆粒剛度和密度敏感性較差，粗粒土顆粒母巖表面摩擦系數和抗拉強度越大其蠕應變越小。

(2) 應力腐蝕模型中，激活應力對粗粒土的蠕變有顯著影響，隨著激活應力的增大，能發生應力腐蝕的膠結鍵減少，顆粒破碎量急劇減少，蠕應變減小；材料參數β1和β2均對粗粒土的蠕變有較大影響，蠕應變隨著β1或β2的增大而增大。

(3) 基于應力腐蝕理論進行考慮顆粒破碎的粗粒土蠕變的數值試驗是可行的，并基于粗粒土在不同應力水平下的蠕變室內試驗結果標定了接觸模型和應力腐蝕模型的細觀參數，為粗粒土顆粒破碎蠕變機理的進一步研究奠定了基礎。