基于離散元法的砂石混合體直剪試驗結果分析

蘇永華，王 棟

（湖南大學土木工程學院，湖南 長沙 410082）

對于地表高程較大甚至懸殊區域，在建、構筑物施工前都要進行一定的或較大、甚至大規模的回填處理，這種狀況下建、構筑物下的地層就稱為重塑地層，砂石地層是常見的重塑地層之一，在基坑、路基等工程中廣泛存在。砂石混合體由力學性質以及結構相差極大的材料組成，是典型的土石混合體。砂石混合體構成的地層由于顆粒粒度變化較大、分布不均勻、塊石形態各異等原因，在實際工程中易出現沉陷或塌陷等問題[1?2]，因此，砂石混合體力學性質研究是亟需進行的研究課題之一。

室內試驗是研究土石混合體力學性質的基本方法，考慮不同影響因素對宏觀力學性能的影響是研究的重點。但室內試驗往往受尺寸效應和試驗條件等因素的影響，存在可重復性差、耗時較久、不易對顆粒形狀進行分析的限制。隨著計算機的發展，顆粒流模擬（PFC）成為研究土石混合體力學性質的一種方法，它能夠模擬微觀層面的結構，充分考慮砂石混合體離散性的特點，從不同角度探討力學性質，如：顆粒實際形狀、顆粒轉動、力鏈演化和孔隙變化等[3]。

礫石形狀是砂礫石力學特性研究的重要屬性參數，目前許多學者采用規則圖像對礫石進行描述，如圓形、正多邊形，沒有考慮礫石的真實形狀，很難反映出礫石真實的力學性質。許多學者利用數字圖像處理技術獲取真實的細觀結構模型，如徐文杰[4]、廖秋林等[5]通過實際拍攝的數碼圖像生成土石混合體結構模型。還有學者通過識別自然塊石的輪廓構建塊石數據庫，實現塊石二維[6]、三維[7]隨機模型建立，數字圖像處理技術可以在不忽略礫石真實形狀的前提下，對形狀參數進行細致化分析，探究其影響作用。

粒度分布影響著土石混合體的變形、強度等物理力學特性。在現有的巖土工程散體材料的粒度分布描述中，大多是采用特征粒徑、不均勻系數和曲率系數等參數來描述粒度分布，這些指標只能表示某些關鍵點的信息，缺乏對粒度整體分布的描述。針對上述情況，杜修力等[8]通過對多組試樣的研究，發現土石混合體的粒度分布均滿足分形結構，可以利用粒度分維值定量描述。舒志樂等[9]建立了土石混合體的二重分形結構模型，對土石混合體的粒度分維值與抗剪強度關系進行分析。

本文通過數字圖像處理技術獲取礫石邊界信息，建立礫石形狀數據庫。結合分形理論建立了砂石混合體二重分形結構模型。通過顆粒流軟件對不同粒度分布、顆粒形狀的砂石混合體進行直剪試驗，并對結果進行分析，為進一步研究砂石混合體重塑地層塌陷機理奠定基礎。

1 砂石混合體隨機模型構建關鍵步驟

1.1 基于礫石形狀數據庫構建離散元模型

通過數字圖像處理技術獲取礫石的邊界坐標和形狀參數，并導入PFC 中生成顆粒簇模板以供模型構建時調用，詳細步驟如下：

（1）圖像獲取及預處理：通過高精度數碼相機獲取大量礫石圖像，但是礫石因光線等環境因素圖像會出現大量噪音，需要利用圖像處理軟件對圖像進行消噪，并使每張圖片僅顯示1 個礫石（圖1a）。

圖1 數字圖像處理技術獲取邊界信息流程示意圖Fig.1 Flow chart showing the digital image processing technology to obtain boundary information

（2）邊緣提取：預處理后的圖像需進行二值化處理，可通過設定合適閾值，使圖像變得簡單，凸顯感興趣目標輪廓，處理后圖像呈現明顯的黑白效果（圖1b）。二值化圖像可通過邊緣檢測獲取礫石邊界，本文利用MATLAB 軟件中內嵌的Canny 算法提取礫石邊界（圖1c）。

（3）邊界擬合：通過上述方法提取的邊界點數目過多，會對后續工作產生較大負擔，需要對邊界進行擬合，擬合方法如下[6]：選定初始點，然后順時針獲取所有邊界點，依次計算邊界點和初始點的距離，如果該距離滿足式（1）則保留該點為擬合后的邊界點，并作為下一輪的初始點，依此循環直至回到初始點。

式中：d—擬合后兩個相鄰邊界點的距離；

D—與礫石顆粒面積相等的等效圓直徑；

α—控制邊界擬合精度的參數。

通過調整 α可以控制礫石邊界的擬合精度。α越大所需邊界點越少，但精度越低。本文選取 α為0.1，此時礫石的形狀和原圖像基本一致（圖1d），坐標點的數量也滿足后續的數值模擬需求。

（4）圖像標準化處理：為了方便后續對礫石大小、旋轉角度等參數的控制，需進行標準化處理，如粒徑統一化、形心歸零化、長軸水平化。

（5）PFC 顆粒簇模板生成：獲取的礫石邊界信息不能直接被PFC 識別，需要將其轉換為矢量數據格式文件（dxf）。本文利用Python 第三方庫dxfwrite 完成批量轉換，導入PFC 并利用內置的顆粒填充算法，就可生成礫石的顆粒簇模板。

1.2 形狀參數分析

顆粒形狀是研究砂石混合體宏觀力學行為的一個重要參數。根據相關研究[10]，顆粒可以通過軸向系數和圓度進行描述，其中軸向系數S為：

式中：LA,max，LA,min—礫石顆粒等效橢圓的長軸和短軸。

礫石顆粒的圓度公式定義為R：

式中：L，A—顆粒的周長和面積。

通過礫石形狀數據庫可導出形狀參數，并利用數據分析軟件得出圓度和軸向系數的頻率分布直方圖（圖2），圓度分布在1.0～1.35 之間，均值為1.15，圓度指標可有效地表征顆粒整體形狀接近圓的程度。圓度越大，其形狀越偏離圓形。軸向系數分布在1.0～1.8之間，均值為1.28，軸向系數越大，礫石的長軸與短軸比值越大，條狀性越明顯。

圖2 圓度頻率（a）和軸向系數頻率（b）分布直方圖Fig.2 Distribution histogram of(a)roundness frequency and(b)axial coefficient frequency

1.3 砂石混合體的二重分形結構模型

Tyler 等[11]提出了基于分形理論的分形結構模型：

式中：r—顆粒粒徑；

R—某一顆粒粒徑；

MT—總質量；

RL—最大顆粒尺寸；

M(r

D—粒度分維值。

根據式（4）得出，如果在雙對數坐標下lg[M(r

砂石混合體的粒度分布具有良好的分形結構，但是由于其復雜的顆粒組成，通常具有多重分維。本文結合砂石混合體的組成以及土石閾值構建二重分形結構模型:

式中：Ds—砂粒度分維值；

Dr—礫石粒度分維值；

RS/RT—土石閾值。

土石閾值，即“土”和“石”的界限粒徑，本文結合以前學者的研究結果[12?13]，選取5 mm 作為砂和礫石的界限粒徑。令式（5）中R=RS/RT即可確定礫石的最大粒徑，公式如下[14]：

式中：P—含石量。

徐文杰等[15]指出含石量在25%～75%之間時，巖土體的強度取決于其中的“土體”與“塊石”，其力學強度是兩者相互作用的共同反映。含石量在這范圍之外，整體的強度會主要取決于“土體”或者“塊石”。為了分別研究砂和礫石級配變化對土石混合體強度的影響，需考慮兩者共同作用，因此本文選取含石量為50%。

根據式（6）得出最大粒徑、含石量和礫石粒度分維值的關系（圖3），可以看出：隨著含石量和礫石粒度分維值的增加，最大粒徑變大；最大粒徑保持不變時，礫石粒度分維值的增加，會導致含石量減少，兩者關系曲線呈近似拋物線形狀。

圖3 最大粒徑、含石量和礫石粒度分維值的關系Fig.3 Relationship among the maximum grain size,rock content and fractal dimension of gravel size

在確定含石量和土石閾值后即可根據不同粒度分維值確定砂石混合體的完整級配分布。首先需要確定每組粒徑區間，取區間的粒徑上限為該組粒徑，并且粒組按粒徑由大到小的順序，即第1 個粒組的粒徑為礫石最大粒徑，則第i個粒組含量占總試樣的百分數：

式中：R(i)—礫石第i組粒徑上限。

將粒度分維值Dr或Ds代入式（7）就可確定不同粒度分維值下砂石混合體的試樣級配。

2 基于PFC 的直剪試驗模擬

2.1 直剪試驗模擬流程

（1）剪切盒建立：砂和礫石顆粒單元均需要在剪切盒內生成。剪切盒尺寸和室內大型直剪試驗截面一致，為300 mm×150 mm。模擬過程中剪切盒被認為是剛性體，剪切盒的剛度需大于砂和礫石的剛度，保證顆粒不飛出，并且避免模型內部產生較大初始接觸力。通過試算，本文剪切盒剛度選取為1×1010N/m 。

（2）試樣生成：砂采用ball 模擬，礫石采用前述構建的顆粒簇模板生成，接觸模型采取線性接觸模型。顆粒級配利用確定的二重分形結構模型，計算每個粒組需要生成顆粒的數量，并采用粒徑擴大法生成特定級配的試樣。需要注意的是：如果顆粒大小按照1∶1 的比例進行模擬，會對計算機提出較高要求，造成效率的降低。本文通過試算，設置顆粒放大系數為1.6，并忽略0.5 mm 以下顆粒，減少生成顆粒的數量，從而平衡PFC 模擬的精度和效率。

（3）施加豎向荷載：顆粒流軟件中墻體無法直接施加力的作用，需要利用墻體的伺服機制，即在每一時步調整墻體的速度來控制墻體的受力大小。通過施加豎向荷載形成均勻試樣，待穩定后，便可進行試樣剪切。本文選取的主要正應力為100，200，300 kPa。

（4）試樣剪切：剪切過程中維持上剪切盒不動，下剪切盒以一定速度向右移動，當剪切位移達到30 mm時停止剪切。剪切過程中監測剪應力和剪切位移的變化，記錄數據。

2.2 細觀參數標定

PFC 中顆粒的細觀參數不能直接通過宏觀參數獲得，需要通過室內試驗獲取應力-應變曲線，并通過不斷改變細觀參數，使得數值模擬的結果與試驗結果相近，完成細觀參數的標定。本文采用Geotest S2450 型摩擦剪切儀進行直剪試驗，礫石和粗砂作為試驗材料。

由于砂石混合體參數較多，因此本文首先通過對單一砂進行直剪試驗，獲取砂的細觀參數。并以此為基準進行砂石混合體直剪試驗，完成其余參數的標定。經過大量參數試算，最終模擬曲線和實際試驗曲線良好（圖4），具體參數見表1。

表1 基于離散元模擬的直剪試驗主要計算參數Table 1 Main computational parameters in DEM simulation of direct shear

圖4 砂（a）和砂石混合體（b）剪應力-剪切位移模擬和實際試驗結果對比曲線Fig.4 Shear stress-shear displacement comparison curve between the simulation and actual test results of(a)sand and(b)sand-gravel mixture

2.3 剪應力-剪切位移曲線特征分析

砂石混合體的抗剪強度明顯高于單一砂時的抗剪強度。隨著正應力從100 kPa 增加至300 kPa，顆粒間受到擠壓排列更加緊密，剪應力逐漸增加。不同正應力下剪應力-剪切位移曲線呈明顯的非線性特征，可分為3 個階段[16]：線彈性階段、初始屈服階段以及應變硬化階段，剪切過程中沒有出現明顯的峰值。

3 砂石混合體直剪試驗結果分析

3.1 粒度分維值的影響分析

砂石混合體的抗剪強度分別受砂粒度分維值和礫石粒度分維值影響。為單獨研究礫石粒度分維值的影響，采用控制變量法進行研究，首先保持砂的粒度分維值為2.55 不變，礫石的粒度分維值分別取：

2.45，2 .50，2.55，2.60，2.65。

砂石混合體的粒度分布曲線按照前文構建的二重分形結構模型計算，圖5 為砂石混合體級配示意圖，為了方便觀察僅列出3 條典型曲線（下同），可以看出隨著Dr增大，最大粒徑也增大，在Dr=2.65 時，最大粒徑達到32.23 mm。

圖5 不同礫石粒度分維值的級配曲線Fig.5 Grading curves of fractal dimension of different gravel sizes

圖6 是相應的剪應力-剪切位移曲線，采用了逐步平均法進行平滑，可以發現3 條曲線都近似呈應變硬化型，在Dr=2.55 時，曲線末尾其剪應力上升趨勢還十分明顯。抗剪強度與礫石的粒度分維值不是簡單的線性變化，Dr為2.55 時抗剪強度最大，其值為220.05 kPa，而在Dr=2.45 和2.65 時，其值分別為196.07 kPa 和195.48 kPa，均明顯小于Dr=2.55 的抗剪強度。

圖6 300 kPa 正應力下不同Dr 剪應力-剪切位移曲線Fig.6 Shear stress-shear displacement curves of different Dr under the normal stress of 300 kPa

依據摩爾庫倫準則τ=σtanφ+c，可以進行抗剪強度參數的計算。砂石混合體由砂和礫石組成，兩者的黏聚力非常小，所以本文不考慮黏聚力，取值為0 MPa，此時抗剪強度主要由內摩擦角產生（表2）。

表2 不同礫石粒度分維值在直剪試驗中獲取的相關參數Table 2 Correlation parameters of fractal dimension of different gravel sizes obtained in the direct shear test

采用同樣方法研究砂粒度分維值的影響，保持礫石的粒度分維值為2.55 不變，改變砂的粒度分維值，最終得出不同正應力下抗剪強度及內摩擦角的變化（表3）。

表3 不同砂粒度分維值在直剪試驗中獲取的相關參數Table 3 Relevant parameters obtained in the direct shear test of the fractal dimension of different sand particle sizes

為了方便觀察，以砂、礫石粒度分維差值作為指標進行分析（圖7），明顯看出內摩擦角與Dr、Ds不是簡單的線性變化，而是呈現先增大后減小的趨勢。隨著砂、礫石粒度分維值的接近，內摩擦角逐漸增加，當Dr?Ds=0 即砂石混合體具有單一分維時，內摩擦角達到最高值36.34°。

圖7 粒度分維差值與內摩擦角關系Fig.7 Relationship between the difference in particle size fractal dimension and the angle of internal friction

為了探究影響原因，通過PFC 繪制出剪切結束時顆粒位移圖（圖8），顆粒的位移主要由兩部分造成：第一部分是由下剪切盒的移動造成的，剪切過程中，上剪切盒保持不動，下剪切盒向右移動了30 mm；第二部分是由于剪切過程中顆粒的相互錯動造成的位移。剪切盒的右下角顆粒位移最大，有的超過了30 mm，這主要是由于剪切過程中顆粒的移動并不能與剪切盒同步，在短時間里顆粒右側形成一段空隙，使其變得較為松散，顆粒接觸減少，易產生錯動位移。由于上剪切盒在剪切過程中沒有移動，位移主要由顆粒之間相互作用影響，因此具體分析上剪切盒顆粒位移。分析10～15 mm 剪切帶附近的顆粒位移可以看出，隨著砂、礫石粒度分維值的接近，10～15 mm 位移的顆粒范圍在不斷擴大，未受影響的顆粒范圍在減少。這說明隨著砂、礫石粒度分維值的接近，可觀察到顆粒排列的越緊密，因剪切盒移動產生的相互錯動越明顯，顆粒位移增加，導致抗剪強度增加，內摩擦角也增大。

圖8 不同粒度分維差值的位移變化圖Fig.8 Displacement change diagram of difference of different fractal dimensions

綜上，當砂、礫石粒度分維值越接近，其抗剪強度越大，內摩擦角也相應越大。當兩者相等時，即砂石混合體只具有單一分維，此時試樣的均一性達到最好，具有最大的抗剪強度和內摩擦角。

3.2 軸向系數的影響分析

由于本文在邊界平滑和顆粒填充算法中對礫石細部邊界信息有所缺失，因此只對軸向系數進行具體分析。從礫石數據庫選取顆粒圓度大致相同，軸向系數分別為1.0，1.4，1.8 的礫石進行模擬（圖9）。為了減少其他方面影響，模擬時選取單一顆粒形狀。由圖9可以看出，軸向系數越大，顆粒條狀性越明顯，反之，越接近于圓形。

圖9 不同軸向系數的礫石Fig.9 Gravel with different axial coefficients

對不同正應力得到的剪應力峰值進行線性擬合（圖10），可以發現擬合相關系數R2都在0.99 以上（表4），表明通過摩爾庫倫強度準則得到的內摩擦角比較可靠。隨著軸向系數的增加，抗剪強度和內摩擦角逐漸增加，S=1.8 時內摩擦角達到最高值35.92°。產生這種情況主要是因為，礫石軸向系數越大，顆粒條狀性越明顯，從而有更多機會接觸不同顆粒，顆粒的抗轉動能力也得到了提高[17]。

表4 不同軸向系數在直剪試驗中獲取的強度參數Table 4 Strength parameters obtained in direct shear test with different axial coefficients

圖10 不同軸向系數下剪應力峰值-正應力曲線Fig.10 Shear stress peak-normal stress curve under different axial coefficients

由圖11 可以發現，顆粒在剪切面附近旋轉角度偏高，并且越靠近左右邊緣旋轉角度越大。由于顆粒移動的延遲性以及不均勻導致剪切盒右下方礫石周圍的接觸數量也明顯降低。

圖11 礫石顆粒旋轉角度及周圍接觸示意圖Fig.11 Schematic diagram of rotation angle and surrounding contact of gravel particles

為具體了解相關變化，通過Fish 語言獲取300 kPa正應力下礫石和周圍顆粒的接觸數量以及旋轉角度（表5）。由表5 可以看出，隨著軸向系數從1.0 至1.8變化，接觸數量從1 581 增加至1 742，礫石的平均旋轉角度從9.43°減少至6.42°。這說明隨著軸向系數增加，礫石顆粒的條狀性越強，和其它顆粒的接觸可能性越多，也導致其在旋轉時遇到的阻力越多，抗轉動能力較強，在剪切過程中旋轉的角度總體上偏低，從而導致抗剪強度提高。

表5 不同軸向系數下礫石接觸數量和旋轉角度Table 5 Gravel contact number and rotation angle under different axial coefficients

4 結論

（1）利用數字圖像處理技術能夠自動化完成礫石真實圖像至顆粒簇模板的轉換，構建礫石的形狀數據庫，并可對軸向系數、圓度等形狀參數進行分析。

（2）基于分形理論建立了砂石混合體的二重分形結構模型，得出完整的粒度曲線。土石閾值一定時，砂石混合體的最大粒徑隨含石量、礫石粒度分維值的增加而變大。最大粒徑保持不變時，礫石粒度分維值的增加，會導致含石量減少，兩者關系曲線呈近似拋物線形狀。

（3）砂、礫石粒度分維值越接近，抗剪強度越大，內摩擦角也相應越大，當兩者相等時，即砂石混合體只具有單一分維時，此時試樣的均一性達到最好，具有最大的抗剪強度和內摩擦角。

（4）軸向系數S是形容礫石形狀的一個重要參數，其分布在1.0～1.8 之間。隨著軸向系數S的增加，礫石條狀性越強，在直剪試驗中抗轉動能力越強、與周圍接觸可能性越大，導致抗剪強度和內摩擦角不斷增加。