楚錫華

（1. 武漢大學(xué) 工程力學(xué)系，武漢 430072；2. 武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

1 引 言

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

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

2 顆粒材料數(shù)值樣本生成的RSA模型

2.1 經(jīng)典RSA模型

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

2.2 RSA模型的若干改進(jìn)

經(jīng)典RSA模型雖然簡(jiǎn)單，但只能生成非常稀疏的顆粒集合體，一旦要求生成的顆粒體的孔隙度較小時(shí)，會(huì)造成相當(dāng)多的顆粒因重疊而被舍棄，使得顆粒集合的粒度分布與目標(biāo)粒度分布往往有較大的差異（主要是由于在生成非單粒度顆粒群的過(guò)程中，較小的顆?？傆休^多的機(jī)會(huì)在給定區(qū)域內(nèi)找到位置）。為此，在應(yīng)用時(shí)需對(duì)RSA模型進(jìn)行改進(jìn)。

（1）改進(jìn)方案1

該改進(jìn)方案主要針對(duì)當(dāng)前顆粒i與已生成的顆粒存在一定重疊時(shí)，并不直接舍棄顆粒i，而是根據(jù)重疊量作如下處理：

設(shè)與其重疊的顆粒為a，其半徑為ra，坐標(biāo)為重疊量為Δ，則存在如下關(guān)系

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

應(yīng)用新的坐標(biāo)判斷顆粒i是否與已存在的顆粒重疊，若不重疊，則保留該顆粒，否則，繼續(xù)修正，若修正k次后，仍與已存在的顆粒重疊，則舍棄該顆粒。

（2）改進(jìn)方案2

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

（3）改進(jìn)方案3

該方案在方案1與方案2的基礎(chǔ)上，若顆粒重疊時(shí)，并不直接舍棄當(dāng)前坐標(biāo)，而是根據(jù)重疊量應(yīng)用方案1調(diào)整當(dāng)前顆粒的坐標(biāo)，若調(diào)整k次后，仍與已存在的顆粒重疊，則保留半徑ri，重新選擇下一個(gè)隨機(jī)坐標(biāo)，該方案得到的樣本一般仍比較松散。

（4）改進(jìn)方案4

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

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

3 數(shù)值算例

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

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

表1 各個(gè)RSA模型生成顆粒數(shù)量的對(duì)比Table 1 Comparison of number of grains generated by different RSA methods

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

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

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

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

4 結(jié)論與評(píng)述

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

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

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

[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]魏群. 巖土工程中散體單元的基本原理、數(shù)值方法及實(shí)驗(yàn)研究[D]. 北京: 清華大學(xué), 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.