本發(fā)明涉及一種模擬顆粒和堆積試驗模擬方法，特別是一種離散元組合顆粒及其離散元堆積試驗模擬方法。

背景技術：

離散元方法與有限元法等連續(xù)介質力學方法相比可以考慮土顆粒的散體特征，克服了宏觀連續(xù)性假設，能夠很好的用于求解土體非線性等問題，通常用于模擬散體介質的力學特性。對于散體介質——砂土的物理力學特性研究，現有的研究主要集中在顆粒粒徑和級配對砂土宏觀特性的影響，而顆粒粒徑、級配及形狀是影響砂土力學及變形特性的三個極其重要的參數，傳統(tǒng)的離散元方法中通常采用圓球形顆粒模擬砂土，與實際砂土的不規(guī)則形狀差別較大，離散元的圓球形顆粒很難實現顆粒的互鎖現象，即圓球顆粒難以實現不規(guī)則形狀天然砂土的咬合摩擦特性，通常采用增大摩擦系數使最后的模擬結果逼近室內試驗，最終導致設置的摩擦系數遠遠超過實際值。此外，采用圓球顆粒模擬實際的砂土力學特性，受力過程中顆粒的咬合力與旋轉角度明顯小于非圓顆粒，與實際各種形狀組合的土體顆粒受力機理差別較為顯著，難以準確的從顆粒細觀角度揭示砂土的受力機理。

技術實現要素：

本發(fā)明要解決的技術問題是：針對目前離散元純圓顆粒模擬方法的不足，提供一種新的組合離散元組合顆粒及堆積試驗模擬方法，能夠更準確的模擬實際試驗砂土顆粒堆積試驗，有效反映砂土力學特性及實際顆粒特征，減少現有模擬試驗工作量，彌補傳統(tǒng)圓形顆粒模擬的不足，提高模擬的準確度。

解決上述技術問題的技術方案是：一種離散元組合顆粒，其特征在于：包括各種球形度的Clump顆粒，所述的Clump顆粒為由單元球A、單元球B、單元球C和單元球D四個單元球構成，所述的單元球A、單元球B、單元球C的半徑相等，所述的單元球A、單元球B、單元球C以圓周分布方式排列于單元球D的四周，且單元球A、單元球B、單元球C均與單元球D相切或相交，所述Clump顆粒的球形度n的取值為0.33≤n≤1。

本發(fā)明的進一步技術方案是：所述的離散元組合顆粒為組合Clump顆?；騿我籆lump顆粒，所述的組合Clump顆粒為包括多種不同球形度的Clump顆粒組合成的有多種外形的顆粒組合，所述的單一Clump顆粒為由1種球形度的Clump顆粒組合成的只有一種形狀的顆粒組合。

所述Clump顆粒的球形度n為0.33，0.4，0.5，0.58，0.70，0.77，0.9或1.0中的1個、2個或多個。

本發(fā)明的進一步技術方案是：一種離散元組合顆粒的堆積試驗模擬方法，包括以下步驟：（1）根據光學顯微鏡拍攝的實際試驗砂土顆粒照片進行黑白二值化處理并提取出該實際試驗砂土顆粒單元的幾何參數，該幾何參數包括球形度；

（2）根據步驟（1）中得到的實際試驗砂土顆粒的幾何參數為基礎，建立一種離散元組合顆粒數值模型，該離散元組合顆粒包括各種球形度的Clump顆粒，所述的Clump顆粒為由單元球A、單元球B、單元球C和單元球D四個單元球構成，所述的單元球A、單元球B、單元球C的半徑相等，所述的單元球A、單元球B、單元球C以圓周分布方式排列于單元球D的四周，且單元球A、單元球B、單元球C均與單元球D相切或相交，所述Clump顆粒的球形度n的取值為0.33≤n≤1；

（3）用離散元PFC3D方法建立上端開口下端閉合的漏斗模型，在漏斗下方處生成底部墻面模擬地面，并設置相應的細觀參數；

（4）根據步驟（1）的參數首先生成對應球形度的離散元組合顆粒模板，然后根據體積等效與質量等效原則在漏斗上方生成與實際試驗砂土顆粒體積一致的離散元組合顆粒，設置離散元組合顆粒細觀參數并施加重力加速度為g＝-9.81m/ s²，然后開始循環(huán)計算，此時離散元組合顆粒在重力作用下逐漸沉積在漏斗內部，設置平衡條件并計算至初始平衡狀態(tài)，在漏斗內部布置測量球測出離散元組合顆粒的平均接觸數及孔隙率；

（5）刪除步驟（4）中漏斗底部的閉合墻體使漏斗中的離散元組合顆粒逐漸下落直至漏斗下方的底部墻面上并形成穩(wěn)定的顆粒堆積體；

（6）通過在步驟（5）中的顆粒堆積體中布置測量球可測得測量球范圍內離散元組合顆粒的平均接觸數及孔隙率，編寫fish函數輸出測量顆粒堆積體表面顆粒坐標以計算自然休止角；

（7）根據步驟（6）中布置在顆粒堆積體中的測量球，可實時獲取離散元組合顆粒堆積體的孔隙率，通過改變顆粒的細觀參數值，可實現求取顆粒堆積體的最大及最小孔隙率，亦可分別研究球形度n、顆粒摩擦系數、孔隙率、接觸數、自然休止角等各參數的相互關系。

通過下述步驟計算自然休止角：（1）在完成步驟（5）中的堆積過程后將離散元組合顆粒球心坐標導出，對非錐頂和非錐腳區(qū)域用最小二乘法對顆粒堆積體表面曲線進行線性擬合，根據擬合直線的斜率求出砂土堆積體的自然休止角；

（2）顆粒堆積體通常呈現不對稱狀態(tài)，為消除方向選擇對休止角測定結果的影響，從顆粒堆積體的x、y坐標軸正負四個方向選取表面坐標，并將擬合結果取均值作為顆粒堆積體的自然休止角。

所述的漏斗通過離散元PFC3D方法分兩段建立漏斗模型，該漏斗的上部采用兩端開口的圓筒墻生成的圓筒結構，圓筒直徑為50mm，高度為55mm；漏斗下部采用上端開口下端閉口的圓錐墻生成的圓錐結構，圓錐上端直徑為50mm，下端直徑為20mm，圓錐高度為15mm，在距漏斗下端50mm處生成底部墻面模擬地面。

所述的離散元組合顆粒為組合Clump顆?；騿我籆lump顆粒，所述的組合Clump顆粒為包括多種不同球形度的Clump顆粒組合成的有多種外形的顆粒組合，所述的單一Clump顆粒為由1種球形度的Clump顆粒組合成的只有一種形狀的顆粒組合。

所述Clump顆粒的球形度n為0.33，0.4，0.5，0.58，0.70，0.77，0.9或1.0中的1個、2個或多個。

由于采用上述技術方案，本發(fā)明之一種離散元組合顆粒及其離散元堆積試驗模擬方法，具有以下有益效果：

（1）本發(fā)明采用離散元方法開發(fā)組合顆粒并進行了堆積試驗模擬，開發(fā)的離散元組合顆粒采取球形度n作為控制顆粒形狀參數，不僅在外形上逼近實際砂土非圓顆粒特征，避免了傳統(tǒng)圓球顆粒無法真實模擬砂土咬合的缺陷，而且可根據顆粒球形度統(tǒng)計結果直接采用實際顆粒的摩擦系數在離散元中模擬砂土顆粒的力學特征，避免了傳統(tǒng)圓球顆粒使用嚴重偏大摩擦系數的缺陷，更逼近實際砂土的力學特性。

（2）本發(fā)明可以方便追蹤砂土顆粒在堆積過程中的孔隙率及接觸數變化，模擬的堆積體自然休止角可與室內實際試驗砂土顆粒自然休止角相吻合，對傳統(tǒng)散體顆粒的細觀模擬研究方法進行了改進，能夠彌補離散元中用圓球顆粒模擬不規(guī)則顆粒的不足，并且還能測出通過室內實際試驗砂土顆粒堆積試驗無法測出的其他參數，例如顆粒摩擦系數、孔隙率、接觸數、等各參數的相互關系，為研究砂土特性的數值模擬提供了新的研究思路。

（3）本發(fā)明根據球形度生成的Clump顆粒試樣最大及最小孔隙比均與室內實際試驗砂土顆粒實測值基本一致，且能同時反應出顆粒間的摩擦作用及嵌鎖效應，對實際試驗或工程進行宏細觀力學機理研究時，可采用球形度作為顆粒形狀統(tǒng)計參數進行數值模擬。

（4）本發(fā)明開發(fā)的離散元組合顆粒，不僅可用于真實模擬砂土堆積試驗，同時亦可用于砂土三軸試驗和直剪試驗，研究砂土在各種工況下的力學特性。利用組合顆粒進行離散元數值試驗不僅能得到室內試驗的宏觀結果，而且還能通過分析組合顆粒的細觀參數變化特征，從而可進一步探討砂土受力特性的宏細觀機理。

下面，結合說明書附圖和具體實施例對本發(fā)明之一種離散元組合顆粒及其離散元堆積試驗模擬方法的技術特征作進一步的說明。

附圖說明

圖1：本發(fā)明之Clump顆粒構造示意圖及多種球形度不同的Clump顆粒對比示意圖。

圖2：本發(fā)明中離散元組合顆粒的堆積試驗數值模擬過程示意圖。

圖2中（a）圖為顆粒堆積試驗圖，（b）圖為測量球布置圖，（c）圖為顆粒堆積體的自然休止角圖。

圖3：本發(fā)明的實施流程圖。

圖4：本發(fā)明的Clump顆粒試樣的自然休止角統(tǒng)計結果。

圖5：本發(fā)明的Clump顆粒試樣的孔隙率統(tǒng)計結果。

圖6：本發(fā)明的Clump顆粒試樣的平均接觸數。

圖7：本實施例中組合Clump顆粒和單一Clump顆粒驗證土樣G的結果比較。

圖7中（a）圖為孔隙率對比圖，（b）圖為平均接觸數對比圖，（c）圖為自然休止角對比圖。

圖8：本實施例中組合Clump顆粒結果擬合曲線。

圖 9：土樣G球形度分布圖

在上述附圖中，各標號說明如下：

1-1為單元球A，1-2為單元球B，1-3為單元球C，1-4為單元球D，1-5最大內切圓，1-6為最小外接圓，1-7為球形度n=0.33的Clump顆粒，1-8為球形度n=0.4的Clump顆粒，1-9為球形度n=0.5的Clump顆粒，1-10為球形度n=0.58的Clump顆粒，1-11為球形度n=0.7的Clump顆粒，1-12為球形度n=0.77的Clump顆粒，1-13為球形度n=0.9的Clump顆粒，1-14為球形度n=1的Clump顆粒。

2-1為漏斗模型，2-2為漏斗頂部開口，2-3為漏斗底部開口墻體，2-4為自然堆積Clump砂土顆粒，2-5為測量球，2-6為砂土自然休止角。

具體實施方式

一種離散元組合顆粒，包括各種球形度的Clump顆粒，所述的Clump顆粒為由單元球A、單元球B、單元球C和單元球D四個單元球構成，所述的單元球A、單元球B、單元球C的半徑相等，所述的單元球A、單元球B、單元球C以圓周分布方式排列于單元球D的四周，且單元球A、單元球B、單元球C均與單元球D相切或相交，所述Clump顆粒的球形度n的取值為0.33≤n≤1。所述的離散元組合顆粒為組合Clump顆?；騿我籆lump顆粒，所述的組合Clump顆粒為包括多種不同球形度的Clump顆粒組合成的有多種外形的顆粒組合，所述的單一Clump顆粒為由1種球形度的Clump顆粒組合成的只有一種形狀的顆粒組合，所述離散元組合顆粒中的Clump顆粒的球形度n為0.33，0.4，0.5，0.58，0.70，0.77，0.9或1.0中的1個、2個或多個[如圖1]，也可以為0.33≤n≤1范圍內任意1個、2個或者多個球形度組合。當球形度n為0.33時所述的單元球A、單元球B和單元球C分別與單元球D相切，當球形度n為1.0時則所述的單元球A、單元球B、單元球C和單元球D四個單元球相互結合成一個圓球形顆粒。

一種離散元組合顆粒的堆積試驗模擬方法，包括以下步驟：（1）根據光學顯微鏡拍攝的實際試驗砂土顆粒照片進行黑白二值化處理并提取出該實際試驗砂土顆粒單元的幾何參數，該幾何參數包括球形度。通過實際試驗砂土顆粒堆積試驗可測得實際試驗砂土顆粒的自然休止角，室內砂土自然休止角測試方法主要采用沖擊法，即將試樣放置在漏斗中然后使其自由下落在方格紙上形成小錐體，通過測量其錐底面直徑和錐頂高度以計算出自然休止角。

為對比驗證本發(fā)明，室內砂土顆粒堆積試驗中采用4種散粒為驗證試驗對象，材料E為純圓玻璃顆粒，材料F、G、H為中-粗石英砂顆粒。將4種試樣在光學顯微鏡下拍攝的照片用ImageJ軟件進行黑白二值化處理，經分析后提取出顆粒單元的幾何參數，具體形狀參數如表1所示。

表1 各種顆?；境叽鐓?/p>

（2）根據步驟（1）中得到的實際試驗砂土顆粒的幾何參數為基礎，建立一種離散元組合顆粒數值模型，而球形度為離散元組合顆粒外形參數，根據體積等效及質量等效原則以離散元組合顆粒代替實際試驗砂土顆粒，以模擬不同球形度的試驗砂土，顆粒間的接觸選擇線性接觸模型。

（2.1）所述的離散元組合顆粒包括各種球形度的Clump顆粒，所述的Clump顆粒由單元球A、單元球B、單元球C和單元球D四個單元球構成，所述的單元球A、單元球B、單元球C的半徑相等，所述的單元球A、單元球B、單元球C以圓周分布方式排列于單元球D的四周，且單元球A、單元球B、單元球C均與單元球D相切或相交，所述Clump顆粒的球形度n的取值為0.33≤n≤1；其中所述的球形度n（亦稱球度）的定義為顆粒的形狀與球體相似的程度，且n=最大內切圓半徑/最小外接圓半徑[如圖1]。

（2.2）在生成Clump顆粒時，單元球A、單元球B、單元球C半徑相等且不變，改變單元球D的半徑大小可獲得不同的球形度的Clump顆粒，根據室內實際試驗砂土顆粒堆積試驗的參數可選取Clump顆粒的球形度為0.33，0.4，0.5，0.58，0.70，0.77，0.9或1.0中的1個、2個或多個，分析顆粒球形度對散粒試樣的宏細觀物理力學特性的影響。所述的Clump顆粒的球形度n與各個單元球的半徑取值關系可參見表2，設與Clump顆粒具有相等體積的圓球顆粒半徑為1，對應的Clump顆粒的各個單元球半徑可按表2選取。

表2 Clump顆粒的各個單元球的相對球半徑

（3）用離散元PFC3D方法分兩段建立漏斗模型。該漏斗的上部采用兩端開口的圓筒墻生成的圓筒結構，圓筒直徑為50mm，高度為55mm；漏斗下部采用上端開口下端閉口的圓錐墻生成的圓錐結構，圓錐上端直徑為50mm，下端直徑為20mm，圓錐高度為15mm。在距漏斗下端50mm處生成底部墻面模擬地面，并設置相應的細觀參數；在漏斗下端開口生成的閉合墻體可防止漏斗中的Clump顆粒下落。

（4）根據步驟（1）的參數首先生成對應球形度的離散元組合顆粒模板，然后根據體積等效與質量等效原則在漏斗上方生成與實際試驗砂土顆粒體積一致的各種球形度的離散元組合顆粒，設置離散元組合顆粒細觀參數并施加重力加速度為g＝-9.81m/ s²，然后開始循環(huán)計算，此時離散元組合顆粒在重力作用下逐漸沉積在漏斗內部，設置平衡條件并計算至初始平衡狀態(tài)，在漏斗內部布置測量球測出離散元組合顆粒的平均接觸數及孔隙率；

（4.1）所述的離散元組合顆粒為組合Clump顆?；騿我籆lump顆粒，所述的組合Clump顆粒為包括多種不同球形度的Clump顆粒組合成的有多種外形的顆粒組合，所述的單一Clump顆粒為由1種球形度的Clump顆粒組合成的只有一種形狀的顆粒組合。

（5）刪除步驟（4）中漏斗底部的閉合墻體使漏斗中的離散元組合顆粒逐漸下落直至漏斗下方的底部墻面上并形成穩(wěn)定的顆粒堆積體[如圖2(a)]。

（6）通過在步驟（5）中的顆粒堆積體中布置測量球[如圖2(b)]可測得測量球范圍內離散元組合顆粒的平均接觸數及孔隙率，編寫fish函數輸出測量顆粒堆積體表面顆粒坐標以計算自然休止角[如圖2(c)]。

（6.1）在完成步驟（5）中的堆積過程后將離散元組合顆粒球心坐標導出，對非錐頂和非錐腳區(qū)域用最小二乘法對顆粒堆積體表面曲線進行線性擬合，根據擬合直線的斜率求出砂土堆積體的自然休止角；

（6.2）顆粒堆積體通常呈現不對稱狀態(tài)，為消除方向選擇對休止角測定結果的影響，從顆粒堆積體的x、y坐標軸正負四個方向選取表面坐標，并將擬合結果取均值作為顆粒堆積體的自然休止角。

（7）根據步驟（6）中布置在顆粒堆積體中的測量球，可實時獲取離散元組合顆粒堆積體的孔隙率，通過改變顆粒的細觀參數值，可實現求取顆粒堆積體的最大及最小孔隙率，亦可分別研究球形度n、顆粒摩擦系數、孔隙率、接觸數、自然休止角等各參數的相互關系，見圖4、圖5、圖6，可以發(fā)現堆積體顆粒球形度相同時，自然休止角和孔隙率隨摩擦系數的增大而增加，平均接觸數隨摩擦系數的增大而減小；相同摩擦系數條件下，自然休止角隨著顆粒球形度的減小而增加，平均接觸數則與孔隙率呈現負相關關系。

表3為室內實際試驗砂土顆粒堆積試驗與離散元組合顆粒的堆積模擬試驗的結果對比，根據表3結果發(fā)現，單一Clump顆粒（即單一球形度Clump顆粒，指的是球形度大小相同的clump顆粒）的最大、最小孔隙比均與室內試驗結果接近。值得注意的是，數值模擬中最小孔隙比為理想光滑狀態(tài)，顆粒與顆粒之間不存在摩擦，實際顆粒之間必然存在摩擦力，因而導致模擬的最小孔隙比均小于實際值，E組試樣為光滑玻璃球顆粒，其摩擦力接近理想狀態(tài)，故與模擬值最接近。球形度較大的E和F試樣，模擬試驗得到的最大孔隙比均略大于實際驗證試驗，球形度較小的G和H兩組，模擬試驗的最大孔隙比均略小于實際驗證試驗得到的值，此現象說明如采用單一Clump顆粒無法完全模擬實際顆粒。根據圖9中對土樣G的顆粒形狀統(tǒng)計結果，將不同球形度的Clump顆粒進行組合構成組合Clump顆粒，組合Clump顆粒的最小孔隙比相對單一Clump顆粒有提高，最大孔隙比與實際最大孔隙比吻合，說明根據顆粒球形度統(tǒng)計方法能有效模擬實際砂土顆粒，即組合Clump顆粒更能有效模擬實際砂土顆粒。

表3室內試驗與離散元組合顆粒數值模擬結果對比

顆粒形狀及摩擦系數對堆積體自然休止角的影響：由圖4可發(fā)現自然休止角隨著顆粒球形度的減小而增加，其增加速率隨球形度的減小而減??；自然休止角與摩擦系數存在正相關關系。

顆粒形狀及摩擦系數對堆積體孔隙率的影響：由圖5發(fā)現堆積體顆粒球形度相同時，孔隙率隨摩擦系數的增大而增加，其增加速率與摩擦系數呈現負相關關系；堆積體在相同摩擦系數條件下，顆粒球形度的增加對孔隙率的影響表現出先減小后增加的趨勢。當球形度小于0.5時，孔隙率隨球形度增加而減小的主要原因為Clump顆粒單元球顆粒體積增加；球形度在0.5-0.7時，孔隙率的變化不大，此時單元球顆粒體積增加導致的孔隙率減小剛好與拱架結構孔隙的增加導致的孔隙率增加效應抵消；當球形度大于0.7后，孔隙率增加的主要原因為顆粒形成的拱架結構中孔隙增加；當球形度為0.9和1.0時，Clump顆粒形狀已經類似球顆粒，下落后顆粒四處分散且難以形成穩(wěn)定的堆積體，從而導致孔隙率的突變。

顆粒形狀及摩擦系數對堆積體平均接觸數的影響：在圖6中，堆積體顆粒球形度相同時，平均接觸數隨摩擦系數的增大而減小，其減小速率與摩擦系數呈現負相關關系；堆積體在相同摩擦系數條件下，顆粒球形度的增加對平均接觸數的影響表現出先增加后減少的趨勢。對照圖5和圖6可以發(fā)現，顆粒的平均接觸數與試樣的密實程度有關，試樣越松散顆粒與顆粒之間接觸的數量越少，反之則平均接觸數越多。

組合Clump顆粒對試驗結果的影響：圖7為組合Clump顆粒與單一Clump顆粒模擬驗證土樣G的對比結果（目的在于比較多種外形組合的砂土顆粒比單一形狀組合的顆粒更符合實際砂土情況）。從圖7中發(fā)現，組合Clump顆粒在相同試驗條件下孔隙率、平均接觸數和自然休止角均大于單一Clump顆粒，其中較為明顯的是組合Clump顆粒平均接觸數大于單一Clump顆粒的10%-20%。根據圖9室內試驗土樣G球形度統(tǒng)計結果，球形度小于0.7的顆粒頻數相對球形度大于0.7的顆粒較多，且小球形度對接觸數的敏感度大于大球形度；圖5已顯示球形度為0.7時孔隙率為最小值，從而出現組合多種球形度顆粒后孔隙率增大；組合Clump顆粒由于顆粒形狀上具備多樣性，球形度較小的顆粒加強了顆粒之間的嵌固效應，從而增強了顆粒的咬合摩擦力，因而組合Clump顆粒的自然休止角增加。

從圖8可以得到，采用組合Clump顆粒模擬土樣G的自然休止角α與摩擦系數的關系可擬合為：

α = -34.492f ² + 60.648f + 9.4278 （1）

式中α為自然休止角，f為摩擦系數，擬合相關系數R²為0.9788。

根據室內直剪試驗結果土樣G的穩(wěn)定狀態(tài)內摩擦角θ為31°，砂土的摩擦系數通常為f=tan(θ)，此時土樣G的室內摩擦系數為0.6，代入擬合公式得模擬的自然休止角為α=33.41°。文獻(范智杰, 屈建軍, 周煥. 沙土內摩擦角與粒徑、含水率及天然坡角的關系[J]. 中國沙漠, 2015, 35(2):301-305.)分別對7種粒徑的砂土進行了堆積試驗和直剪試驗，發(fā)現砂土的自然休止角均略大于內摩擦角3°左右。該結果與本發(fā)明采用組合Clump顆粒進行數值模擬得出的結果相一致，說明采用組合Clump顆粒模擬方法可直接采用試樣實測摩擦系數，驗證了離散元組合顆粒模擬方法的可行性。