摘要：利用隨機集成模型模擬單一粒徑橢球體顆粒的隨機堆積過程，并分析了邊界條件與顆粒幾何形狀對堆積結構的影響。研究表明：相對邊界尺度（顆粒當量半徑與立方體容器邊長之比）對散體結構的密度、配位數影響顯著；三維尺寸越接近，顆粒堆積越密實；散體顆粒的平均配位數與其對應的密實度隨著縱橫比的變化趨勢較為一致，均隨著橢球體縱橫比的增大而減小。該程序為模擬多相顆粒體（如道砟、瀝青混合料等）的堆積規律、級配等奠定了基礎。
　　關鍵詞：顆粒堆積；空隙率（密度）；配位數
　　
　　0引言
　　目前大多數離散元分析模型是多邊形、圓形，前者常用在巖體模型中，后者常用在土顆粒模型中。用橢圓形代替圓形來模擬顆粒介質似乎已經成為一種提高模擬效果的發展趨勢，已有眾多學者對二維橢圓形離散元模型作過研究[1-2]，并在土顆粒的研究中有所應用[3]，而且Rothenburg等對橢圓形與圓形顆粒模型進行了對比，可知用橢圓形表示顆粒在散體的特性方面要優于用圓形表示的顆粒[4]。然而，二維顆粒堆積在現實中是不存在的，它與連續介質力學中的平面應變和應力的情況是不一樣的。瀝青混合料、道砟等散體的一些重要指標，如空隙率和配位數（單位顆粒與周圍顆粒接觸的數量），很難用二維顆粒加以清楚的描述。
　　對于三維顆粒模型來說，圓球體顆粒堆積的研究較為廣泛而深入，其結論大都已取得學術界的共識，工程實踐中也常將非球體用當量球來簡化。但近年來的研究表明，非球體在堆積特性上與球體有很大的不同，顆粒形狀的微小變化會使堆積密度發生顯著的改變。然而由于模型及計算的復雜性，人們對非圓球體堆積的研究相對較少。而橢球體顆粒由于其三軸在長度上彼此獨立，可以調整，在某種程度上更接近于現實顆粒。本文將建立橢球體顆粒的隨機堆積空間分布模型，在VC++平臺上開發出相應的計算程序，并將其用于分析表征堆積結構特性的重要參數：空隙率、配位數，為今后進一步研究復雜形狀顆粒的微觀結構打下基礎。
　　1模型及程序設計
　　1.1顆粒的生成
　　顆粒空間分布模型的實質是顆粒隨機堆積集成模擬，即模擬顆粒逐個在容器內的堆積，并最終形成穩定的散體結構。堆積開始的第一步就是顆粒的隨機生成，由于本文研究的顆粒粒徑是一定的，所以其隨機性表現為顆粒主軸的初始方向向量和質心在上的初始坐標。顆粒沿方向落入已存在的散體界面上并有第一接觸點，在其下落I階段中，顆粒在直角坐標系平面上的坐標是不變的，豎向坐標是變化的，因此的初始值可設置為容器的高度。
　　MonteCarlo方法亦稱統計模擬方法(StatisticalSimulationMethod)，是一種利用隨機數進行數值模擬的方法。它是以概率統計理論為基礎的一種方法。由于蒙特卡羅方法能夠比較逼真地描述事物的特點及物理試驗過程，解決一些數值方法難以解決的問題，因而該方法的應用領域日趨廣泛。本文主要利用蒙特卡羅法，在立方體容器內，借助隨機數產生粒子中心的平面坐標和主軸的方向余弦，使其分別隨機分布在(為立方體容器的邊長)和內。
　　1.2顆粒的接觸判斷
　　在一個由眾多顆粒組成的體系中，直接判別顆粒是否接觸需要耗費大量的計算時間，為提高計算效率，本文將不直接精確判別任意兩個顆粒間是否存在接觸，而是分兩個步驟判別顆粒間的接觸是否存在。對當前所研究顆粒首先采用橢球體的最大外接球初步的判斷其潛在的鄰居個數，然后采用幾何勢能法[5]精確的確定該顆粒與每個鄰居是否接觸。雖然在確定鄰
　　居數目時也要耗費一定的計算時間，但是比逐個準確判別顆粒間是否存在接觸要節約很多時
　　間。
　　(1)初步接觸判斷
　　對于初步接觸判斷滿足下面的條件：
　　  (1)
　　式(1)中，，為顆粒和質心之間的距離，，分別為顆粒和的中心坐標，為顆粒的長軸半徑，為用戶輸入的常數，調節距離，當顆粒的累積位移超過時，才進行全場搜索，更新其初步接觸鏈。當顆粒通過初步接觸判斷，方可進入細步接觸判斷。
　　(2)細步接觸判斷
　　如圖1所示，設兩相鄰顆粒和在全局坐標系下的曲面方程分別為：
　　  (2)
　　  (3)
　　式(2)與式(3)中,、(…..)是常數。為了判斷顆粒是否與顆粒接觸，則問題轉換
　　為函數：
　　  (4)
　　在式(3)約束下求極值的問題，本文采用拉格朗日條件極值，引入變量，構造如下函數：
　　
　　 。 (5)
　　對式(5)求微分并整理得： 
　　  (6)
　　通過式(6)，在顆粒上可求得離顆粒的最近點的坐標。如果，則兩個橢球顆粒接觸存在，同理可求得點的坐標，接觸點坐標即為點與點連線的中點，如果，則兩個橢球顆粒沒有接觸。
　　
　　圖1基于幾何勢能法的接觸模型
　　1.3程序設計
　　瀝青混合料或道砟等散體在形成過程中通過振動、搗實等外荷載作用才得以密實，因此顆粒在空間的隨機位置應同時滿足勢能最小原理。在顆粒隨機堆積模擬過程中，將顆粒堆積體的構造過程看作是一個順序過程，即顆粒一旦定位，則認為已處于穩定狀態，不再受其它任何顆粒的作用或擾動而產生移動，與邊界條件共同形成待生成顆粒的約束。
　　除少數停留在容器底部的顆粒外，大多數顆粒只有形成三點支撐結構才能穩定。為問題敘述方便，將這三個支撐點按其先后形成順序分別叫做第一、第二和第三落點。相應地，根據這三個接觸點的形成，將當前顆粒開始下落直到三點支撐結構形成的全過程分為三個階段：下落I、II、III階段。本文以MicrosoftVisualStudio6.0為平臺進行三維橢球體顆粒隨機堆積的程序開發，并實現堆積結構的可視化。
　　2.實例模擬及結果分析
　　2.1模擬條件
　　為便于研究分析，橢球體的長、中、短軸半徑滿足：。本文對單一粒徑橢球體縱橫比值（長短軸比值）分別為1.0，1.3，1.5，1.7，2.0五個樣本的隨機堆積過程進行模擬，圖2即為程序隨機集成的五組不同縱橫比的可視化散體結構，同時對堆積體的空隙率和配位數隨相對邊界尺度的變化情況進行研究分析，相對邊界尺度定義為：橢球體的當量半徑與立方體容器邊長的比值，橢球體的當量半徑定義為：與橢球體等體積時所對應的圓球體半徑。為了消除邊界效應，利用回歸法得到零邊界條件下堆積體的空隙率和配位數。
　　
　　a)1.0b)1.3c)1.5d)1.7e)2.0
　　圖2五個樣本的散體結構
　　2.2結果分析
　　2.2.1空隙率
　　空隙率定義為顆粒間空隙體積與堆積體總體積的比值：，它是表征顆粒堆積特性的基本參數，其大小直接反映堆積的密實度，堆積密實度可表示為。圖3為散體的不同縱橫比對應的空隙率隨相對邊界尺度的變化情況，由圖可知，在五種情況下，散體的空隙率均隨相對邊界尺度的減小而減小[6]；為避免邊界效應，將五條擬合直線分別線性外推至零邊界尺度，得到無邊界影響下，五個樣本對應的散體空隙率分別為0.43180，0.46157，0.47847，0.49060，0.51291，就本文所研究的條件而言，空隙率隨著縱橫比的增大而呈增大趨勢，與Ng研究成果[7]基本一致，由于受力條件及堆積方式不同，所以大小范圍有一點差距，但是總體變化趨勢是一致的。
　　
　　a)1.0b)1.3
　　
　　c)1.5d)1.7e)2.0
　　圖3橢球粒子形成的堆積體在五種縱橫比下的空隙率隨相對邊界尺度的變化
　　2.2.2配位數
　　配位數定義為與研究顆粒接觸的相鄰顆粒的個數，反映散體材料內部力的傳遞、結構穩定度、強度等特性。配位數的大小和分布與接觸定義的極限距離（簡稱極距）有很大的關系，極距越大，配位數越大[8]，本文認為凈距小于0.1倍橢球體的當量半徑相互接觸[9]。散體結構內部平均配位數計算域的選取方法通常有兩種[10]，一種是wholecontainer，即容器框架內所有的粒子，另一種是smallbox，即容器中心部分一半邊界長度范圍的粒子。本文采用前者，對容器內所有的粒子進行配位數的計算，并取其平均值。圖4為橢球粒子在五個樣本下的平均配位數隨相對邊界尺度的變化情況。由圖可知，五組不同縱橫比的散體結構對應的配位數均隨著相對邊界尺度的減小而增大，為了避免邊界因素的影響，將五條擬合直線分別線性外推至零邊界尺度，得到無邊界效應下，五組不同縱橫比的配位數分別是5.98374，5.74890，5.36023，5.57413，5.01903。為了便于分析，繪制如圖5，由圖可知，就本文所研究的條件而言，配位數隨著縱橫比的增大而呈減小趨勢。
　　
　　a)1.0b)1.3
　　
　　c)1.5d)1.7e)2.0
　　圖4橢球粒子在五個樣本下的平均配位數隨相對邊界尺度的變化
　　為了便于分析比較散體密實度與平均配位數的關系，將五個樣本在零邊界效應下的密實度隨縱橫比值的變化情況也繪制在圖5中。由圖可知，散體粒子的平均配位數與其對應的密實度隨著縱橫比的變化趨勢基本一致，均隨著橢球體縱橫比的增加而減小，即散體粒子堆積的越密實，其對應的平均配位數也越大，這與已有文獻的研究成果相吻合。歐陽鴻武等介紹了單一尺寸球形粉末顆粒堆積密度與配位數和排列方式的關系，即顆粒的堆積密度隨配位數的增大而提高[6]；Bernal等對鋼球進行了隨機松散和密實堆積，得到隨機密實堆積的配位數大于隨機松散堆積的配位數，并指出隨機松散堆積粒子之間有更大的空隙，而且粒子之間的配位數是不對稱的，比如在某一側有其它粒子與之接觸，在另一側則是空隙[11]。
　　
　　圖5橢球粒子形成的散體密實度及其平均配位數隨縱橫比的變化
　　3結論
　　利用三維橢球體顆粒的隨機集成模型模擬顆粒的隨機堆積過程，以能夠表征散體顆粒微觀結構特征的重要參數空隙率和配位數為研究指標，就顆粒所在的容器邊界條件及其幾何形狀這兩個方面對單一粒徑顆粒堆積結構的影響進行討論分析，得到如下結論：散體結構的空隙率和配位數隨著顆粒相對容器邊界尺度的減小而分別減小和增大；顆粒的三維尺寸越接近，散體結構的空隙率越小，堆積越密實；散體粒子堆積的密實度與其平均配位數有密切的聯系，兩者均隨著縱橫比的增大而減小，堆積越密實，其對應的平均配位數越大。
　　本文通過三維數值模擬得到的以上結果與前人的研究成果較吻合，該成果是研究雙粒徑混合橢球體以及不同大小分布的多粒徑混合橢球體堆積規律的基礎，同時從圖3可以清晰看出顆粒骨架結構（其內部結構可操作鼠標觀察）的實際效果，對今后瀝青混合料、道砟等散體的可視化級配設計有較強的指導意義。