李歌天,曾援,张波,张国华

(北京科技大学数理学院,北京 100083)

在自然界中,存在着大量的颗粒物质,它有着与气体、固体和液体所不同的特点和性质,颗粒物质作为一种复杂的离散体系,无序性是它的重要特点[1]。想用现今的理论知识对颗粒物质进行完整描述是非常困难的,这也是研究颗粒体系的意义所在。理论表明,晶体材料的低频振动为长波平面波,其低频态密度(density of states,DOS)满足D(ω)～ωd,其中,d为维度,D(ω)为低频态密度,ω为振动圆频率。但是非晶材料的低频态密度出现大量偏离德拜预测的额外振动模式,这些额外模式对非晶材料的力学性质有非常大的影响。因此,研究无序材料振动态密度对其力学行为的影响,一直是材料力学的研究热点之一。近年来,许多学者对颗粒体系进行了试验和数值模拟,以探究它的振动特性、堵塞、耗散和颗粒链的尺寸效应[2-5]。如颗粒体系的振动态密度、态密度图像、堵塞相图以及临界点J处的特征等。文献[6-8]对颗粒物质的振动态密度进行了研究,将Hessian矩阵的简正模式分析法用于研究无序颗粒体系(如液体、玻璃介质等)的振动态密度。Zhang等[9]研究了无序程度对玻璃状固体振动态密度的影响,发现当二维颗粒填料从晶体结构逐渐调整为多晶结构乃至非晶结构时,玻色峰和van-Hove singularity奇异性在多晶体中很好地分离,且当无序性足够高时van-Hove singularity奇异性消失。Krishnan等[10]用Hessian矩阵法研究了剪切无序固体的非德拜低频振动,发现构型选择对于低频态密度的重要影响,并且发现当接近塑性事件时,Hessian矩阵的最小本征值到达一个普适的分布。Wang等[11]研究了三维玻璃材料的低频态密度,发现当频率低于第一声模式时三维玻璃材料的低频态密度D(ω)～ωβ,β<4,β为幂指数,对于不同体系取不同的数值。Dale等[12]研究了超均匀(hyperuniform)三维球形颗粒填料的振动态密度,发现超均匀球形颗粒填料在低频下不存在态密度随频率的标度,而是在有限频率下的突然低频截止,形成了大约零频率的声子带隙。在这个间隙上方,存在一个离散的模式带。然而,目前关于不同维度和压强下颗粒材料低频态密度的系统研究尚鲜见报道。

为此,运用离散元法(discrete element method,DEM),借助PFC软件生成不同维度下,压强P处于0.1～10 000 Pa的范围内,由2 048个无摩擦颗粒组成的颗粒体系数值模型,再运用MATLAB软件,通过Hessian矩阵的简正模式分析法,对数值模型进行计算并总结其规律。为进一步理解维度、压强对颗粒材料的低频力学性质的影响有帮助。

1 颗粒物质

颗粒物质指的是离散且半径大于0.5 μm,间隙液体的黏性较低,并且饱和度在1以下的大量固体颗粒,在相互作用下组成的复杂体系。它广泛存在于自然界中,在日常生活、生产等各个方面中都有所涉及[13]。

颗粒物质的性质及其丰富[14],如在没有外力作用时,颗粒体系的状态与固体的状态相似,有外力作用时,颗粒物质会发生流动,此时它的状态会更接近气体和流体[15]。该体系属于宏观,但又和微观相联系,属于固体,但又和流体和气体相关,具有较强的综合性,对它的研究可以涉及生活的各个方面。作为一种具有代表性的无序系统,颗粒体系的振动态密度与非晶体材料具有类似的特点,比如说存在许多低频软模式[16]、Boson峰和局域化模式等[17]。

2 数值模拟

2.1 样品制备

用PFC2D分别生成包括2 048个光滑圆盘颗粒的一维颗粒链和处于2 m×2 m方盒子中的二维颗粒圆盘填料,采用周期性边界条件并忽略重力的影响。采用PFC3D软件在2 m×2 m×2 m立方盒子中生成了包括2 048个光滑球形颗粒的三维颗粒填料。

具体颗粒填料的产生方式如下:首先,在选定的边界内随机产生2 048个半径很小的颗粒(此时颗粒之间重叠量很小,可认为处于松弛状态)。然后,增大颗粒半径,使得体系压强达到2×104Pa,采用共轭梯度法使得体系处于能量极小的状态,生成初始压强的颗粒填料。为了得到不同压强下的颗粒填料,通过逐渐减小颗粒半径的方法进行卸载,使得体系压强分别达到10 000、1 000、100、10、1、0.1 Pa。填料稳定的判据是经过1 000个循环前后系统的能量差的比例小于1×10-15填料稳定后,将各个颗粒填料的信息保存到文件。对每个压强都随机生成了10个位形,所涉及的物理量都是这10个位形的平均值。

2.2 态密度的求解

体系振动模式及振动态密度可通过求解下列久期方程来实现。

|F-ω2I|=0

(1)

式(1)中:I为单位矩阵;ω为振动圆频率;F为Hessian矩阵。

(2)

光滑圆盘颗粒之间的相互作用为单边线性弹簧势,一般情况下,颗粒i和j的相互作用势可表示为

(3)

“区块链”已成为信息技术应用行业的热门词，区块链技术应用也得到了国家在信息技术领域的关注和重视，一方面加大对首次币发行（Initial Coin Offering，ICO）项目的监管，另一方面积极推动国内区块链在相关领域的标准化制定及产业化发展等研究工作［1］。从“互联网+”逐步上升到“区块链+”。无论从社会价值还是投资价值投资考虑，在这股热潮之中，最终需要沉淀和落地的是“区块链+实际场景”的应用，只有先进的技术应用到具体项目才能拥有可持续的生命力，才能创造出真正的价值［2］。

(4)

式(4)中:ki和kj分别为颗粒i和j的法向接触刚度系数。

光滑球形颗粒i和j的相互作用为Hertz势,可表示为

(5)

(6)

对于各向同性材料,E=2G(1+ν),其中G为剪切模量。

3 结果与讨论

3.1 一维颗粒链的振动态密度

在一维条件下,选取若干个一维颗粒链,并对所有结果进行了平均计算,运用Hessian矩阵求解稳定状态下harmonic接触的一维颗粒链的振动态密度,并绘制态密度D(ω)～ω曲线,如图1所示。

图1 Harmonic接触的一维颗粒链态密度图Fig.1 Density of states curve of one-dimensional particle chain with harmonic contact

从图1可以看出,在不同压强条件下,由2 048个颗粒组成的一维无摩擦颗粒链体系的振动态密度完全一致。一维颗粒体系的振动态密度在0.1～10 000 Pa的压强下时,其振动态密度曲线不受压强影响。且曲线在ω的值较小时,曲线振荡频率较小且较为稳定,在ω逐渐增大后,曲线振荡幅度变大,且D(ω)的平均值随着的增大呈逐渐上升趋势,将ω的值和D(ω)都用对数坐标进行表示时,更能体现这一特点。

如图2所示,不同压强的曲线完全重合,即压强变化与态密度曲线的变化无关,并且在低频区域,态密度曲线趋近于一条直线。根据德拜理论,D(ω)～ωγ,γ=d-1。其中,d为体系维度,当d=1时,γ=d-1=0。因此,在一维条件下,D(ω)～ω0,反映在曲线上即在低频区域,态密度曲线趋近于直线。可以看出,图2的曲线符合德拜理论的预测。

图2 双对数坐标下的harmonic接触的一维颗粒链态密度图Fig.2 One-dimensional particle chain of the density of states for harmonic contact in double logarithmic coordinates

3.2 二维颗粒体系的振动态密度

首先用PFC2D软件在一个2 m×2 m方盒子中生成颗粒半径在0.009～0.011 m平坦分布的2 048个光滑二维颗粒填料。x、y方向均采用周期性边界条件。此时体系中颗粒所处位置完全是随机的,且颗粒间没有任何重叠,体系处于完全松弛的状态。接着,使颗粒半径逐渐增大,直至体系的压强到达一个固定值(P=1 000 Pa)。然后采用共轭梯度法使颗粒体系卸载到实验所需要的一系列压强值。随机改变初始样品中颗粒的位置,重复上述过程,从而制备若干压强相同二维颗粒填料。通过Hessian矩阵的简正模式分析法得到该压强下二维颗粒填料的平均振动态密度。图3为计算得到的压强分别为0.1、1、10、100、1 000、10 000 Pa的二维颗粒体系的态密度曲线。

图3 二维无摩擦颗粒体系的振动态密度Fig.3 Vibronic density of states a two-dimensional frictionless particle system

图4 双对数坐标下二维颗粒体系的态密度Fig.4 Vibronic density of states a two-dimensional frictionless particle system in double logarithmic coordinates

图5 双对数坐标下二维颗粒体系过渡频率随压强P的变化关系Fig.5 Variation of the two-dimensional system with pressure P in double logarithmic coordinates

3.3 三维颗粒体系的振动态密度

首先用PFC3D软件在2 m×2 m×2 m立方体盒子中生成颗粒半径在0.009～0.011 m平坦分布的2 048个光滑三维球形颗粒填料。x、y、z方向均采用周期性边界条件。接着,使颗粒半径逐渐增大,直至体系的压强到达一个固定值(P=1 000 Pa)。然后采用共轭梯度法使颗粒体系卸载到实验所需要的一系列压强值。随机改变初始样品中颗粒的位置,重复上述过程,从而制备若干压强相同三维颗粒填料。通过Hessian矩阵的简正模式分析法得到该压强下的平均振动态密度。图5为计算得到的压强分别为0.1、1、10、100、1 000、10 000 Pa的三维颗粒体系的态密度曲线。

由图6可知,对于三维无摩擦颗粒体系,每个一压强下的振动态密度曲线都呈现一个明显的峰值,D(ω)先随着ω的增大而增大,在到达峰值后随着ω的增大而减小;而在不同压强条件下,压强越大,其态密度曲线相比于压强较小条件下的态密度曲线,显得更为“矮胖”,曲线底部的宽度增大,D(ω)的最高峰峰值越小,且峰值越靠右,即在取到峰值时,ω的值越大。

图6 双对数坐标下三维颗粒体系的态密度Fig.6 Vibronic density of states a three-dimensional frictionless particle system

图7 双对数坐标下三维颗粒体系的态密度Fig.7 Vibronic density of states a three-dimensional frictionless particle system in double logarithmic coordinate

图8 双对数坐标下三维颗粒体系随压强P的变化关系Fig.8 Variation of the three-dimensional system with pressure P in double logarithmic coordinates

4 结论

采用离散元法生成了不同维度和压强下的颗粒体系,通过Hessian矩阵简正模式分析法计算得到颗粒体系的态密度,进一步研究了维度和压强对颗粒体系振动态密度的影响,得出以下结论。

(1)一维条件下,一维无摩擦颗粒链的态密度曲线不随压强的变化。在低频处,态密度曲线几乎不随频率变化;在高频率处,态密度曲线随着频率增大呈现振荡上升趋势。