綜合考慮堆石體復雜特性的顆粒流高精度高效率模擬新方法

趙佩臣，劉玉鳳

（1.中國電建集團中南勘測設計研究院有限公司，長沙 410007；2.42082119××××295525，長沙 410007）

1 引言

基于離散元的顆粒流方法（PFC）能從細觀層面對堆石體的力學性質進行研究，在堆石體的模擬上應用日益廣泛。近年來，國內外學者運用PFC 對堆石體的力學特性展開了大量的研究工作，但這些研究仍存在一定的不足之處。首先，模擬的精度有待進一步提高，PFC 將顆粒簡化為理想化的剛性圓形，未考慮堆石體的不規則形狀和破碎效應。另外，模擬的效率也有待進一步提高，主要表現為建模十分耗時。

本文利用PFC 5.0 程序，針對當前堆石體模擬存在的精度和效率問題，提出了一種高效率高精度方法，可以快速建立計算模型，并且真實再現堆石體在加載條件下的宏細觀力學行為。

2 綜合考慮堆石體復雜特性的顆粒流模擬方法

堆石料具有典型的級配特性，其中粗粒料為隨機不規則形狀，加載時會發生破裂。本文利用PFC 5.0 中的2D 模塊進行模擬計算，并以糯扎渡I 區[1]堆石料為例，對在顆粒流模型中如何模擬堆石體的破裂性質進行具體說明。

在PFC 5.0 中的2D 模塊，可根據級配在模型區域直接生成初始模型，模型平衡后記錄所有顆粒的形心坐標和半徑，隨后將生成的顆粒全部刪除；在相同的區域以較細的粒徑再次生成計算模型，平衡后由上述記錄的形心坐標和半徑確定相關區域，將落入該區域的小顆粒黏結在一起構成顆粒簇。假設某顆粒簇由N 個細顆粒組成，細顆粒的粒徑為ri（i=1，…，N），計算模型的空隙率為n，則該顆粒簇的等效粒徑R 的計算方法如下：

圖1 是截斷級配中的細小粒徑后，根據上述方法建立的初始模型和計算模型，并反映了二者的級配關系。從圖1 可以看出，除細粒料的含量一定差別外，計算模型的級配與堆石體的原始級配基本一致。

按照本文提出的方法，可直接在指定區域內生成密實的顆粒模型，這是該方法高效性的第一個體現。在建立模型的過程中，雖然需要2 次生成顆粒并使之平衡，但是可采用時步放大法可使模型快速達到平衡狀態。

圖1 PFC 2D模塊數值模型及其級配曲線

初始模型中的剛性圓形顆粒在計算模型中由若干細顆粒黏結在一起進行等效模擬。從計算模型的視圖可以看出，雖然模擬粗粒料的顆粒簇是將某圓形區域內的細小顆粒黏結在一起構成的，但是該顆粒簇的輪廓形狀呈現出不規則性；由于落入圓形區域內的小顆粒是隨機的，因此，這些小顆粒所形成的顆粒簇的幾何形狀具有隨機性；顆粒簇由若干小顆粒黏結在一起，當黏結力達到接觸的強度時，黏結會發生斷裂，顆粒簇出現局部或整體的破碎，可模擬堆石料的破碎特性。

從上述粗粒料等效替換的過程中可以看出，將初始模型中的顆粒位置和粒徑等信息提取出來，用在計算模型中建立對應的等效顆粒簇，可以使二者的粒徑基本一致，還可以直接指定等效顆粒簇的位置，這是該方法高效性的第二個重要體現。

對于顆粒流模型而言，無論顆粒簇的形狀如何復雜，其相互間的接觸關系始終為簡單的顆粒-顆粒接觸，因而在計算過程中無須對接觸類型進行識別，這是顆粒流方法相對于SGDD方法[2]最大的優勢，可顯著提高計算速度。

綜上所述，本文提出的模擬方法，可以綜合體現堆石體的級配、隨機不規則形狀及破碎效應，真實地反映堆石體的細觀特性，具有很高的效率。

3 綜合考慮堆石體復雜特性的顆粒流建模步驟

綜合考慮堆石體特性的顆粒流建模步驟主要包括以下階段：

1）根據顆粒級配在模型區域內生成初始模型，記錄所有顆粒的形心坐標和半徑，并刪除全部顆粒；

2）設計一組合適的粒徑dmin和粒徑比Rratio，再次在模型區域生成顆粒，孔隙率與步驟1）中保持一致，同時在模型中摻入一定含量粒徑小于dmin的細小顆粒；

3）根據步驟1）中記錄的顆粒形心坐標和粒徑，在步驟2）中生成的計算模型里確定對應的區域，將落入該區域的顆粒黏結在一起，構成顆粒簇，并賦予模型細觀力學參數；

4）刪除步驟2）中摻入的細小顆粒，生成最終的計算模型。

4 基于該方法的糯扎渡I 區堆石體壓縮模擬

為了驗證本文提出方法的精度，本文針對糯扎渡I 區堆石料進行壓縮模擬，并將模擬的數值與物理試驗進行對比分析。

4.1 計算模型及力學參數

本文將模型的尺寸設計為3.25m（寬）×6.50m（高），同時將級配中小于0.1m 的粒料按照面積等效的原則由d=0.1m 的粒料進行替換，計算模型中的最小粒徑為dmin=0.3d，粒徑比Rratio=1.66，空隙率0.21，摻入的細小顆粒粒徑為dsmall=0.6dmin，最終建立計算模型中顆粒總數為3 215，計算模型及其級配如圖1 所示。

計算模型中粗粒料由顆粒簇模擬，簇內部顆粒間的接觸關系為平行黏結接觸，細顆粒間、細顆粒與顆粒簇間、相鄰的顆粒簇間以及顆粒與墻體的接觸關系均為線性接觸。

根據物理試驗，堆石料E-B 模型計算參數如表1 所示，反算出堆石料在不同圍壓下的應力應變曲線。

表1 糯扎渡I 區堆石料E-B模型計算參數

經過大量試算，最終確定了一組比較合適的細觀力學參數，如表2 和表3 所示。

表2 線性接觸細觀參數

表3 平性接觸細觀參數

模型的頂部墻體加載速度為0.01m/s，底部墻體固定，計算時步取8×10-7s/step，則模型實際的加載速度為8×10-6mm/step。

4.2 模擬結果對比及分析

根據建立的計算模型和設置的接觸模型及細觀力學參數，模擬堆石體在不同圍壓下的壓縮試驗，模擬時設置4 組不同的圍壓，分別為 σ3=300kPa、800kPa、1 400kPa 及 2 000kPa。

1）應力應變關系對比。如圖2 所示，經過模擬計算得到的各組應力應變曲線與反算曲線十分接近。

2）體積變形特征分析。如圖3 所示，通過與物理試驗結果[3]進行對比可以發現，各組數值模型與物理試樣的體積變形規律十分接近。

圖2 數值模型與E-B模型的應力應變曲線

圖3 數值模型體應變- 軸應變關系

3）堆石體宏觀破壞特征。在不同的圍壓下，當加載結束時，計算模型中顆粒的位移分布如圖4 所示。通過圖4 可以看出，在低圍壓下，模型內部出現了一定范圍的宏觀剪切帶；隨著圍壓不斷增加，這種剪切破壞特征逐漸減弱，在中高圍壓下，模型中沒有出現明顯剪切帶。

圖4 不同圍壓條件下模型位移分布

物理試驗[4]結果表明，堆石體在加載至破壞時，其宏觀破壞類型為壓鼓，試樣內部沒有出現明顯剪切帶。而文獻[5]對堆石體進行模擬時，中高圍壓下模型內部出現了明顯的X 形剪切破碎帶，顯然這與物理試驗觀察的結果不符。通過本文的模擬結果可以看出，低圍壓條件下模型內部出現一定程度的剪切帶，隨著圍壓升高，模型內部顆粒間的運動特征以破碎壓實為主，沒有發生大量錯動形成宏觀剪切帶。另外，堆石體破裂形成的裂紋隨機分布在模型之中，并非集中在模型的對角區域，這與模型內部沒有出現明顯剪切帶的結論是相互印證的。由此可見，本文的模擬方法能更加真實地體現堆石體在壓縮條件下的宏觀破裂特征。

4.3 模擬結果分析

通過從不同方面進行對比可以看出，本文的模擬結果與物理試驗結果比較接近，真實地反映了堆石體的宏觀力學特性，模擬精度高。

5 主要結論

針對當前顆粒流方法模擬堆石體時面臨的精度和效率問題，本文提出了一種新的模擬方法，并運用該方法對糯扎渡I 區堆石料的壓縮試驗進行模擬。主要結論包括以下幾個方面：

1）本文提出的方法可快速建立計算模型，并綜合考慮堆石體的級配、隨機不規則形狀和破碎效應等特性，極大地提高了顆粒流方法在模擬堆石體時的建模效率；

2）本文的模擬結果真實地再現了堆石體在不同圍壓下的宏觀力學行為，數值模型的應力應變曲線、體積變形規律以及宏觀破壞特征與物理試驗的結果十分接近，充分說明了本文提出的方法具有較好的模擬精度；

3）該模擬方法不僅可以用于模擬小尺度的室內試驗，從而進一步推動顆粒流方法在實際工程問題上的應用。