顆粒尺寸分散度對(duì)顆粒體系力學(xué)和幾何結(jié)構(gòu)特性的影響*

馮旭張國華 孫其誠

1)(北京科技大學(xué)物理系，北京 100083)

2)(清華大學(xué)，水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

(2013年5月4日收到;2013年6月18日收到修改稿)

1 引言

近年來，無序材料的力學(xué)和幾何結(jié)構(gòu)特性受到了較多關(guān)注.Silbert和Silbert[1]對(duì)比分析了三維無摩擦和有摩擦顆粒體系的靜態(tài)結(jié)構(gòu)因子S(k)，發(fā)現(xiàn)光滑顆粒體系在雙對(duì)數(shù)坐標(biāo)下S(k)曲線在低k區(qū)的線性行為，其斜率近似為1.Xu和Ching[2]研究了三維雙分散光滑顆粒體系中，發(fā)現(xiàn)靜態(tài)結(jié)構(gòu)因子強(qiáng)烈的依賴顆粒粒徑比，以前發(fā)現(xiàn)的低k部分的S(k)～k僅僅是顆粒的粒徑比接近1時(shí)的特殊情況.

與此同時(shí)，關(guān)于二維體系力學(xué)和幾何結(jié)構(gòu)特性的研究也取得了很大進(jìn)展.Han等[3]研究了二維膠體晶體的取向序關(guān)聯(lián)函數(shù)，發(fā)現(xiàn)當(dāng)處于Hexatic相-液相共存區(qū)時(shí)，系統(tǒng)的取向序關(guān)聯(lián)函數(shù)g6(r)～e-r/ξ6，且取向序關(guān)聯(lián)長度ξ6隨面密度ρ的變化規(guī)律與Kosterlitz-Thouless-Halperin-Nelson-Young理論的預(yù)測相符[4].Peng等[5]發(fā)現(xiàn)甚至當(dāng)二維膠體晶體融化到液相時(shí)，取向序關(guān)聯(lián)g6(r)曲線仍然以指數(shù)方式衰減.

近年來，關(guān)于二維體系結(jié)構(gòu)因子的研究也取得了很大進(jìn)展.Meyer等[6]發(fā)現(xiàn)在中間波矢區(qū)域二維聚合物溶液結(jié)構(gòu)因子隨k呈冪律標(biāo)度S(k)～kν，且隨著鏈長度N的增加，冪律指數(shù)ν從-2變到-3.Wen等[7]發(fā)現(xiàn)二維顆粒鏈的靜態(tài)結(jié)構(gòu)因子S(k)～k-2，與密集的聚合物溶液結(jié)構(gòu)相似.

在本文中，我們用顆粒離散元法(DEM)生成了由2048個(gè)二維圓形顆粒、邊壁壓強(qiáng)為1000Pa的系統(tǒng)，為了研究顆粒尺寸分散度s對(duì)系統(tǒng)性質(zhì)的影響，每個(gè)分散度下隨機(jī)生成100個(gè)位形.通過研究體系的配位數(shù)、剪切模量、靜態(tài)結(jié)構(gòu)因子、取向序關(guān)聯(lián)、力的累積分布等物理量隨尺寸分散度的變化規(guī)律，進(jìn)一步分析了體系的無序程度對(duì)二維顆粒體系力學(xué)和幾何結(jié)構(gòu)特性造成的影響.

2 數(shù)值模擬

在2 m×2 m的正方形盒子里，隨機(jī)放置質(zhì)量相等的2048個(gè)光滑圓盤顆粒，采用了周期性邊界條件.為了研究粒徑分散度s對(duì)系統(tǒng)力學(xué)性質(zhì)的影響，顆粒半徑在范圍內(nèi)等概率取值，其中d0=0.01 m為顆粒的平均直徑，s的取值范圍為[0，0.5].我們對(duì)每個(gè)粒徑分散度隨機(jī)生成100個(gè)壓強(qiáng)為1000 Pa的位形，本文所涉及的物理量都是這100個(gè)位形的統(tǒng)計(jì)平均值.文中，顆粒與顆粒相互作用為單邊線性彈簧，即當(dāng)顆粒i和j間距rij=ri-rj小于它們的半徑之和時(shí)存在相互作用，其中ri，rj分別為顆粒i和j的位置矢量.接觸力的法向分量由Fnij=-knnij給出，其中kn是法向接觸剛度，nij=rij/|rij|.在本文的模擬中，顆粒法向剛度系數(shù)為1.0×106N/m，切向剛度系數(shù)為0，不考慮重力.

具體產(chǎn)生位形的過程如下：首先，我們?cè)诤凶又须S機(jī)產(chǎn)生2048個(gè)粒徑很小的多分散顆粒，此時(shí)顆粒間沒有任何重疊，體系處于松弛狀態(tài).然后，在保持分散度固定的前提下使顆粒半徑增大，直至體系的體積分?jǐn)?shù)達(dá)到0.88.最后，利用半徑減小的辦法卸載，使得體系達(dá)到目標(biāo)壓強(qiáng)1000 Pa.其中，位形穩(wěn)定的判據(jù)是經(jīng)過1000個(gè)循環(huán)前后系統(tǒng)的能量差的比例小于1.0×10-15.

3 結(jié)果與分析

3.1 配位數(shù)和剪切模量

配位數(shù)是某一顆粒與周圍其他顆粒的接觸數(shù)目，它是表征系統(tǒng)幾何特征的重要物理量.為了研究尺寸分散度對(duì)系統(tǒng)幾何特征的影響，我們分析了平均配位數(shù)隨尺寸分散度的變化規(guī)律，如圖1(a)中所示.當(dāng)s＜0.1時(shí)，平均配位數(shù)的值隨分散度的增加迅速減小;而當(dāng)s＞0.1時(shí)，體系中顆粒的平均配位數(shù)趨于4.1，說明此時(shí)系統(tǒng)的尺寸無序程度對(duì)配位數(shù)基本沒有影響.

剪切模量是彈性材料承受τ剪應(yīng)力與產(chǎn)生的γ剪應(yīng)變的比值，是表征顆粒物質(zhì)抗剪切能力的力學(xué)指標(biāo).圖1(b)是系統(tǒng)平均剪切模量隨尺寸分散度s的變化曲線.可以看出，平均剪切模量和平均配位數(shù)隨s的變化趨勢一致.圖1(c)給出了平均剪切模量隨平均配位數(shù)隨s的線性變化規(guī)律，〈Z〉=3.5+3.3×10-5·G.由圖 1(a)，(b)，(c)可以看出，s=0.1可能是個(gè)分界點(diǎn)：當(dāng)體系的尺寸s＜0.1時(shí)，s越大(系統(tǒng)越無序)，系統(tǒng)的抗切應(yīng)變能力越弱;而當(dāng)s＞0.1時(shí)，系統(tǒng)的無序程度已經(jīng)對(duì)其幾何和力學(xué)性質(zhì)沒有任何影響.

圖1 平均配位數(shù)〈Z〉，剪切模量G隨分散度s的變化及二者的線性關(guān)系 (a)〈Z〉隨s的變化;(b)G隨s的變化;(c)〈Z〉隨G的變化

3.2 靜態(tài)結(jié)構(gòu)因子

靜態(tài)結(jié)構(gòu)因子S(k)是表征顆粒體系細(xì)觀結(jié)構(gòu)的典型參量[8-13]，定義為

我們模擬了s=0，0.001，0.005，0.008，0.02，0.1，0.2，0.3，0.4和0.5時(shí)的靜態(tài)結(jié)構(gòu)因子，結(jié)果如圖2所示.可以看出，在高k區(qū)域，尺寸分散度s＜0.1的S(k)曲線基本重合，并且在k′=2kπ/L=45，80，90附近各出現(xiàn)一個(gè)尖銳峰.隨著s增大，峰值趨于平緩，而且在k′=80和k′=100附近這兩個(gè)峰逐漸合并成一個(gè)平滑的峰.在低k區(qū)域(k＜20)，當(dāng)s＞0.1時(shí)，S(k)幾乎不隨k變化，S(k)的值隨著分散度的增大而增加;當(dāng)s≤0.02時(shí)，二維體系的S(k)曲線在3＜k′＜5區(qū)域遵從冪律標(biāo)度，S(k)～k-4/3，如圖2插圖所示.這一點(diǎn)符合二維線性聚合物鏈體系中在中間波矢范圍的標(biāo)度關(guān)系，S(k)～k-1/υ，與文獻(xiàn)[6，7，15]的結(jié)果類似，這暗示著二維單分散體系中存在顆粒鏈結(jié)構(gòu).

圖2 靜態(tài)結(jié)構(gòu)因子S(k)隨k的變化，插圖為單分散體系下S(k)曲線低k部分?jǐn)?shù)值擬合

為了研究維度對(duì)體系結(jié)構(gòu)因子的影響，我們生成了由10000個(gè)球形顆粒組成、壓強(qiáng)為10-4Pa的三維單分散顆粒體系，并計(jì)算了其結(jié)構(gòu)因子，如圖2插圖所示.可以發(fā)現(xiàn)，三維單分散體系的靜態(tài)結(jié)構(gòu)因子在低k區(qū)域滿足：S(k)～k，與文獻(xiàn)[1]中的結(jié)論一致.造成二維和三維體系S(k)在低k區(qū)域不同的原因可能是，在二維體系中更容易形成長程關(guān)聯(lián).

圖3為不同分散度下二維顆粒體系S(k)的云圖.可以看出，單分散體系(即s=0)表現(xiàn)出晶體的特征，其S(k)呈三角格子排列;而s=0.1的體系則表現(xiàn)出典型的多晶衍射圖樣的特征，衍射圖形呈現(xiàn)出明暗不均勻的同心圓環(huán);隨著粒徑多分散度的增大，圖形中從中心開始的第二個(gè)圓環(huán)由正六邊形演變?yōu)閳A形，圓環(huán)的徑向?qū)挾纫搽S之增加，而且看不到明顯的點(diǎn);當(dāng)s≥0.3時(shí)，同心圓環(huán)變得更寬且彼此合并，其外圍輪廓逐漸變得模糊，表明系統(tǒng)隨著s的增大越來越無序;當(dāng)s=0.5時(shí)，只能看到一個(gè)比較清晰的圓環(huán).

3.3 取向序關(guān)聯(lián)

為了量化顆粒i的取向序，引入鍵取向序參數(shù)Ψ6i[16]：

式中，ni表示顆粒i的最近鄰顆粒的數(shù)目，θij表示顆粒i和其近鄰顆粒j之間的極角.為了進(jìn)一步量化取向序的空間關(guān)聯(lián)，本文計(jì)算了顆粒體系的取向序關(guān)聯(lián)函數(shù)g6(r)[3,17,18]，定義為

圖3 不同分散度下靜態(tài)結(jié)構(gòu)因子S(k)

圖4 取向序關(guān)聯(lián)g6(r)隨r的變化及擬合參數(shù)ξ6隨分散度s的變化規(guī)律 (a)g6(r)隨r的變化;(b)ξ6隨s的變化;實(shí)線為ξ6∝14.2e-s/0.2

圖4(a)給出了不同分散度下的取向序關(guān)聯(lián)函數(shù)g6(r).由圖4(a)可知，取向序關(guān)聯(lián)函數(shù)呈現(xiàn)隨著r的增大幅值逐漸減小的振蕩行為.其中，g6(r)的極小可能與顆粒沒有占據(jù)晶格的位置有關(guān)，而g6(r)振蕩衰減的長度范圍則對(duì)應(yīng)取向關(guān)聯(lián)的長度.為了得到取向序關(guān)聯(lián)長度，我們對(duì)曲線進(jìn)行了e指數(shù)擬合，如圖4(a)所示.虛線為擬合曲線：g6(r)∝ae-r/ξ6，其中ξ6是序關(guān)聯(lián)長度.由圖4(a)可知，s＞0.1的g6(r)曲線均呈e指數(shù)衰減，表明s＞0.1的顆粒體系表現(xiàn)出與液體類似的短程序[3,19].圖4(b)給出了序關(guān)聯(lián)長度隨分散度變化的曲線，顯然，隨著s的增加序關(guān)聯(lián)長度近似以e指數(shù)規(guī)律衰減，如圖4(b)中的實(shí)線所示.總之，隨著粒徑分散度s增大，ξ6減小，而且減小的幅度越來越慢，變化比較連續(xù).說明隨著粒徑多分散度s從0.1到0.5，體系越來越無序，而且它的變化并不是一個(gè)突變，而是一個(gè)結(jié)構(gòu)連續(xù)變化的過程.此時(shí)，體系無序程度對(duì)顆粒間取向序關(guān)聯(lián)的影響比較明顯.

3.4 力的累積分布

在靜態(tài)顆粒排布中，顆粒間接觸力形成高度各向異性的力網(wǎng)絡(luò)，其基本特征可以用接觸力的概率密度函數(shù)p(f)表征，其中f=F/〈F〉為歸一化力.近年來，大量的實(shí)驗(yàn)和模擬關(guān)注大f處p(f)的分布，由于實(shí)驗(yàn)上的困難，對(duì)于小f處p(f)分布的研究還不是很多[20].我們計(jì)算了不同分散度下顆粒體系中的力累積分布G(f)，G(f)=圖5可以看出，當(dāng)f＜2時(shí)，隨著f的增大，G(f)值也增大;當(dāng)f＞2后，G(f)趨近于1.而且，不同分散度下的G(f)曲線幾乎重合，暗示系統(tǒng)的無序程度對(duì)力累積分布幾乎沒有影響.

圖5 力的累積分布G(f)隨 f的變化

4 結(jié)論

采用顆粒離散元方法生成了具有不同尺寸分散度的二維顆粒體系，進(jìn)一步研究了尺寸分散度對(duì)體系力學(xué)和幾何結(jié)構(gòu)特征的影響，得到如下結(jié)論：1)二維顆粒體系，平均配位數(shù)和平均剪切模量隨著s的增大而減小;當(dāng)s＞0.1時(shí)，二者都基本保持為定值，說明隨著s的增大，體系的無序程度增加;2)當(dāng)分散度s≤0.02時(shí)，靜態(tài)結(jié)構(gòu)因子低k部分的曲線基本重合;當(dāng)s＞0.1時(shí)，隨著s的增大，S(k)的值也均勻增加，尤其是s=0時(shí)，二維顆粒體系S(k)的低k區(qū)域行為與聚合物顆粒鏈的類似，S(k)～k-1/υ(υ=3/4)，這一點(diǎn)不同于三維單分散顆粒體系(其S(k)曲線低k部分斜率近似等于1)，暗示二維單分散顆粒體系(包括分散度較小的體系)結(jié)構(gòu)上與顆粒鏈相似;3)不同分散度下，取向序關(guān)聯(lián)函數(shù)g6(r)曲線的峰值均滿足e指數(shù)關(guān)系，序關(guān)聯(lián)長度ξ6也隨尺寸分散度的增大而減小;4)力的累積分布G(f)不隨尺寸分散度的變化而變化，說明它基本不受系統(tǒng)無序程度的影響.

[1]Silbert L E，Silbert M 2009Phys.Rev.E 80 041304

[2]Xu N，Ching E S C 2010Soft Matter6 2944

[3]Han Y，Ha N Y，Alsayed A M，Yodh A G 2008Phys.Rev.E 77 041406

[4]Artoni R，Santomaso A C，Gabrieli F，Tono D，Cola S 2013Phys.Rev.E 87 032205

[5]Peng Y，Wang Z，Alsayed A M，Yodh A G，Han Y 2010Phys.Rev.Lett.104 205703

[6]Meyer H，Schulmann N，Zabel J E，Wittmer J P 2011Comput.Phys.Commun.182 1949

[7]Wen P P，Zheng N，Li L S，Li H，Sun G，Shi Q F 2012Phys.Rev.E 85 031301

[8]Yang J K，Schreck C，Noh H，Liew S F，Guy M I，O’Hern C S，Cao H 2010Phys.Rev.A 82 053838

[9]Xu W S，Sun Z Y，An L J 2012J.Chem.Phys.137 104509

[10]Berthier L，Chaudhuri P，Coulais C，Dauchot O，Sollich P 2011Phys.Rev.Lett.106 120601

[11]Paulus M，Gutt C，Tolan M 2008Phys.Rev.B 78 235419

[12]Donev A，Stillinger F H，Torquato S 2005Phys.Rev.Lett.95 090604

[13]Torquato S，Stillinger F H 2003Phys.Rev.E 68 041113

[14]Warr S，Hansen J P 1996Europhys.Lett.36 589

[15]Maier B，R¨adler J O 1999Phys.Rev.Lett.82 1911

[16]Schreck C F，O’Hern C S，Silbert L E 2011Phys.Rev.E 84 011305

[17]Agarwal U，Escobedo F A 2012Soft Matter8 5916

[18]Prestipino S，Saija F，Giaquinta P V 2011Phys.Rev.Lett.106 235701

[19]Bakker A F，Bruin C，Hilhorst H J 1984Phys.Rev.Lett.52 449

[20]Charbonneau P，Corwin E I，Parisi G，Zamponi F 2012Phys.Rev.Lett.109 205501