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