楊文劍，劉懷忠*，謝紅強，肖明礫，卓 莉

(1.四川大學(xué) 水利水電學(xué)院，四川 成都 610065；2.四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065)

由堆積碎石層組成的坡體廣泛存在于中國各地，若受工程活動或自然影響致滑成災(zāi)，往往會造成巨大的經(jīng)濟損失和人員傷亡。如：2001年，西藏昌都芒康縣由于巖土體凍脹融縮、風(fēng)蝕干燥引發(fā)山體出現(xiàn)持續(xù)的崩塌碎屑流現(xiàn)象；2010年，貴州關(guān)嶺縣崗烏鎮(zhèn)大寨村發(fā)生特大型崩滑碎屑（石）流災(zāi)害，造成99人死亡失蹤；2013年，西藏墨竹工卡縣澤日山東坡發(fā)生滑坡碎屑流，造成83人死亡失蹤[1]。這種崩滑產(chǎn)生的滑坡碎屑流的宏觀力學(xué)性質(zhì)與降雨導(dǎo)致的泥石夾雜的滑坡泥石流差異較大。由于泥的黏結(jié)作用，泥石流通常存在較大的黏聚力，而碎屑流可視為沒有黏聚力的大量堆積碎石類材料在重力作用下的滾落滑動。研究這類堆積碎石材料的力學(xué)性質(zhì)，可增加對滑坡碎屑流等自然災(zāi)害的認識，從而對這類災(zāi)害進行更有效的預(yù)防治理。

在堆積碎石類材料的力學(xué)特性研究方面，相關(guān)領(lǐng)域?qū)W者現(xiàn)階段主要從現(xiàn)場試驗[2-3]和數(shù)值模擬[4-7]兩個方面進行。由于現(xiàn)場試驗費用昂貴且實施復(fù)雜，大量的試驗研究在時間和資源上的投入成本較大，因此，結(jié)合數(shù)值模擬揭示堆積碎石類材料相關(guān)力學(xué)特性也成為一種主要的研究手段。PFC是一種成熟的商用顆粒離散元計算軟件，廣泛運用于巖土工程、災(zāi)害治理等諸多領(lǐng)域。如：Utili等[8]提出一種基于Morh-Coulomb準則的離散元接觸模型，模擬了邊坡的破壞過程； Wang等[9]利用PFC2D對節(jié)理巖體邊坡進行了穩(wěn)定性分析；Hadjigeorgiou等[10]在PFC3D中引入裂縫系統(tǒng)并研究了硬巖中垂直開挖的邊坡穩(wěn)定性。PFC利用牛頓三定律和應(yīng)力-應(yīng)變法則確定各顆粒間的接觸力及運動方式，顆粒間的接觸特性，即接觸本構(gòu)及其參數(shù)值決定了顆粒材料的宏觀力學(xué)性質(zhì)[11-12]；而接觸本構(gòu)由接觸剛度、摩擦系數(shù)、阻尼等細觀參數(shù)描述，只有正確標定細觀參數(shù)，才能保證離散介質(zhì)數(shù)值模擬結(jié)果的科學(xué)性和可靠性。

由于巖土材料的宏觀、細觀參數(shù)眾多，故細觀參數(shù)的標定工作難度較大，許多學(xué)者基于室內(nèi)壓縮試驗開展巖石介質(zhì)的細觀參數(shù)標定工作，如：Chen等[13]基于巖石三軸壓縮試驗，利用試錯法標定了巖樣的細觀力學(xué)參數(shù)，對含微裂紋的巖石壓縮和拉伸力學(xué)行為進行了研究；Shi等[14]基于巖石單軸壓縮試驗，提出一套巖石細觀參數(shù)的標定流程，并根據(jù)室內(nèi)試驗?zāi)M了花崗巖的力學(xué)行為；Yoon[15]基于單軸壓縮試驗，提出一種利用篩選試驗設(shè)計和中心組合設(shè)計的細觀參數(shù)標定優(yōu)化方法。但是，目前離散元方法中的細觀參數(shù)標定方法研究主要聚焦于巖石材料，對于滑坡碎屑流等顆粒介質(zhì)的細觀參數(shù)標定方法缺乏系統(tǒng)性研究。在目前主流滑坡離散元數(shù)值模擬中，細觀參數(shù)均按經(jīng)驗進行取值[16-18]，忽略了巖土材料宏-細觀參數(shù)聯(lián)系；且不同滑坡體材料具有獨特的物理力學(xué)特性，其細觀參數(shù)必然顯著不同，按經(jīng)驗或類比對細觀參數(shù)取值不利于對滑坡碎屑流災(zāi)害進行正確模擬。因此亟待提出能夠快速、正確地標定堆積碎石類材料細觀參數(shù)的方法。

相較于三軸壓縮試驗，直剪試驗具有設(shè)備輕便、操作簡單、省時節(jié)材的優(yōu)點，對于現(xiàn)場環(huán)境有更好的適應(yīng)性。結(jié)合室內(nèi)堆積碎石直剪試驗，采用PFC3D軟件建立等尺度3維數(shù)值模型，基于線性模型揭示堆積碎石宏-細觀力學(xué)參數(shù)之間的關(guān)系，提出堆積碎石類材料細觀參數(shù)的系統(tǒng)性標定方法；基于此，研究直剪過程中堆積碎石的剪切力學(xué)特性并開展堆積碎石的滑坡模擬試驗，以驗證本文的細觀參數(shù)標定方法在堆積碎石類材料滑坡分析中的適用性及合理性。

1 堆積碎石直剪試驗

本文采用堆積碎石開展3組重復(fù)直剪試驗，試樣如圖1所示，呈平均長寬比為1.4（隨機20粒測量均值）的不規(guī)則粒狀；使用標準篩進行篩分，其粒徑分布集中于兩個粒徑范圍，見表1 ，相關(guān)物理參數(shù)也列于表1，其中，粒間摩擦系數(shù)取該堆積碎石試樣自然休止角的正切值[19]；用于直剪試驗的直剪盒如圖2所示。

圖1 直剪試驗堆積碎石試樣Fig. 1 Accumulated debris in direct shear test

圖2 直剪盒Fig. 2 Direct shear box

表1 堆積碎石物理參數(shù)Tab. 1 Physical parameters of accumulated debris

直剪試驗中，為消除顆粒之間的間隙或無效接觸，首先，采用100 kPa的法向力對堆積碎石散粒進行預(yù)壓；然后，對試樣進行200、300、400 kPa的逐級法向加載，并記錄相應(yīng)荷載下的位移值；之后，在400 kPa的法向荷載作用下進行土樣直剪試驗，剪切速度為0.8 mm/min。

法向荷載施加過程中，試件法向應(yīng)力與法向位移之間的關(guān)系曲線如圖3（a）所示，法向應(yīng)力與法向位移之間近似呈線性正相關(guān)，試樣承受的法向荷載越大，其法向位移也越大。剪切過程中，剪應(yīng)力與剪位移之間的關(guān)系曲線如圖3（b）所示，剪應(yīng)力隨剪切位移的增大而增大，但增長速率逐漸減小。

圖3 室內(nèi)直剪試驗結(jié)果Fig. 3 Results of indoor direct shear test

2 直剪試驗宏細觀參數(shù)關(guān)系

為正確把握細觀參數(shù)對于宏觀性質(zhì)的影響，以合理提出堆積碎石細觀參數(shù)標定方法，建立直剪試驗離散元數(shù)值模型，并對直剪試驗中的宏細觀參數(shù)關(guān)系進行探究。

2.1 直剪離散元模型

基于堆積碎石室內(nèi)直剪試驗開展等尺寸的數(shù)值模擬。剪切盒模型尺寸大小為（直徑×高）61.82 mm×45.75 mm，上下盒等高，模型周圍由剛性墻組成邊界，模擬直剪盒的側(cè)向約束作用。剪切時，下盒固定不動，上盒水平向右錯動進行剪切，如圖4所示。

圖4 3維離散元數(shù)值模型Fig. 4 Three dimensional discrete element numerical model

實際堆積碎石為不規(guī)則狀，數(shù)值模擬采用與實際試樣相同級配的球形顆粒進行等效模擬；進行合理的細觀參數(shù)取值后，堆積碎石組的數(shù)值模擬宏觀力學(xué)性質(zhì)能較好地與實際吻合。根據(jù)室內(nèi)試驗試樣的顆粒級配，在剪切盒內(nèi)部生成與實際級配相匹配的球形顆粒。顆粒間的相對密度為2.79 g/cm3，粒間摩擦系數(shù)為0.735。

對于堆積碎石這類粒間無黏性、有摩擦的顆粒材料，線性模型是最適用的接觸本構(gòu)[20]，故本文的數(shù)值模擬采用線性模型，各細觀參數(shù)見表2。

表2 數(shù)值模擬細觀參數(shù)Tab. 2 Mesoscopic parameters in numerical simulation

直剪數(shù)值模擬過程與室內(nèi)直剪試驗一致，在線性模型中，所涉及的接觸參數(shù)均與時間無關(guān)，剪切速率對直剪試驗宏觀力學(xué)性質(zhì)的影響不大，這一點已由劉方成等[21]等進行了計算論證。因此，為提高計算效率，數(shù)值模擬中在逐級施加法向荷載后，采用2 mm/s的剪切速率進行碎石顆粒的剪切試驗。

2.2 接觸剛度對宏觀彈性常數(shù)的影響

為了從堆積碎石的宏觀力學(xué)響應(yīng)出發(fā)確定粒間接觸剛度kn、ks，基于室內(nèi)直剪試驗結(jié)果，采用Kn和Ks兩個宏觀力學(xué)參數(shù)表征試樣的彈性變形。其中，Kn為法向荷載施加過程中法向應(yīng)力-法向位移曲線的斜率，Ks為剪切過程中剪應(yīng)力-剪位移曲線彈性段（近直線段）的斜率。故可建立如式（1）所示的細-宏觀彈性常數(shù)關(guān)系：

為揭示顆粒接觸剛度與直剪試樣宏觀彈性常數(shù)的內(nèi)在聯(lián)系，基于不同的kn和ks開展多組直剪數(shù)值模擬。圖5和6為不同切向接觸剛度ks條件下，試件宏觀彈性常數(shù)Kn、Ks與kn之間的擬合關(guān)系曲線。

圖5 不同ks下Kn-kn線性趨勢曲線Fig. 5 Linear trend curves of the relationship between Kn and kn under different ks

由圖5和6可知，在切向接觸剛度ks不變情況下，隨著法向接觸剛度kn增大，宏觀彈性常數(shù)Kn、Ks隨之增大，但Kn隨法向接觸剛度kn增大的變化率較大。

圖7和8為不同法向接觸剛度kn條件下，試樣宏觀彈性常數(shù)Kn、Ks與切向接觸剛度ks之間的擬合關(guān)系曲線。由圖7和8可知，法向剛度Kn、切向剛度Ks與切向接觸剛度ks的直線擬合關(guān)系近似一條水平線，切向接觸剛度ks對Kn、Ks的影響相對較小。由此可知，試樣宏觀彈性常數(shù)Kn和Ks主要受粒間法向接觸剛度kn的影響，而受粒間切向接觸剛度ks的影響相對較小。

圖7 不同kn下Kn-ks線性趨勢曲線Fig. 7 Linear trend curves of the relationship between Kn and ks under different kn

圖8 不同kn下Ks-ks線性趨勢曲線Fig. 8 Linear trend curves of the relationship between Ks and ks under different kn

綜上可知：堆積碎石試樣的宏觀彈性常數(shù)Kn及Ks均主要受顆粒間法向接觸剛度kn的影響；切向接觸剛度ks對于其宏觀力學(xué)性質(zhì)的影響遠不及kn。故在一定的切向接觸剛度范圍內(nèi)，可忽略ks變化對堆積碎石試樣宏觀彈性常數(shù)的影響，可將ks作為常量，或者將ks與kn通過一定關(guān)系聯(lián)系起來使式（1）只存在唯一變量kn，即簡化為：

由圖5可知，Kn與kn呈現(xiàn)明顯的線性相關(guān)關(guān)系，可假設(shè)Kn由kn的線性單增函數(shù)表示；由圖6可知，Ks總體呈現(xiàn)隨kn增加而增速減緩、略有波動的增加模式，可假設(shè)Ks由kn的斜率減小的非線性單增函數(shù)表示；又考慮到各常數(shù)的物理意義，當(dāng)kn= 0時，有Kn= 0、Ks= 0，故式（2）可表述為：

圖6 不同ks下Ks-kn線性趨勢曲線Fig. 6 Linear trend curves of the relationship between Ks and kn under different ks

式中，a、b、c均為常數(shù)。

3 細觀參數(shù)的標定方法

基于前文，提出一種堆積碎石細觀參數(shù)標定方法，步驟如下：

1）進行堆積碎石室內(nèi)直剪試驗，得到試驗法向位移-法向應(yīng)力曲線、剪位移-剪應(yīng)力曲線。

2）通過試驗手段測得碎石顆粒試樣級配、孔隙率、密度、摩擦系數(shù)，作為直剪試驗數(shù)值模擬建模計算的基礎(chǔ)參數(shù)，本文參數(shù)如表2所示。

3）建立直剪數(shù)值模型，通過試算確定剛度比k*（kn/ks），通常取1～3[22]；收斂原則為隨著kn的增加，Kn、Ks總體呈增加趨勢且無較大波動，本文取k*=1.25。

4）以kn為變量，開展多次數(shù)值模擬，計算并擬合得出宏觀彈性常數(shù)Kn、Ks分別與kn之間的經(jīng)驗公式。

5）分別以Kn、Ks試驗值反演得到法向接觸剛度kn(n)、kn(s)；若兩者相差較?。?0%），則標定成功，kn取kn(n)、kn(s)的均值即可；若兩者相差較大，則需調(diào)整剛度比k*，重新標定。

標定流程如圖9所示。

圖9 細觀參數(shù)標定流程Fig. 9 Flow of mesoscopic parameters calibration method

利用此標定方法對本文所用堆積碎石進行細觀參數(shù)標定，步驟1）、2）在第2節(jié)已完成，步驟3）、4）、5）的具體過程如下。

圖10、11分別為數(shù)值模擬剛度比k*=1.25時的Kn、Ks與法向接觸剛度kn之間的擬合關(guān)系曲線，擬合公式為：

圖10 k＊=1.25時Kn-kn擬合關(guān)系曲線Fig. 10 Fittng curve of the relationship between Kn and kn when k＊=1.25

圖11 k＊=1.25時Ks-kn擬合關(guān)系曲線Fig. 11 Fittng curve of the relationship between Ks and kn when k＊=1.25

對3次室內(nèi)試驗結(jié)果（圖3）取平均值計算可得，堆積碎石的直剪宏觀彈性常數(shù)分別為Kn=501.39 kPa/mm、Ks=229.27 kPa/mm。根據(jù)式（4）可反演得到粒間法向接觸剛度kn(n)=2.283 1×105N/m、kn(s)=2.313 2×105N/m，兩擬合公式所標定的法向接觸剛度kn相差不大，進一步說明在顆粒材料直剪試驗中宏觀力學(xué)參數(shù)Kn、Ks對粒間切向接觸剛度ks的變化不敏感。本文取kn(n)、kn(s)的均值作為粒間切向接觸剛度的標定值，即kn=2.3×105N/m。

4 室內(nèi)直剪試驗的數(shù)值模擬

以表2的細觀參數(shù)及標定后的接觸剛度kn= 2.3×105N/m、ks= 1.84×105N/m（按k*=1.25取值）為最終參數(shù)，進行堆積碎石直剪試驗數(shù)值模擬計算，監(jiān)測并記錄試樣的位移與應(yīng)力變化，同時對試樣直剪過程中的剪切特性進行研究。

4.1 數(shù)值模擬結(jié)果及分析

圖12為數(shù)值模擬與室內(nèi)試驗的法向應(yīng)力-法向位移曲線對比。由圖12可知，Kn模擬值較試驗值大，偏差為21.4%，其原因是室內(nèi)試驗的碎石顆粒為不規(guī)則顆粒且其棱角的細部結(jié)構(gòu)在壓密過程中存在破碎，其粒間空間更易受外力荷載的擾動，故實際不規(guī)則顆粒粒間空隙較數(shù)值模擬純圓顆粒更易壓密。

圖12 數(shù)值模擬與室內(nèi)試驗法向應(yīng)力-法向位移曲線對比Fig. 12 Comparison of normal stress-normal displacement curves between numerical simulation and laboratory test

圖13為數(shù)值模擬與室內(nèi)試驗的剪應(yīng)力-剪位移曲線對比。由圖13可知：除了彈性階段的剪切過程，Ks模擬值與試驗值基本一致；并且，在剪切過程中的塑性階段，模擬值也與試驗值吻合度較高，證明了本文堆積碎石直剪試驗數(shù)值模擬細觀參數(shù)標定方法的合理性及可行性。

圖13 數(shù)值模擬與室內(nèi)試驗的剪應(yīng)力-剪位移曲線對比Fig. 13 Comparison of shear stress-shear displacement curves between numerical simulation and laboratory test

4.2 顆粒接觸力鏈網(wǎng)絡(luò)研究

顆粒體系內(nèi)部的接觸互相連接逐漸形成能夠傳遞荷載的通道，即力鏈，這些通道錯綜復(fù)雜地交織，形成力鏈網(wǎng)絡(luò)[23]。本文針對直剪試驗數(shù)值模擬中的力鏈網(wǎng)絡(luò)演變，從微觀層面研究直剪過程中碎石骨架受力情況。

圖14為在400 kPa法向荷載作用下，剪切過程中試樣內(nèi)部的力鏈網(wǎng)絡(luò)演變過程。由圖14（a）可知，試樣剪切前的力鏈網(wǎng)絡(luò)強力鏈較少，且強力鏈方向受法向荷載影響呈豎向分布的特征。在剪切初期，隨著上部剪切盒向右運動帶動試樣剪切變形，力鏈網(wǎng)絡(luò)發(fā)生重構(gòu)，強力鏈數(shù)量開始增加且力鏈方向開始傾斜，如圖14（b）所示。由圖14（c）可見，隨著剪位移的增加，強力鏈顯著沿對角線方向聚集，試樣左下方及右上方幾乎沒有強力鏈分布。研究表明，隨著剪切位移的增加，力鏈網(wǎng)絡(luò)中強力鏈數(shù)量隨之增加，強力鏈分布由豎直方向逐漸向?qū)蔷€方向演變，而剪應(yīng)力也隨著強力鏈數(shù)量增加和方向偏轉(zhuǎn)而增大。

圖14 試樣內(nèi)部力鏈網(wǎng)絡(luò)演變Fig. 14 Evolution of force chain network in specimen

4.3 剪脹性與顆粒速度場研究

圖15為剪切過程中試樣的體積變化率-剪位移曲線及不同階段試樣顆粒速度場。由圖15可知：試樣剪脹性與顆粒速度場之間存在著密切聯(lián)系；在曲線點A處，法向加載施加完成，但試樣還未開始剪切，體積變化率為0，顆粒速度場雜亂分布；在曲線點B處，試樣正處于剪縮階段，從顆粒速度場也可見顆粒從剪切盒左上方向右下方擠壓，上部顆粒整體從左往右運動的趨勢，故顆粒間彼此貼緊，體積有收縮趨勢；隨著剪切的進行，位移增大至曲線點C時，試樣處于回脹階段，從顆粒速度場可見，除左下角外，顆粒體系的整體運動趨勢為斜向右上運動，此時由于顆粒間咬合摩擦，造成顆?？缭蕉a(chǎn)生相對錯動，故體積有膨脹趨勢。

圖15 試樣體積變化率-剪位移曲線及特征時刻試樣速度場Fig. 15 Curve of volume change rate-shear displacement and velocity field of specimen at specific time

5 堆積碎石滑坡模擬試驗

為進一步驗證本文細觀參數(shù)標定在堆積碎石滑坡分析時的適用性，進行堆積碎石滑坡模擬試驗并進行數(shù)值模擬對比。試驗滑槽模型如圖16所示，室內(nèi)試驗采用的碎石材料與圖1中材料一致，重200 g；滑坡數(shù)值模擬的顆粒為圓球形，其級配、摩擦系數(shù)、密度、接觸剛度等各參數(shù)均與上文直剪試驗數(shù)值模擬所確定的參數(shù)一致。

圖16 滑坡模擬試驗滑槽模型Fig. 16 Slide model of landslide simulation test

試驗時，以擋板將碎石顆粒組擋在滑槽中離地高度20 cm處，隨即釋放擋板，使碎石組沿滑槽自由下滑至堆積槽內(nèi)堆積成型。數(shù)值模擬過程與室內(nèi)試驗保持完全一致。

碎石滑坡堆積形態(tài)對比如圖17所示。由圖17（a）可以看出：室內(nèi)試驗及數(shù)值模擬結(jié)果的俯視堆積形態(tài)均近似呈半橢圓形，橢圓中心均在遠離滑板一側(cè)的立板中軸線上，短半軸均為堆積槽寬度（10 cm）；但是，數(shù)值模擬結(jié)果的長軸為24.0 cm，較實際的22.5 cm偏寬，相差6.67%。由圖17（b）可以看出：室內(nèi)試驗及數(shù)值模擬結(jié)果的正視堆積形態(tài)均近似呈正態(tài)分布，峰值點均位于滑槽中心線上，數(shù)值模擬與室內(nèi)試驗最大堆積高度分別為2.3、 2.7 cm，相差14.8%；通過輪廓線對比可以看出，數(shù)值模擬堆積高度整體較實際偏低。由圖17（c）可以看出：室內(nèi)試驗及數(shù)值模擬的側(cè)視堆積形態(tài)均近似呈半正態(tài)分布，遠離滑板一側(cè)為最高點；通過輪廓線對比可以看出，左側(cè)部分數(shù)值模擬與實際擬合良好，右側(cè)部分數(shù)值模擬結(jié)果較實際偏低。

圖17 碎石堆積形態(tài)對比Fig. 17 Comparison of gravel accumulated morphology

分析認為，數(shù)值模擬所得碎石堆積形態(tài)與室內(nèi)試驗基本一致，但數(shù)值模擬的堆積范圍較實際偏大（6.67%），堆積高度較實際偏低（14.8%）。這些誤差是由數(shù)值模擬采用的球型顆粒較實際不規(guī)則顆粒更易滾動所導(dǎo)致的。數(shù)值模擬的滑坡體外輪廓線與室內(nèi)試驗整體吻合較好，進一步驗證了本文細觀參數(shù)標定在模擬堆積碎石滑坡時的適用性。

6 結(jié)論與展望

基于線性接觸模型開展堆積碎石直剪試驗，揭示了宏觀力學(xué)性質(zhì)與粒間細觀接觸參數(shù)之間的相關(guān)關(guān)系，提出顆粒材料細觀參數(shù)標定方法并進行驗證，主要結(jié)論如下：

1）在碎石直剪數(shù)模試驗中，碎石的宏觀彈性參數(shù)主要受粒間法向接觸剛度kn的影響，受粒間切向接觸剛度ks的影響相對較小。

2）提出一種基于室內(nèi)直剪試驗的堆積碎石細觀參數(shù)標定方法，該方法利用兩個彈性常數(shù)對細觀參數(shù)進行標定，操作簡便，耗材耗時均較小，可重復(fù)性高；并通過直剪試驗數(shù)值模擬與室內(nèi)試驗對比驗證了該方法的適用性和準確性。

3）通過直剪數(shù)模試驗分析了碎石剪切過程中顆粒粒間力鏈網(wǎng)絡(luò)演變過程，隨著剪位移和剪應(yīng)力的增加，力鏈網(wǎng)絡(luò)中的強力鏈數(shù)量增加，且強力鏈分布由豎直方向逐漸向?qū)蔷€方向偏轉(zhuǎn)。

4）土樣在剪切過程中，顆粒速度場由最初的雜亂分布，逐漸呈現(xiàn)斜向下近似水平的運動趨勢，顆粒間彼此貼緊，體積收縮；隨著剪切位移的增大，顆粒間咬合摩擦造成顆?？缭蕉a(chǎn)生相對錯動，速度場呈現(xiàn)斜向上運動趨勢，發(fā)生剪脹現(xiàn)象。

5）通過本文直剪試驗細觀參數(shù)標定方法所得到的堆積碎石細觀參數(shù)，可較好地應(yīng)用于滑坡堆積碎石的數(shù)值模擬，其模擬堆積形態(tài)、堆積范圍均與實際相符。