基于離散元法研究顆粒球度對粗粒土抗剪強度的影響

李遠征， 魏玉峰， 雷 壯， 何 亮

(成都理工大學 地質災害防治與地質環境保護國家重點實驗室， 四川 成都 610059)

1 研究背景

粗粒土是一種松散的顆粒集合體，由不同粒徑，不同性質的顆粒互相充填而成，作為一種散體材料，粗粒土粒間作用力、幾何排列方式與組成等對其宏觀力學性質有重要影響[1-2]。土的結構性被譽為21世紀土力學的核心問題。眾多學者從宏、細、微觀等不同尺度在粗粒土結構、強度、整體穩定性等方面開展了深入研究[3]，顆粒形態、排列、破碎等因素與宏觀強度特性的聯系也已經被研究者所關注。常規的室內試驗只能從宏觀層面反映粗粒土應力與應變的關系，而無法從顆粒層面探究其細觀機理。且由于試驗設備的限制，試驗過程中顆粒的破裂、旋轉、位移和重排列等都難以直接觀測。自20世紀70年代提出了離散元法后,隨著該理論的發展，越來越多的國內外學者運用離散元法模擬室內三軸試驗和直剪試驗，研究粗粒土組構與其宏觀力學之間的關系[4-7]，顆粒形狀是影響最直接的因素之一。

為了研究顆粒形狀對粗粒土各種宏觀力學特性的影響，首先要探討如何對顆粒形狀進行量化。有關顆粒形狀量化方面的研究方法較多，劉廣等[8]定義球度作為顆粒形狀的量化指標，證明了球度指標的簡便和可行性。劉清秉等[9]提出可以從3個層次對不規則顆粒形狀進行描述：球形度、磨圓度、粗糙度。石崇等[10]和Mattrand等[11]提出應用數學方法對顆粒形狀進行描述，如分形理論、傅里葉分析等。同時一些學者通過數值模擬和室內試驗探究顆粒形狀對不同材料力學特性的影響。史旦達等[12-13]通過模擬砂土直剪試驗研究了顆粒形狀對砂土強度和變形的影響以及與剪切帶厚度的關系。楊貴等[14]采用水泥凈漿澆筑方法制備試樣，通過三軸剪切試驗研究了不同圍壓條件下顆粒形狀對材料強度和變形特性的影響，揭示了隨著顆粒形狀的變化其內摩擦角、極限應力比等的變化規律。以往的研究很大程度上表明了顆粒形狀對材料力學性質有很大影響，但由于顆粒形狀的復雜性，很多學者對顆粒形狀的量化指標缺乏科學的定義。

鑒于離散元法在粗粒土細觀力學機理研究上的巨大優勢，本文基于PFC3D程序的簇單元(clump),由基本圓球型顆粒按不同排列方式組成其余兩種不規則顆粒并生成試樣，對3種不同形狀的顆粒單元進行量化，模擬室內直剪試驗，探究顆粒形狀的變化對粗粒土抗剪強度特性的影響及其細觀機理，同時揭示黏聚力c和內摩擦角φ的變化規律。最后根據室內常規直剪試驗得出的物理結果進行雙向驗證。

2 數值直剪試驗

2.1 顆粒形狀指標量化

顆粒形狀的量化分析，目前主要為采用分形理論和傅里葉分析等數學方法進行描述以及對顆粒形狀參數進行定義。以往研究中，或者涉及的數學計算異常繁雜，或者對顆粒形狀量化指標缺乏科學的定義。為此作者在前人研究的基礎上，從二維層面描述顆粒形狀的指標圓度出發，在三維層面上定義顆粒的球度，最大程度來描述顆粒整體外形與圓球的相似度以及顆粒表面的凹凸程度。

二維層面顆粒圓度F的定義如下[15]：

F=Fs/Fp

(1)

式中：Fp為與實際顆粒等周長的圓面積;Fs為實際顆粒面積;F為顆粒的圓度。

基于二維層面對圓度的定義，首先對顆粒從空間六種方向(前后左右上下)進行拍攝，運用Matlab軟件對拍攝的二維圖像進行處理后提取圓度尺寸信息，對提取到的數據進行加權平均，得出在三維層面定義顆粒形狀的指標球度公式。

(2)

式中：S為顆粒的球度;F1～F6分別為顆粒6種方向(前后左右上下)的投影圓度；α1、α2、β1、β2、γ1、γ2分別為顆粒6種方向投影圖中兩正交軸長軸L與短軸I的比值。

2.2 簇單元的建立和模型的生成

卵狀石英砂巖顆粒表面光滑，顆粒局部并無棱角度變化，本文基于選取并處理后的3種真實顆粒(分別為:Ⅰ型顆粒為表面光滑無棱角近似圓形；Ⅱ型顆粒為表面光滑無棱角呈長條型；Ⅲ型顆粒表面光滑無棱角呈橢圓型)，利用公式(2)計算出3種顆粒的球度后，再通過PFC3D中的“clump”命令，讓基本圓球型顆粒通過一定的排列方式黏結成其余2種球度的顆粒(分別為類長條型和類橢圓型)，生成的顆粒表面粗糙度均相同，如圖1所示。通過clump命令生成的顆粒都可視為剛性體。

球度數值由Matlab軟件處理圖像信息并通過公式(2)得到，真實顆粒與數值模型球度值見表1。因篇幅有限，具體求解過程在此不再贅述。

本文生成顆粒試樣的順序，先生成圓球型試樣，然后按照“體積相等”、“質量中心等效”、“任意旋轉”原則[16]將基本圓球型顆粒單元替換成類長條型和類橢圓型顆粒單元，最后生成其余兩種試樣。

圖1 真實顆粒與數值模型顆粒

數值模型圓球型顆粒類橢圓型顆粒類長條型顆粒真實顆粒圓球型顆粒類橢圓型顆粒類長條型顆粒1.00.7260.4460.9670.7430.435

(1)體積相等。即在替換過程中保持類長條型和類橢圓型顆粒的體積與圓球型顆粒相等。需要計算出縮放系數ζ，通過ζ對組成兩種不規則顆粒的各個從屬球體的大小進行縮放。

(3)

式中：V1為不規則顆粒的體積;V2為被替換的圓球體積。

(2)質量中心等效。不同形狀顆粒替換時保持兩種不規則顆粒的質量中心與基本圓球型顆粒的質量中心一致。不規則顆粒單元的質量中心由組成其從屬顆粒的質量中心經計算得到。

(4)

式中：mi為不規則顆粒中第i個附屬顆粒的質量；xi為不規則顆粒中第i個附屬顆粒的質心坐標；xc為質量中心坐標。

(3)任意旋轉。在替換過程中對顆粒進行任意旋轉，使生成的顆粒避免朝著同一個方向，保證試樣在各個方向的力學性質一致。

假設某一顆粒單元在三維中的坐標為(x，y，z),同時繞x、y、z坐標軸進行旋轉，旋轉角度分別為(θx,θy,θz)，旋轉過后的坐標為(x′，y′，z′)，其轉換關系可用矩陣表示為：

(5)

式中:R1，R2，R3可表示為：

本試驗中模擬的材料為無黏性粗粒土，黏聚力應當按零考慮，因此選用線性接觸剛度模型。

數值模擬生成的3種試樣尺寸均為500 mm×500 mm×500 mm，其中圓球型顆粒的最大粒徑為20 mm，最小粒徑為15 mm，按指定級配生成試樣，試樣如圖2所示。

圖2 3種不同類型顆粒試樣

2.3 細觀參數確定

目前通過室內試驗方法很難直接獲得PFC3D所使用的細觀參數。通常情況下，當室內試驗結果與數值模擬結果相近時，可以認為得到的細觀參數是可行的，數值模擬結果也是正確的。本文參考已有文獻和前人研究成果，經過反復計算、測試與調整，選定了一組較為理想的PFC3D細觀參數，見表2。

表2 數值模擬顆粒力學參數

2.4 數值直剪試驗

數值試驗直剪儀尺寸為500 mm×500 mm×500 mm，由10面剛性墻和2面擋墻組成(圖3)，1,2墻分別構成頂板墻和底板墻，3,4,5,6墻為上剪切盒的側壁墻，7,8,9,10墻為下剪切盒的側壁墻，11,12墻為防止顆粒溢出的擋墻。在剪切過程中，設定豎直方向的壓力分別為100 kPa、200 kPa、300 kPa、400 kPa，經過循環后達到了給定的法向應力狀態。

試樣在穩定的垂直壓力下進行剪切，下剪切盒不動，上盒沿水平方向對試樣進行剪切，剪切速率設為0.1 mm/s，剪切位移設置200 mm。為獲得剪切過程中的全面數據，在試驗過程中設定歷史(history)進行了監測。剪切過程如圖4所示。

圖3 數值模型 (S=1.0)

圖4 直剪試驗的剪切過程(S=0.726)

3 顆粒形狀對抗剪強度的影響

3.1 顆粒形狀對粗粒土直剪試驗應力應變的影響

如圖5為數值模擬試驗得到的3種球度的粗粒土試樣在4種法向荷載作用下剪切位移-剪應力曲線圖。

由圖5可得出，不同球度粗粒土數值試樣的剪切位移-剪應力曲線規律相似，均是先硬化后逐漸軟化，呈單峰曲線型，且球度越小的粗粒土，越過峰值應力后的軟化現象愈明顯，隨著法向荷載的增大，這一現象更為突出。

將數值模擬得到的剪切位移-剪應力曲線進行數據整理，提取3種球度的粗粒土在4種法向荷載下的峰值應力、殘余強度以及峰值位移數值，如表3所示。

表3 數值模擬峰值應力-位移與殘余強度表

由表3數據，得到隨著球度的變化不同法向荷載下，粗粒土峰值應力、殘余強度以及峰值位移三者的變化規律，如圖6所示。

圖6可以看出，在相同法向荷載下，隨著形狀指標球度S的減小，粗粒土峰值應力、殘余強度和峰值位移均顯著增加，但增長趨勢趨于平緩，如在400 kPa法向荷載下，球度S由1.0減小到0.726時，峰值應力、殘余強度和峰值位移分別增加84.0、41.4 kPa和55.8 mm，但由0.726減小到0.446時,則三者分別增加65.0、15.0 kPa和25.2 mm。同時，隨著法向荷載減小，增長趨勢也趨于平緩，如在400 kPa法向荷載下，當球度S由1.0減小到0.446時，峰值應力、殘余強度和峰值位移共分別增加149.0、56.1 kPa和81.0 mm；但在100 kPa法向荷載下，則共分別增加45.0、17.5 kPa和26.3 mm。

另外，當球度相同時，隨著法向荷載的增加，峰值應力、殘余強度和峰值位移均逐漸增加，增加趨勢越來越明顯且呈現出有規律的遞增現象，如圖7所示。當球度S為0.446，法向荷載由100 kPa增加到400 kPa時，峰值應力分別增加48.3、59.6、65.3 kPa，殘余強度分別增加22.6、23.9、27.5 kPa，峰值位移分別增加18.6、27.7、30.6 mm。但當球度為1.0時，法向荷載對峰值應力、殘余強度與峰值位移的影響均較小，雖然隨著法向荷載增大三者均增加，但增加趨勢較小且呈無規律遞增。如S為1.0時，當法向荷載由100 kPa增加到400 kPa，峰值應力分別增加22.9、22.8、23.5 kPa，殘余強度分別增加10.8、7.9、16.7 kPa，峰值位移分別增加10.8、3.7、7.7 mm，增加趨勢缺乏明顯規律性。

圖5 不同法向荷載下3種球度試樣的剪切位移-剪應力曲線

圖6 不同法向荷載下3種球度試樣的峰值應力、殘余強度和峰值位移變化規律

以上分析說明法向荷載和球度均對粗粒土抗剪強度特性有一定的影響，且隨著球度和法向荷載的變化，峰值應力和殘余強度以及峰值位移均呈現出有規律的變化。

3.2 顆粒形狀對粗粒土黏聚力和內摩擦角的影響

由上述結論，得到3種不同球度的粗粒土峰值應力和法向荷載的強度擬合直線，如圖8所示。由圖8可以看出，不同球度粗粒土峰值應力和法向荷載之間存在很好的線性關系，結合摩爾-庫倫定律,并應用直線擬合方法，可以得到試樣的內摩擦角φ和黏聚力c，其中φ和c的值為通過擬合直線的截距y和斜率α得到，如表4。黏聚力c和內摩擦角φ與球度S的關系曲線見圖9。

表4 不同球度的內摩擦角和黏聚力

從圖9可以看出，隨著球度的變化，黏聚力和內摩擦角都呈現出有規律的變化。當顆粒球度減小時，內摩擦角φ和黏聚力c均增加，內摩擦角φ增加速度較均勻，如球度從1.0到0.726，內摩擦角增加7.82°；從0.726到0.446，增加8.84°，兩者遞增幅度相差較小。但黏聚力c隨著球度的減小而遞增的趨勢趨于平緩，球度從1.0到0.726，黏聚力增加7.57 kPa;但從0.726到0.446，黏聚力增加0.47 kPa，類長條型和類橢圓型顆粒試樣的c值與圓球型顆粒試樣的c值相比相差較大。

4 室內直剪試驗

在進行數值模擬的同時，設置了4種法向荷載下3種球度的粗粒土室內直剪試驗，進行雙向驗證。

4.1 試驗材料及試驗方法

為明確反應顆粒形狀對粗粒土抗剪強度的影響，并同時對數值模擬結果進行驗證，基于上文中選取的3種真實顆粒(圖1)，通過拋光打磨，并篩分后選取粒徑為15～20 mm的石英砂巖顆粒。為了更好地與數值模擬結果進行對比，自制了尺寸大小為500 mm×500 mm×500 mm的可視化剪切盒，并按數值模擬指定的級配制作試樣。試驗過程中，下剪切盒固定不動，上剪切盒在水平面內移動，試驗剪切速度控制在0.1 mm/s，剪切至200 mm停止。通過對真實粗粒土的顆粒大小、形狀、巖性及級配以及剪切盒大小、剪切速率、位移等進行一系列選取，保證室內試驗中顆粒表面粗糙度、試樣大小和密度以及試樣空隙度等和數值模擬基本一致，球度雖然會存在差異，但平均差值在0.02以內，對整體結論無任何影響。本文設置了3組不同顆粒形狀的室內直剪試驗，每組試驗包含的法向荷載為100、200、300、400 kPa。

4.2 試驗結果與數值模擬結果對比

如圖10所示為室內試驗得到的剪切位移-剪應力曲線與數值模擬結果兩者之間的對比。由圖10可知，在4種不同法向荷載下，室內試驗曲線達到峰值時的剪切位移較模擬試驗均有所提前，且室內試驗得到的強度也均略低于模擬試驗得到的結果。

圖7 球度相同時不同法向荷載下的峰值應力、殘余強度和峰值位移變化規律

圖8不同球度粗粒土的強度擬合直線 圖9粗粒土黏聚力、內摩擦角和球度的關系

這主要是由于粗粒土在室內試驗的剪切過程中會產生破碎，且試樣的夯填密度和均勻度都不如模擬狀態好，因此與數值模擬結果對比有一定的偏差。但整體而言，室內試驗曲線在形態上與模擬曲線擬合度較好，尤其是在應力上升及應力下降階段，其上升與下降幅度相當，證明數值模型與試樣在抗剪強度的形成和破壞上具有相同的變化規律，驗證了模型的正確性，由此也說明了數值模擬采用的細觀參數是合理的，生成3種不同類型的顆粒單元來模擬粗粒土形狀的模擬方法也是可行的。表5為不同法向荷載下室內試驗得到的抗剪強度參數。

由表5得到室內試驗中隨著球度的變化粗粒土峰值應力、殘余強度和峰值位移的變化規律，以及球度一定時，法向荷載發生變化時三者的變化規律，同時得出內摩擦角φ和黏聚力c的變化規律，并和數值模擬得到的結果進行對比。如圖11～13所示，具體過程前面已有描述，在此不再贅述。

表5 室內試驗峰值應力-位移與殘余強度表

圖10 不同法向荷載下3種球度粗粒土室內試驗剪切位移-剪應力曲線和數值模擬對比圖

圖11 室內試驗和數值模擬峰值應力、殘余強度和峰值位移變化規律對比圖

圖12 球度相同時不同法向荷載下峰值應力、殘余強度和峰值位移變化規律對比

由圖11～13可以看出，隨著球度以及法向荷載的變化，由室內試驗得出的參數變化趨勢與數值模擬的結果基本相同，均呈現出有規律的變化：當粗粒土球度減小時，其抗剪強度參數(峰值應力、峰值位移和殘余強度)和細觀力學參數(黏聚力c和內摩擦角φ)均增大，且增大趨勢和相互之間的差異與數值模擬的結果基本相似。同時，隨著法向荷載的變化，其參數的變化規律，相互之間的差異，與數值模擬得到的結果擬合度也較好。整體而言，由室內試驗得到的一系列結果充分證明了數值模擬方法的可行性以及結論的正確性。

4.3 細觀機制探討

通過數值模擬和室內試驗得出的結果以及相互之間的對比，從細觀機理角度對以上現象進行分析：形狀不規則的粗粒土，在直剪試驗中當剪應力越過峰值后，由于顆粒發生移動、轉動甚至翻越相接觸的顆粒都比較強烈，所以應變軟化現象較明顯。隨著球度S的減小，粗粒土的顆粒形狀越不規則，顆粒間的咬合自鎖作用越強烈，顆粒運動所需克服的阻力也越大，從而粗粒土可以獲得更高的峰值應力和更大的峰值位移。同時，因為球度的減小使顆粒之間的相互運動比較困難，導致顆粒之間的咬合承擔了部分能量，所以殘余強度也較大。當法向荷載增大時，由于外力作用使顆粒間得咬合更加強烈，顆粒運動需要克服更大的外在阻力，所以以上現象更加明顯。

圖13 室內試驗和數值模擬黏聚力、內摩擦角變化規律對比

內摩擦角反映的是粗粒土的摩擦特性與抗剪強度，粗粒土顆粒之間嚴格來講并沒有黏聚力c的存在，這里的c值更多的是反映顆粒之間的咬合自鎖作用。隨著球度的減小，顆粒形狀變的更加不規則，顆粒之間相互咬合自鎖作用愈徹底，因此內摩擦角和黏聚力愈大。

5 結 論

(1)粗粒土的顆粒形狀對其抗剪強度特性有很大影響。相同法向荷載下，球度越小的粗粒土，剪應力越過峰值應力后的應變軟化現象越明顯。同時隨著球度的減小，粗粒土的峰值應力和殘余強度以及峰值位移均逐漸增加，但增長趨勢趨于平緩。

(2)顆粒形狀對粗粒土抗剪強度的影響與法向荷載有關。球度相同時，隨著法向荷載的增大，粗粒土的峰值應力和殘余強度以及峰值位移均呈規律性遞增,而球度為1.0時增加趨勢無明顯規律。

(3)顆粒形狀對粗粒土的黏聚力和內摩擦角有重要影響。隨著球度的減小，內摩擦角和黏聚力均顯著增加，內摩擦角增長速度較為均勻，而黏聚力增長速度越來越小并趨向于平緩，圓球型顆粒的黏聚力明顯小于其他兩種不規則顆粒。

(4)總的來說，3種顆粒試樣中，2種非圓球型顆粒由于其形狀不規則，在其他細觀參數相同的前提下，比圓球型顆粒展現出了更強的咬合自鎖作用，從而使其抗剪強度增大，參數變化規律最明顯，其中類長條型顆粒的表現最為顯著。