冯旭张国华 孙其诚

1)(北京科技大学物理系，北京 100083)

2)(清华大学，水沙科学与水利水电工程国家重点实验室，北京 100084)

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

1 引言

近年来，无序材料的力学和几何结构特性受到了较多关注.Silbert和Silbert[1]对比分析了三维无摩擦和有摩擦颗粒体系的静态结构因子S(k)，发现光滑颗粒体系在双对数坐标下S(k)曲线在低k区的线性行为，其斜率近似为1.Xu和Ching[2]研究了三维双分散光滑颗粒体系中，发现静态结构因子强烈的依赖颗粒粒径比，以前发现的低k部分的S(k)∼k仅仅是颗粒的粒径比接近1时的特殊情况.

与此同时，关于二维体系力学和几何结构特性的研究也取得了很大进展.Han等[3]研究了二维胶体晶体的取向序关联函数，发现当处于Hexatic相-液相共存区时，系统的取向序关联函数g6(r)∼e-r/ξ6，且取向序关联长度ξ6随面密度ρ的变化规律与Kosterlitz-Thouless-Halperin-Nelson-Young理论的预测相符[4].Peng等[5]发现甚至当二维胶体晶体融化到液相时，取向序关联g6(r)曲线仍然以指数方式衰减.

近年来，关于二维体系结构因子的研究也取得了很大进展.Meyer等[6]发现在中间波矢区域二维聚合物溶液结构因子随k呈幂律标度S(k)∼kν，且随着链长度N的增加，幂律指数ν从-2变到-3.Wen等[7]发现二维颗粒链的静态结构因子S(k)∼k-2，与密集的聚合物溶液结构相似.

在本文中，我们用颗粒离散元法(DEM)生成了由2048个二维圆形颗粒、边壁压强为1000Pa的系统，为了研究颗粒尺寸分散度s对系统性质的影响，每个分散度下随机生成100个位形.通过研究体系的配位数、剪切模量、静态结构因子、取向序关联、力的累积分布等物理量随尺寸分散度的变化规律，进一步分析了体系的无序程度对二维颗粒体系力学和几何结构特性造成的影响.

2 数值模拟

在2 m×2 m的正方形盒子里，随机放置质量相等的2048个光滑圆盘颗粒，采用了周期性边界条件.为了研究粒径分散度s对系统力学性质的影响，颗粒半径在范围内等概率取值，其中d0=0.01 m为颗粒的平均直径，s的取值范围为[0，0.5].我们对每个粒径分散度随机生成100个压强为1000 Pa的位形，本文所涉及的物理量都是这100个位形的统计平均值.文中，颗粒与颗粒相互作用为单边线性弹簧，即当颗粒i和j间距rij=ri-rj小于它们的半径之和时存在相互作用，其中ri，rj分别为颗粒i和j的位置矢量.接触力的法向分量由Fnij=-knnij给出，其中kn是法向接触刚度，nij=rij/|rij|.在本文的模拟中，颗粒法向刚度系数为1.0×106N/m，切向刚度系数为0，不考虑重力.

具体产生位形的过程如下：首先，我们在盒子中随机产生2048个粒径很小的多分散颗粒，此时颗粒间没有任何重叠，体系处于松弛状态.然后，在保持分散度固定的前提下使颗粒半径增大，直至体系的体积分数达到0.88.最后，利用半径减小的办法卸载，使得体系达到目标压强1000 Pa.其中，位形稳定的判据是经过1000个循环前后系统的能量差的比例小于1.0×10-15.

3 结果与分析

3.1 配位数和剪切模量

配位数是某一颗粒与周围其他颗粒的接触数目，它是表征系统几何特征的重要物理量.为了研究尺寸分散度对系统几何特征的影响，我们分析了平均配位数随尺寸分散度的变化规律，如图1(a)中所示.当s＜0.1时，平均配位数的值随分散度的增加迅速减小;而当s＞0.1时，体系中颗粒的平均配位数趋于4.1，说明此时系统的尺寸无序程度对配位数基本没有影响.

剪切模量是弹性材料承受τ剪应力与产生的γ剪应变的比值，是表征颗粒物质抗剪切能力的力学指标.图1(b)是系统平均剪切模量随尺寸分散度s的变化曲线.可以看出，平均剪切模量和平均配位数随s的变化趋势一致.图1(c)给出了平均剪切模量随平均配位数随s的线性变化规律，〈Z〉=3.5+3.3×10-5·G.由图 1(a)，(b)，(c)可以看出，s=0.1可能是个分界点：当体系的尺寸s＜0.1时，s越大(系统越无序)，系统的抗切应变能力越弱;而当s＞0.1时，系统的无序程度已经对其几何和力学性质没有任何影响.

图1 平均配位数〈Z〉，剪切模量G随分散度s的变化及二者的线性关系 (a)〈Z〉随s的变化;(b)G随s的变化;(c)〈Z〉随G的变化

3.2 静态结构因子

静态结构因子S(k)是表征颗粒体系细观结构的典型参量[8-13]，定义为

我们模拟了s=0，0.001，0.005，0.008，0.02，0.1，0.2，0.3，0.4和0.5时的静态结构因子，结果如图2所示.可以看出，在高k区域，尺寸分散度s＜0.1的S(k)曲线基本重合，并且在k′=2kπ/L=45，80，90附近各出现一个尖锐峰.随着s增大，峰值趋于平缓，而且在k′=80和k′=100附近这两个峰逐渐合并成一个平滑的峰.在低k区域(k＜20)，当s＞0.1时，S(k)几乎不随k变化，S(k)的值随着分散度的增大而增加;当s≤0.02时，二维体系的S(k)曲线在3＜k′＜5区域遵从幂律标度，S(k)∼k-4/3，如图2插图所示.这一点符合二维线性聚合物链体系中在中间波矢范围的标度关系，S(k)∼k-1/υ，与文献[6，7，15]的结果类似，这暗示着二维单分散体系中存在颗粒链结构.

图2 静态结构因子S(k)随k的变化，插图为单分散体系下S(k)曲线低k部分数值拟合

为了研究维度对体系结构因子的影响，我们生成了由10000个球形颗粒组成、压强为10-4Pa的三维单分散颗粒体系，并计算了其结构因子，如图2插图所示.可以发现，三维单分散体系的静态结构因子在低k区域满足：S(k)∼k，与文献[1]中的结论一致.造成二维和三维体系S(k)在低k区域不同的原因可能是，在二维体系中更容易形成长程关联.

图3为不同分散度下二维颗粒体系S(k)的云图.可以看出，单分散体系(即s=0)表现出晶体的特征，其S(k)呈三角格子排列;而s=0.1的体系则表现出典型的多晶衍射图样的特征，衍射图形呈现出明暗不均匀的同心圆环;随着粒径多分散度的增大，图形中从中心开始的第二个圆环由正六边形演变为圆形，圆环的径向宽度也随之增加，而且看不到明显的点;当s≥0.3时，同心圆环变得更宽且彼此合并，其外围轮廓逐渐变得模糊，表明系统随着s的增大越来越无序;当s=0.5时，只能看到一个比较清晰的圆环.

3.3 取向序关联

为了量化颗粒i的取向序，引入键取向序参数Ψ6i[16]：

式中，ni表示颗粒i的最近邻颗粒的数目，θij表示颗粒i和其近邻颗粒j之间的极角.为了进一步量化取向序的空间关联，本文计算了颗粒体系的取向序关联函数g6(r)[3,17,18]，定义为

图3 不同分散度下静态结构因子S(k)

图4 取向序关联g6(r)随r的变化及拟合参数ξ6随分散度s的变化规律 (a)g6(r)随r的变化;(b)ξ6随s的变化;实线为ξ6∝14.2e-s/0.2

图4(a)给出了不同分散度下的取向序关联函数g6(r).由图4(a)可知，取向序关联函数呈现随着r的增大幅值逐渐减小的振荡行为.其中，g6(r)的极小可能与颗粒没有占据晶格的位置有关，而g6(r)振荡衰减的长度范围则对应取向关联的长度.为了得到取向序关联长度，我们对曲线进行了e指数拟合，如图4(a)所示.虚线为拟合曲线：g6(r)∝ae-r/ξ6，其中ξ6是序关联长度.由图4(a)可知，s＞0.1的g6(r)曲线均呈e指数衰减，表明s＞0.1的颗粒体系表现出与液体类似的短程序[3,19].图4(b)给出了序关联长度随分散度变化的曲线，显然，随着s的增加序关联长度近似以e指数规律衰减，如图4(b)中的实线所示.总之，随着粒径分散度s增大，ξ6减小，而且减小的幅度越来越慢，变化比较连续.说明随着粒径多分散度s从0.1到0.5，体系越来越无序，而且它的变化并不是一个突变，而是一个结构连续变化的过程.此时，体系无序程度对颗粒间取向序关联的影响比较明显.

3.4 力的累积分布

在静态颗粒排布中，颗粒间接触力形成高度各向异性的力网络，其基本特征可以用接触力的概率密度函数p(f)表征，其中f=F/〈F〉为归一化力.近年来，大量的实验和模拟关注大f处p(f)的分布，由于实验上的困难，对于小f处p(f)分布的研究还不是很多[20].我们计算了不同分散度下颗粒体系中的力累积分布G(f)，G(f)=图5可以看出，当f＜2时，随着f的增大，G(f)值也增大;当f＞2后，G(f)趋近于1.而且，不同分散度下的G(f)曲线几乎重合，暗示系统的无序程度对力累积分布几乎没有影响.

图5 力的累积分布G(f)随 f的变化

4 结论

采用颗粒离散元方法生成了具有不同尺寸分散度的二维颗粒体系，进一步研究了尺寸分散度对体系力学和几何结构特征的影响，得到如下结论：1)二维颗粒体系，平均配位数和平均剪切模量随着s的增大而减小;当s＞0.1时，二者都基本保持为定值，说明随着s的增大，体系的无序程度增加;2)当分散度s≤0.02时，静态结构因子低k部分的曲线基本重合;当s＞0.1时，随着s的增大，S(k)的值也均匀增加，尤其是s=0时，二维颗粒体系S(k)的低k区域行为与聚合物颗粒链的类似，S(k)∼k-1/υ(υ=3/4)，这一点不同于三维单分散颗粒体系(其S(k)曲线低k部分斜率近似等于1)，暗示二维单分散颗粒体系(包括分散度较小的体系)结构上与颗粒链相似;3)不同分散度下，取向序关联函数g6(r)曲线的峰值均满足e指数关系，序关联长度ξ6也随尺寸分散度的增大而减小;4)力的累积分布G(f)不随尺寸分散度的变化而变化，说明它基本不受系统无序程度的影响.

[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