顆粒材料數值樣本的坐標排序生成技術

楚錫華

（1. 武漢大學 工程力學系，武漢 430072；2. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072)

1 引 言

顆粒材料的宏觀性質與其微觀結構（顆粒形狀、粒度分布以及配位數等）密切相關，特別是顆粒材料的孔隙度，或者說結構密度關系到顆粒材料的諸多性質，如力學響應、滲透性、溶解質和熱量的擴散性能等。因此，Herrmann[1]認為，結構密度是描述顆粒系統行為的最重要參數，它實際上控制了顆粒材料的力學行為。孔隙度是顆粒集合中孔隙體積與顆粒集合總體積的比值，在二維情況下，孔隙度是孔隙面積與顆粒集合總面積的比值。它是顆粒材料的宏觀參數，表征了顆粒材料的相對密度。與其相應的微觀參數是顆粒材料的平均配位數（average coordinate number）。一個顆粒的配位數是指與之相接觸的顆粒的個數（主要是針對圓形顆粒或橢圓形顆粒定義的參數），平均配位數是指顆粒集合內所有顆粒配位數的算術平均值。孔隙度與顆粒材料的顆粒形狀、粒度分布及排列方式等有密切的關系。顯然顆粒材料的孔隙度不能任意取值，對特定粒度分布的顆粒材料，一方面在顆粒相互不重合的情況下孔隙度有一個下限值，如對單粒度圓形顆粒（二維情況）組成的顆粒集合，其孔隙度最小值為0.093 1，單粒度球形顆粒（三維）組成的顆粒集合，其最小孔隙度為0.259 5；另一方面，雖然顆粒材料的孔隙度理論上可趨近于 1，但顆粒稀疏到一定程度后，顆粒之間將不接觸，此時顆粒集合的平均配位數將趨于0，顆粒集合失去承載能力。

此外，以離散顆粒模型研究巖土顆粒材料時，由于不同類型的巖土材料在孔隙度以及粒度分布上存在一定的差別，為了更真實地模擬不同巖土材料的力學行為，在數值模擬時有必要研究顆粒材料的隨機生成技術，以使生成的顆粒材料數值樣本具有相應的孔隙度及粒度分布。因此，顆粒材料數值樣本生成及顆粒排列技術日益受到重視[2－8]，研究者提出了多種樣本生成方案，然而至今仍無一條公認的優化途徑[7]。文獻[2]指出，生成顆粒材料的困難在于隨機生成時布點的分散性和生成的顆粒材料的真實性。魏群[9]基于RSA模型給出了隨機布點及挑選滿足要求的顆粒的方法。本文在此基礎上改進了其隨機布點方式，即先通過對隨機生成的x坐標序列或y坐標序列進行排序，然后再形成坐標對進行布點，從而使在給定的區域內能夠生成更多符合要求的顆粒。

2 顆粒材料數值樣本生成的RSA模型

2.1 經典RSA模型

RSA模型也稱為 SSI模型（simple sequential inhibition model）。其主要思想如下：首先按粒度分布要求生成備選粒徑序列{ri}，然后在給定二維或者三維區域內按照隨機方式逐一給定r1,r2......的中心位置設當前顆粒為i，若顆粒i與已存在的顆粒重合，則將i舍棄，繼續用i+1試探，該過程直至給定區域填充滿，亦即任何顆粒無法放入該區域為止。為簡便起見，以二維為例，在所考慮區域Ω內隨機布點i（ xi,yi），且滿足如下條件 xi?[Xmin,Xmax]yi?[Ymin, Ymax]，其中Xmin、Xmax、Ymin、Ymax為區域Ω最小外接矩形的角點坐標值。xi，yi通常以隨機數列的方式給出，即

2.2 RSA模型的若干改進

經典RSA模型雖然簡單，但只能生成非常稀疏的顆粒集合體，一旦要求生成的顆粒體的孔隙度較小時，會造成相當多的顆粒因重疊而被舍棄，使得顆粒集合的粒度分布與目標粒度分布往往有較大的差異（主要是由于在生成非單粒度顆粒群的過程中，較小的顆粒總有較多的機會在給定區域內找到位置）。為此，在應用時需對RSA模型進行改進。

（1）改進方案1

該改進方案主要針對當前顆粒i與已生成的顆粒存在一定重疊時，并不直接舍棄顆粒i，而是根據重疊量作如下處理：

設與其重疊的顆粒為a，其半徑為ra，坐標為重疊量為Δ，則存在如下關系

式中：d為兩顆粒中心之間的距離，顆粒i的坐標按如下方式修正

應用新的坐標判斷顆粒i是否與已存在的顆粒重疊，若不重疊，則保留該顆粒，否則，繼續修正，若修正k次后，仍與已存在的顆粒重疊，則舍棄該顆粒。

（2）改進方案2

該改進模型為：若當前顆粒i與已生成顆粒重疊時，僅僅舍棄坐標而保留ri，并將以ri為半徑，以坐標為圓心的顆粒視為當前顆粒。從理論上看，這一方案能較好地維持粒度分布的要求，但要求備選的坐標序列要遠遠大于目標顆粒數目。若將能夠得以保留的坐標稱為存活坐標，將存活坐標的數目與備選坐標的數目之比稱為坐標存活率，則上述改進方案的坐標存活率很低。對同一粒度序列，不難理解對備選坐標數目相同的坐標序列，坐標存活率的提高在某種程度上可使生成密實顆粒集合的可能性更大。這是由于當坐標存活率較低時，往往會發生粒度序列仍未循環結束，而備選坐標序列已用完，但此時生成的顆粒數目仍未滿足要求，因而在商用顆粒流軟件 PFC2D中，若保持粒度大小不變的情況下常會發生無法生成目標顆粒數目的情況。

（3）改進方案3

該方案在方案1與方案2的基礎上，若顆粒重疊時，并不直接舍棄當前坐標，而是根據重疊量應用方案1調整當前顆粒的坐標，若調整k次后，仍與已存在的顆粒重疊，則保留半徑ri，重新選擇下一個隨機坐標，該方案得到的樣本一般仍比較松散。

（4）改進方案4

增加生成顆粒數目的一種最直接的方式是增大粒度序列與坐標序列。但在這里筆者將尋求另外一種方式，即保持現有粒度序列與坐標序列大小不變的情況下，如何提高生成顆粒的數量。注意到經典RSA模型及改進方案2中對坐標序列及種子序列的編號是隨機數列進行的，為了簡單起見，假定粒度序列為等粒度序列，由于對顆粒的挑選是按照編號進行的，則可能會出現圖 1(a)所示的情況，在這種情況下，顆粒1首先生成，那么顆粒2（坐標2）和顆粒3（坐標3）都將因為和顆粒1重疊而被舍棄。但如果將圖1（a）的3個坐標序列排序后按照圖1（b）進行編號，那么只有顆粒2因重疊被舍棄，從而提高坐標的存活率，因此，首先對進行排序，然后再形成坐標序列（ xi, yi），結合經典的RSA模型可稱為改進方案4a，結合改進方案1稱為改進方案4b，結合改進方案2稱為改進方案4c，結合改進方案3成為4d。下面將以一個簡單的算例來說明經過坐標排序后的方案能夠生成更多的顆粒。

圖1 顆粒隨機生成示意圖Fig.1 Diagram of random packing

3 數值算例

在寬為20 cm、高為40 cm的矩形區域內生成半徑為0.3、0.4、0.5 cm的顆粒集合（其顆粒數量之比為 4:3:3）。其中，隨機數列采用混同余法，其遞推公式可表示為

式中：{ri}為待生成的隨機數列；通常選取D為正奇數； M=2k； C=4q+ 1；X0為任意非負整數，通常隨機數列的種子，在本算例中 D= 13 849，C=2 053，M= 65 536。的種子為20 001，的種子為60 001。表1給出了采用不同方案時，生成的顆粒總數量隨粒度序列的變化，為了簡便起見，這里取粒度序列的大小與坐標序列相等。表1中N表示粒度序列（坐標序列）的大小，從表可以看出，各個方案隨粒度序列（坐標序列）的增大，生成的顆粒數目在遞增。對比各個方案可看出改進后的方案均優于經典RSA模型，包含坐標修正的方案優于不包含修正的方案，即方案1、方案3優于方案2，方案4b與方案4d優于4a及4c。并且經過坐標排序的方案均優于與其對應的未經排序的方案，即方案4a優于經典RSA，方案4b優于方案1，方案4c優于方案2，方案4d優于方案3。

表1 各個RSA模型生成顆粒數量的對比Table 1 Comparison of number of grains generated by different RSA methods

下面將著重探討當N =5 000時生成顆粒數在1 000以上的3種方案，即方案1、方案3、方案4b及4d。圖2(a)～(d)給出了其對應生成的顆粒集合，統計可知，其生成的半徑為0.3、0.4、0.5 cm的顆粒的數量之比分別為（547：271：190）；（505：306：264）；（521：314：268）；（453：339：327），對比可看出，方案4d生成的顆粒級配最接近要求的級配（4：3：3）。下面從生成的顆粒集合的接觸分布探討顆粒集合的各向異性情況，圖3給出了方案1、方案3、方案4b及方案4d生成的顆粒集合的接觸法向分布。

可以看到未經坐標排序的方案生成的顆粒集合較接近各向同性，而以x坐標排序的方案生成的顆粒集合在x方向上的接觸較為密集。對圓形顆粒集，在理論上接觸分布應是關于原點對稱的，但對稀疏顆粒集合接觸統計時誤差的影響比較顯著，因此，圖 3(a)，3(b)在一定程度上喪失了對稱性。此外，需指出坐標排序方案雖然能生成較多的顆粒，但通常其顆粒集合的密實程度仍不能滿足顆粒材料二維離散單元數值模擬的樣本要求，如方案4 d雖生成了1 119個顆粒，但此時的孔隙度為0.305，顆粒的平均配位數（接觸點個數）為1.89。通常而言，RSA模型可通過顆粒膨脹并結合顆粒樣本生成的動力方案[8]，從而使顆粒材料樣本達到離散單元數值模擬的要求。

圖2 不同方案生成顆粒集合Fig.2 Granular assemblies generated by different methods

圖3 不同方案生成的顆粒集合的接觸法向分布圖Fig.3 Distributions of contact normals of granular assemblies generated by different methods

4 結論與評述

（1）本文在顆粒材料數值樣本生成的經典RSA模型基礎上，討論了幾種修正方案，并建議了一種基于坐標排序的樣本生成技術。數值算例表明，坐標序列經過排序后，再結合相應的RSA模型能夠在一定程度上提高坐標序列的存活率，從而使生成的顆粒集合趨于密實。坐標排序屬于坐標序列的規則化，這與 Sobolev[10]所建議的技術在一定程度上相類似。需要指出，顆粒的規則排列及隨機排列技術，特別是如何隨機排列使生成的顆粒材料趨于更密實，至今仍是一個活躍的課題[11－12]。本文僅以RSA模型為基礎從幾何構造的角度進行了初步的探討。文中所給例題雖為二維算例，但該方法不難推廣至三維情況。

（2）目前顆粒材料數值樣本的生成方法還處在百家爭鳴的階段，總體而言，可分為動力方法與構造方法[4－8]。構造方法在生成樣本過程中不涉及顆粒的物性，因此，所生成的初始數值樣本可用于模擬不同物性的顆粒材料；動力方法在生成樣本的過程中需用到顆粒間的相互物理作用，所生成的顆粒樣本與顆粒間的相互作用密切相關。此外，需特別值得關注的是，基于網格的顆粒樣本幾何構造算法[7－13]因能較好地應用有限單元法的網格剖分技術而具有一定的優勢。

（3）顆粒材料數值樣本生成問題本質上屬于顆粒排列堆積問題。顆粒排列特別是圓形（球形）顆粒的排列問題在學術（包括數學，物理、材料，巖土、化工等領域）及工業應用中都具有重要的意義，因而成為經久不衰的研究主題[1－17]。經典顆粒排列問題主要關心在不同形狀的容器內如何用等尺寸圓盤得到最致密排列。此外比較關注生成結構的幾何性質，如顆粒尺寸分布與指定尺寸分布的一致性、接觸的各向同性特性、配位數等。對顆粒材料數值樣本的要求及各類生成技術的詳細比較筆者將在另文表述。

[1]HERRMANN H J. Granular matter[J]. Physical A., 2003,313: 188－210.

[2]LAURENT T, PIERRE P, AURELE P. Generation of granular media[J]. Transport in Porous Media, 1997, 26:99－107.

[3]WANG C Y, VEI CHUNG L. A packing generation scheme for the granular assemblies with planar elliptical particles[J]. Int. J. Numer. Anal. Methods Geomech.,1997, 21: 327－358.

[4]CHENG Y F, GUO S J, LAI H Y. Dynamic simulation of random packing of spherical particles[J]. Powder Technology, 2000, 107: 123－130.

[5]FENG Y T, HAN K, OWEN D R J. Filling domains with disks: An advancing front approach[J]. Int. J. Numer.Meth. Engng., 2003, 56: 699－713.

[6]JIANG M J, KONRAD J M, LEROUEI S. An efficient technique for generating homogenous specimens for DEM studies[J]. Computers and Geotechnics, 2003,30(7): 579－597.

[7]CUI L, SULLIVAN O. Analysis of a triangulation based approach for specimen generation for discrete element simulations[J]. Granular Matter, 2003, 5: 135－145.

[8]BAGI K. An algorithm to generate random dense arrangements for discrete element simulations of granular assemblies[J]. Granular Matter, 2005, 7(1): 31－43.

[9]魏群. 巖土工程中散體單元的基本原理、數值方法及實驗研究[D]. 北京: 清華大學, 1990.

[10]SOBOLEV K, AMIRJANOV A. A simulation model of the dense packing of particulate materials[J]. Advanced Powder Technology, 2004, 15(3): 365－376.

[11]TORQUATO S, TRUSKETT T M, DEBENEDETTI P G.Is random close packing of spheres well defined?[J].Physical Review Letters, 2000, 84(10): 2064－2067.

[12]ZAMPONI F. Packings close and loose[J]. Nature, 2008,53: 606－607.

[13]JERIER J F, RICHEFEU V, IMBAULT D, et al. Packing spherical discrete elements for large scale simulations[J].Computer Methods in Applied Mechanics and Engineering, 2010, 199: 1668－1676.

[14]SCOTT G D. Packing of spheres[J]. Nature, 1960, 188:908－909.

[15]ZASKI A M. An efficient method for packing polygonal domains with disks for 2D discrete element simulation[J].Computers and Geotechnics, 2009, 36(4): 568－576.

[16]HERN C S, LANGER S A, LIU A J, et al. Random packing of frictionless particles[J]. Physical Review Letters, 2002, 88(7): 1－4.

[17]BENABBOU A, BOROUCHAKI H, LAUG P, et al.Numerical modeling of nanostructured materials[J].Finite Elements in Analysis and Design, 2010, 46(1－2):165－180.