振動形成的單一尺寸硬球晶的結構表征及其動力學的離散元模擬

安希忠，李長興，楊潤宇，余艾冰，曲德毅

(1.東北大學 材料與冶金學院，沈陽 110004;2.新南威爾士大學 材料科學與工程學院顆粒系統模擬中心，澳大利亞 悉尼新南威爾士 2052;3.沈陽鼓風機集團股份有限公司，沈陽 110869)

振動形成的單一尺寸硬球晶的結構表征及其動力學的離散元模擬

安希忠1，李長興1，楊潤宇2，余艾冰2，曲德毅3

(1.東北大學 材料與冶金學院，沈陽 110004;2.新南威爾士大學 材料科學與工程學院顆粒系統模擬中心，澳大利亞 悉尼新南威爾士 2052;3.沈陽鼓風機集團股份有限公司，沈陽 110869)

采用離散元(DEM)數值方法結合物理實驗，分別實現了單一尺寸硬球在三維(3D)間歇振動及批量加料條件下堆積的有序化(結晶)，通過對晶體的結構表征、數值仿真中的微觀性能如配位數(CN)、徑向分布函數(RDF)、角分布函數(ADF)、Voronoi/Delaunay孔的尺寸分布、力的網絡分析以及晶體形成過程中的動力學分析，得出如下結論:(1)數值仿真和物理實驗中所獲得的均為{111}取向FCC晶體，晶體中粒子的最大堆積密度可達0.74;(2)晶體中粒子的配位數分布在CN=12處出現峰值，表明每個粒子具有12個近鄰;RDF和ADF分布表明堆積粒子之間分別在徑向具有長程相關性及角度上具有相關性，這是晶態結構的典型特征;Voronoi/Delaunay孔的尺寸分布與非晶態相比表現為高且窄的曲線，且向左移，表明所獲得的硬球晶中孔隙很小且分布均勻;粒子間力的網絡更直觀地表明了晶體的結構，同時還指出在FCC硬球晶中分別存在著面內的力和面間的力，其中后者小于前者;(3)晶體形成過程中的速度場和力場等動力學信息的變化表明粒子在堆積過程中一旦形成有序，它將作為一個整體運動，并為后續粒子的有序化提供模板，起到了結晶的晶核的作用.

顆粒堆積;致密化;振動;離散元模擬;硬球晶

無論是在工業生產還是在科學研究中，對硬球堆積致密化的研究一直是一個重要的課題［1，2］.因為對于許多物理學家和材料學家而言，硬球的堆積可用于研究多種系統結構(直接研究這些系統的結構是相當困難的)的有效出發點.在過 去 的 幾 十 年 里，從 Bernal［3］，Scott［4］和Finney［5］等人首次利用硬球的堆積模擬簡單液體結構取得成功以來，許多科學工作者又成功實現了硬球堆積非晶態(隨機堆積)的轉變［6～13］，即從初始的松裝密度約0.60轉變為非晶態結構中硬球的最大堆積密度約0.64，在這個過程中，外部施加的振動起到了至關重要的作用.采用振動及合理地控制工藝參數也可以實現粒子的晶態堆積結構.這一點無論是從物理實驗上還是從數值計算中都得到證實［14～23］.但由于數值及物理實驗中的特殊性，如現有的物理實驗中對硬球結晶的動力學分析的難度，大多數數值仿真都是基于幾何的模擬而忽略了動力學的范疇［21］，從而導致所獲得的結果很難再現實際的過程，故此對于完美硬球晶的形成過程及其動力學還有待于進一步的研究.

使用基于分子動力學的離散元(DEM)數值仿真，可以克服其他數值仿真方法的不足，從而從動力學的角度真實再現粒子堆積的致密化，DEM方法的有效性已經在許多科學研究中得到了證實［12，21，23～25］.在以前的工作中［21］，使用 DEM 模擬，我們已經成功實現了振動條件下完美硬球晶的獲得.本文將在此基礎上，對該硬球晶的宏觀和微觀性能及其形成過程的動力學進行分析，以期找出硬球振動堆積過程中有序化的規律，為實際完美硬球晶的獲得提供參照.

1 模擬方法及條件

1.1 數值仿真方法

在DEM模擬中，每個粒子有兩種運動形式，即平動和轉動.根據牛頓第二定律，平動和轉動方程可分別表示為

其中 mi，vi，ωi和 Ii分別表示粒子 i的質量、平動速度、角速度及轉動慣量.Ri是由粒子i的中心指向接觸點的矢量，其大小等于粒子i的半徑Ri.法向接觸力和切向接觸力分別表示為

1.2 模擬條件

模擬由在一個長方體容器中隨機產生一批無重疊的等同球形粒子開始，容器底面為邊長為12 d(d表示粒子直徑，模擬中所用的顆粒子直徑d=1 cm)的正方形，為了避免容器壁的影響，在水平方向上采用了周期性的邊界條件.初始時，隨機產生的粒子在重力作用下自由下落，向下運動的過程中粒子與其近鄰之間互相碰撞并來回彈跳躍.該動力學過程一直持續到所有粒子都達到穩定狀態，此時由于能量完全耗散的結果，所有粒子的運動速度為零，這一過程在一秒鐘內完成.一秒鐘后，容器在相互垂直的3個方向上進行三維(3D)周期性正弦振動，每個方向上振動的控制方程為:Z(t)=Asin［ω(t－t0)］，其中 Z 表示容器的位移，A和ω分別表示振動的振幅和頻率.t0表示振動開始時間(本工作中t0=1.0 s)，如果是1D振動，則只在豎直方向上施加周期性正弦振動.從位移方程很容易推導出振動的加速度方程:a(t)= －Aω2sin［ω(t－t0)］.本研究中采用的是間歇振動、批量加料，即從1 s開始，粒子先振動1 s，此過程中以一定的批量加料，然后停止1 s使粒子達到穩定狀態，接著再振動1 s，此時再加入下一批料，之后再停止1 s，……，這樣的過程一直持續到振動結束.實驗中所用的粒子為玻璃球，為了簡化起見，取3D振動中3個方向上的振幅和頻率分別相同.表1中列出了模擬中所用的參數及其所對應的取值，其中假定粒子與容器具有相同的性能.

表1 模擬中所用的參數Table 1 The parameters used in the simulation

2 研究結果及分析

2.1 硬球晶的實現

圖1給出了硬球在3D間歇振動、批量加料的條件下，結晶過程中所對應的宏觀性能如堆積密度ρ隨振動時間t的演變及最終所獲得的局部完美硬球晶體結構.圖中的內插圖是對密度變化曲線的局部放大，其中振幅A=0.1 d，振動頻率ω=200 rad/s，每次加料的批量為Nb=98個粒子/批［21］.可以看出在ρ－t曲線上存在著兩個區域，第一個區域一直持續到28.0 s，在這期間堆積密度隨時間增加.其中曲線上的起伏表示振動過程中對粒子堆積周期性壓縮與松馳的結果.在第二個區域，振動使堆積密度沿著一個相對穩定的值變化.當振動停止時，球逐漸停下來形成穩定堆積，最終的堆積密度可達約0.74，結構表征發現其對應的是{111}取向的FCC晶體.圖1中的DEM數值仿真結果由圖2(其中的內插圖為去除表面堆積粒子后所得到的內部有序結構)中的物理實驗結果所驗證［22］.圖2中的結果是在5個不同尺寸容器中展開的單一尺寸玻璃球在3D間歇振動(振幅A=0.4 d，振動頻率ω=55 rad/s)及批量加料(Nb=1層/批)條件下得到的，為了消除容器壁的影響，我們對每種情況下不同容器中的實驗結果進行了外推，結果得到在無限大容器中粒子的最大堆積密度可達約0.74，且對最終的堆積結構表征發現為{111}取向的FCC晶體，即在豎直的{111}方向上六方密排面的堆垛順序為ABCABC……，體現了數值仿真結果和物理實驗結果很好的一致性.

圖1 3D間歇振動、批量加料時，晶態結構(內插圖為原圖的局部放大)的實現(左)及最終所對應的部分晶體(右)，振動條件為:A=0.1 d，ω=200 rad/s，加料批量Nb=98粒子/批，圖中A1為每一個振動周期內的振動開始點(在此處開始加料)，A2為振動結束點［21］Fig.1 The realization of crystalline structure(left，the inset figure is the partial enlargement)and part of the final crystal(right)obtained under 3D interval vibration and batch－wised feeding，where the vibration condition is:A=0.1 d，ω =200 rad/s，batch of feeding Nb=98 particles/batch;A1represents the vibration starting point in each vibration period where the particles are fed at this moment，and A2is the end of vibration［21］

2.2 硬球晶的微觀結構表征

圖2 物理實驗中，3D間歇振動、批量加料且A=0.4 d，ω=55 rad/s，Nb=1層/批時容器尺寸D(容器直徑)對堆積密度的影響，內插圖為去除表面粒子而得到的有序結構［22］Fig.2 The effects of container size on the packing density in physical experiments with 3D interval vibration and batch－wised feeding when A=0.4 d，ω =55 rad/s，and Nb=1 layer/batch，where the inset figures are the ordered structure obtained when the particles on the packing surface are removed［22］

圖3 晶態結構的配位數分布Fig.3 Coordination number distribution of the crystalline structure

除了硬球晶宏觀性能的演變外，我們也對所獲得的硬球晶的微觀性能進行了表征，這些微觀性能包括:配位數(CN)、徑向分布函數(RDF)、角分布函數(ADF)、Voronoi/Delaunay孔的尺寸分布以及晶體內力的網絡結構，它們都是形成穩態的硬球晶后得到的.圖3給出了硬球晶的配位數的分布，CN被定義為與給定粒子“接觸”的粒子數，本工作中取臨界距離為1.005 d，即如果兩個粒子之間的距離小于該值，就被認為是接觸.與非晶態結構CN分布的峰值為7不同地是［5，12，23，24］，圖3 中的 CN 分布的峰值對應于12，這是FCC晶體對應的配位數，表明我們所獲得的結構是有序的.

與配位數CN相比，徑向分布函數RDF和角分布函數ADF是表征晶體結構的更直接的方式.RDF定義為堆積中與特定粒子中心給定距離范圍內發現一個粒子中心的概率，它包含了堆積粒子間徑向相關性的有用信息;相對地，ADF定義為堆積中與特定粒子中心給定角度范圍內發現一個粒子中心的概率，它包含了堆積粒子間角度相關性的有用信息.圖4給出了我們所獲得的最終硬球晶的RDF和ADF分布，與非晶態(隨機密排)RDF 分布上的分叉第二峰不同［5，12，23，24］，晶態結構的RDF分布還在r＞2d的許多相關位置上出現了各種孤立的峰，這些在較遠位置上的峰的存在表明了堆積體內粒子之間位置的長程相關性，這是晶態結構的典型特征［21］.同樣，我們所獲得的晶態結構的ADF也與非晶態結構的ADF分布完全不同，后者在圖上表現在角度在60(°)～180(°)之間的一條連續的曲線，而前者的ADF只在60(°)、90(°)、120(°)和180(°)的角度上存在明顯的峰，這表明堆積體內粒子之間存在著角度的相關性，這進一步表明了我們所獲得的硬球堆積結構為晶態而不是非晶態［23］.

對所獲得硬球晶局部結構更具體的分析來自于 Voronoi和 Delaunay 對空間的劃分［29～33］，它比圖4中的一維RDF和ADF能夠提供更豐富的信息［34］.Voronoi空間劃分是通過 Voronoi多面體來實現的，它定義為堆積體內垂直平分給定粒子中心與其他粒子中心之間線段的面所圍成的最小封閉凸多面體［30］，如果連接堆積體內每個粒子與其近鄰粒子中心，則會構成一個空間網絡，該網絡則稱為Delaunay的空間劃分.網絡中的四面體單體形成了對空間的完全填充.圖5中給出了Voronoi(內插圖為Delaunay)的孔的尺寸分布，我們將Voronoi或Delaunay單體內孔的體積V除以粒子的體積Vp而使其無量綱化.由圖可見，無論是Voronoi還是Delaunay孔的尺寸分布，都呈現出高且窄的趨勢，并且明顯小于粒子的體積.這表明晶態結構中孔的尺寸很小且分布很均勻;這一點與非晶態結構中呈現右移的寬且矮的Voronoi/Delaunay孔的尺寸分布完全不同［12］.

為了更直觀的表征硬球晶的內部結構，我們也給出了最終部分硬球晶的力的網絡結構，如圖6所示.圖中的“球”表示晶體中堆積的粒子，“棍”表示粒子之間的力，“棍”的粗細與力的大小成正比.從圖中除了可以很清晰地看出所形成的是{111}取向FCC晶體以外，更重要地是我們可以得到力在晶體中傳播的信息，即在晶體中存在著面內的力和面間的力.在堆積中的每個球都可以看作是一個“外部載荷”，在這些“外部載荷”的作用下，力可以沿著三角架的方向傳遞，從而形成面間的網絡，這與Mueggenburg等人［35］的實驗結果相一致，所不同的是，在我們的結果有存在著較大的面內的力.

2.3 硬球晶形成的動力學分析

除了前面對最終硬球晶的宏觀性能和微觀性能進行表征外，我們還從動力學的角度研究了硬球晶的形成過程.硬球結晶過程中的速度場和力場的變化分別如圖7和8所示.這兩個圖中顯示地是三維矢量場中的二維視場，其中視角的方向為沿著X軸的方向.從這兩個圖中我們可以發現堆積過程中出現的幾個特征.首先，是振動過程所帶來的對粒子的激發及松馳過程，其中對粒子的激發發生在t=53.318 6s，在這個過程中粒子的運動速度和力都很大.而粒子的松馳發生在t=47.044 2 s，此時粒子的速度和力逐漸減弱，漸漸形成穩定結構.不管是激發還是松馳，粒子運動最劇烈地是在表面，這可以從速度和力的變化反映出來，這是因為表面粒子在振動堆積致密化過程中由于缺少約束而處于流化狀態，隨著外部輸入能量的停止，它們將找到各自的穩定位置從而形成有序的結構.另一方面，在有序堆積過程中，先堆積的粒子的運動狀態基本相同，也就是說在底層形成有序結構后，其中粒子的速度和力的變化不大，因此，在振動過程中其是作為一個整體運動的，它作為后續粒子有序化堆積的模板，起到了結晶的晶核的作用，這一點已經在我們的數值仿真［23］和物理實驗［22］中得到了證實.

圖4 晶態結構的徑向分布函數(左)和角分布函數(右)Fig.4 Radial distribution function(left)and angular distribution function(Right)for the crystalline structure

3 結論

采用離散元(DEM)數值方法模擬了單一尺寸球形粒子在三維(3D)間歇振動及批量加料條件下堆積的有序化(結晶)，并對所獲得的硬球晶的宏觀性能如堆積密度及微觀性能如配位數CN、徑向分布函數 RDF、角分布函數 ADF、Voronoi/Delaunay孔的尺寸分布及力的網絡結構進行了表征，同時從動力學角度分析了硬球堆積有序化過程中的速度場和力場的變化，得出結論如下:

(1)對單一尺寸粒子在3D間歇振動及批量加料條件下堆積的DEM模擬可以獲得最大堆積密度達0.74的硬球晶，同時，該過程已得到了物理實驗的驗證.

(2)硬球晶的配位數(CN)分布在CN=12處出現峰值，結構表征發現對應的為FCC晶體;徑向分布函數(RDF)和角分布函數(ADF)的表征發現晶體中的粒子分別存在徑向的長程相關性及角度的相關性，這是晶態結構區別于非晶態結構的典型特征;Voronoi/Delaunay孔的尺寸分布表現為左移的高且窄的峰，表明堆積結構中的孔的尺寸很小且分布均勻;力的網絡直觀的表現了晶體的結構，同時表明所形成的{111}取向FCC硬球晶中，在面內和面間均有力的存在，且面內的力大于面間的力.

(3)通過對粒子振動堆積有序化過程中的動力學信息如速度場和力場變化的分析，發現硬球在結晶的過程中，一旦有序結構形成，該有序結構中的粒子在后面的振動過程中保持相同的運動狀態，即作為一個整體運動，它為后面有序化的繼續提供了模板，起到了晶核的作用.

［1］German R M. Particle packing characteristics［M］.Princeton，New Jersey:Metal Powder Industries，1989:1－18.

［2］Bideau D，Hansen A.Disorder and granular media，random materials and process［M］.Edited by Stanley H E，Guyon E，North － Holland.Amsterdam:Elsevier Science Publishers，1993:1－25.

［3］Bernal J D.A geometrical approach to structure of liquids［J］.Nature，1959，183(4655):141－147.

［4］Scott G D.Packing of spheres［J］.Nature，1960，188(4754):908－909.

［5］Finney J L.Random packings and the structure of simple liquids:I the geometry of random close packing ［J］.Proc Roy Soc Lond A，1970，319:479－493.

［6］Woodcock L V.Glass transition in the hard－sphere mode［J］.J Chem Soc Faraday II，1976，72:1667－1672.

［7］Finney J L.Modelling the structures of amorphous metals and alloys［J］.Nature，1977，266:309 －314.

［8］Zallen R.The physics of amorphous solids［M］.New York:Wiley，1983:1－15.

［9］Knight J B，Fandrich C G，Lau C N，et al.Density relaxation in a vibrated granular material［J］.Phys Rev E，1995，51(5):3957－3963.

［10］Torquato S.Glass transition:hard knock for thermodynamics［J］.Nature，2000，405:521 －523.

［11］Stachurski Z H.Definition and properties of ideal amorphous solids［J］.Phys Rev Lett，2003，90(15):155502.

［12］An X Z，Yang R Y，Dong K J，et al.Micromechanical simulation and analysis of one－dimensional vibratory sphere packing［J］.Phys Rev Lett，2005，95:205502 －1 －4.

［13］An X Z，Li C X，Yang R Y，et al.Experimental study of the packing of mono－sized spheres subjected to one－dimensional vibration［J］.Powder Technol，2009，196:50 －55.

［14］Owe Berg T G，Mcdonald R L，Trainor R J.The packing of spheres［J］.Powder Technol，1969/1970，3:183 －188.

［15］Pouliquen O，Nicolas M，Weidman P D.Crystallization of non－Brownian spheres under horizontal shaking［J］.Phys Rev Lett，1997，79:3640 －3643.

［16］Van Blaaderen A，Ruel R，Wiltzius P.Template－directed colloidal crystallization［J］.Nature，1997，385:321 －324.

［17］Blair D L，Mueggenburg N W，Marshall A H，et al.Force distributions in three－dimensional granular assemblies:effects of packing order and interparticle friction［J］.Phys Rev E，2001，63:041304.

［18］Nahmad－Molinari Y，Ruiz－Suárez J C.Epitaxial growth of granular single crystals ［J］. Phys Rev Lett， 2002，89:264302.

［19］Kansal A R，Torquato S，Stillinger F H.Diversity of order and densities in jammed hard － particle packings［J］.Phys Rev E，2002，66:041109－1－8.

［20］Spannuth M J，Mueggenburg N W，Jaeger H M，et al.Stress transmission through three－dimensional granular crystals with stacking faults［J］.Granular Matter，2004，6:215 －219.

［21］Yu A B，An X Z，Zou R P，et al.Self－assembly of particles for densest packing by mechanical vibration［J］.Phys Rev Lett，2006，97:265501.

［22］Li C X，An X Z，Yang R Y，et al.Experimental study on the packing of uniform spheres under three－dimensional vibration［J］.Powder Technol，2011，208(3):617 － 622.

［23］An X Z，Yang R Y，Dong K J，et al.DEM study of crystallization of monosized spheres under mechanical vibrations［J］.Computer Physics Communications，2011，182:1989 －1994.

［24］Yang R Y，Zou R P，Yu A B.Computer simulation of the packing of fine particles［J］.Phys Rev E，2000，62(3):3900－3908.

［25］An X Z，Yang R Y，Zou R P，et al.Effect of vibration condition and inter－particle frictions on the packing of uniform spheres［J］.Powder Technol，2008，188(2):102－109.

［26］Martin C L， Bouvard D， Shima S. Studyofparticle rearrangement during powder compaction bythe discrete element method［J］.J Mech Phys Solids，2003，51(4):667－693.

［27］Schwager T，P?schel T.Coefficient of normal restitution of viscous particles and cooling rate of granular gases［J］.Phys Rev E，1998，57(1):650－654.

［28］Langston P A，Tüzün U，Heyes D M.Discrete element simulation of granular flow in 2D and 3D hoppers:dependence of discharge rate and wall stress on particle interactions［J］.Chem Eng Sci，1995，50(6):967 －987.

［29］Oger L，Gervois A，Troadec J P，et al.Voronoi tessellation of packing of spheres:topological correlation and statistics［J］.Phil Mag B，1996，74:177－197.

［30］Yang R Y，Zou R P，Yu A B.Voronoi tessellation of the packing of fine uniform spheres［J］.Phys Rev E，2002，65:041302.

［31］Aste T，Saadatfar M，Senden T J.Geometrical structure of disorderedspherepackings ［J］. Phys Rev E，2005，71:061302.

［32］Yang R Y，Zou R P，Yu A B，et al.Pore structure of the packing of fine particles［J］.J Colloid Interface Sci，2006，299:719－725.

［33］An X Z.Evolution of Voronoi/Delaunay characterized micro structure with transition from loose to dense sphere packing［J］.Chin Phys Lett，2007，24(8):2327 －2330.

［34］Montoro J C G，Abascal J L F.The Voronoi polyhedra as tools for structure determination in simple disordered systems［J］.J Phys Chem，1993，97:4211－4215.

［35］Mueggenburg N W，Jaeger H M，Nagel S R.Stress transmission through three－dimensional ordered granular arrays［J］.Phys Rev E，2002，66:031304－1.

Structure characterization on the vibration induced mono－sized hard sphere crystal and its dynamics by DEM modeling

AN Xi－zhong1，LI Chang－xing1，YANG Run－yu2，YU Ai－bing2，QU De－yi3

(1.School of Materials and Metallurgy，Northeastern University，Shenyang 110004，China;2.Laboratory for Simulation and Modeling of Particulate System，School of Materials Science and Engineering，University of New South Wales，Sydney 2052，Australia;3.Shenyang Blower Works Group，Co.，Ltd.，Shenyang 110869，China)

In this paper，the ordering(Crystallization)of mono－sized sphere packing under three－dimensional(3D)interval vibration and batch－wised feeding was studied with numerical discrete element method(DEM)combined with the physical experiments.Through the characterization on the crystalline structure and the analysis on the micro－properties obtained numerically such as coordination number(CN)，radial distribution function(RDF)，angular distribution function(ADF)，Voronoi/Delaunay void size distribution，force network and the dynamics formed during ordering，following conclusions can be drawn:(1)The crystalline obtained both numerically and physically all indicates{111} －oriented FCC structure，where the maximum packing density can reach 0.74;(2)There is a peak at CN=12 in the coordination number distribution curve，which implies that each particle in the hard sphere crystal has 12 neighbors in contact;the distributions for RDF and ADF show the long range correlation in radial distance and angle correlation among particles，which are the typical characteristics of the crystalline structure;the Voronoi/Delaunay void size distribution indicates the higher and narrower curve shifting to the left compared with non－crystal state，which means that the voids in crystalline structure are small and uniform;the force network among particles can identify the crystalline structure more intuitively，meanwhile，it’s also indicated that the in － plane forces and inter－ plane forces coexist in the crystal，and the latter is smaller than the former;(3)The variation of dynamic information such as the velocity field and the force field all indicates that once the ordered structure formed during crystallization，it will move as a whole part，which served as the template for the other particles，i.e.as the nucleus inthe crystallization.

particle packing;densification;vibration;DEM simulation;hard sphere crystal

TF 122;TB 44;O 756

A

1671－6620(2011)04－0318－07

2011－09－20.

國家自然基金 (50974040)，中央高校基本科研業務費 (N100402009)，澳大利亞研究委員會基金 (DP0665276).

安希忠 (1973—)，男，遼寧蓋州人，東北大學副教授，E－mail:anxz@mail.neu.edu.cn