王蓬 孔平 李然 华云松 厚美瑛 孙其诚
1) (上海理工大学医疗器械与食品学院, 上海 200093)
2) (上海健康医学院分子影像学重点实验室, 上海 201318)
3) (上海理工大学光电信息与计算机工程学院, 上海 215300)
4) (中国科学院物理研究所, 北京 100190)
5) (清华大学土木水利学院, 北京 100084)
研究颗粒体系中的结构与缺陷对于研究固-液融化的物理机制具有重要的价值.本文实验研究了垂直振动下单层湿颗粒在固-液融化过程中的结构与缺陷.根据实验及理论分析构建了湿颗粒体系的接触模型, 量化了准二维湿颗粒体系融化过程中颗粒的结构变化.然后以颗粒为点建立Voronoi图对颗粒体系的“相”转变进行研究, 并引入了局部体积分数来确定融化过程中缺陷变化的临界状态.实验结果表明, 颗粒系统在团簇的边缘开始发生缺陷, 并呈现链状的缺陷对向中心蔓延的现象.并且, 颗粒发生缺陷时七相缺陷颗粒的局部体积分数显著减小, 明显小于五相缺陷和六角相颗粒的局部体积分数.对局部体积分数的分析表明, 当最小局部体积分数 φ ≤0.6652 时发生缺陷, 当 φ ≤0.4872 时颗粒系统发生从固体到液体的转变.
研究湿颗粒在融化过程中结构与缺陷的演变对了解固体和液体之间的转化发挥重要作用[1,2].缺陷的产生与蔓延导致结构发生显著变化, 是固-液融化“相”变的物理机制[3,4].在准二维颗粒体系中, 粒子的结构可用来表征“相”状态[5].因此, 可根据粒子在空间中缺陷的变化来描述“相”变化的过程.虽然目前在模拟研究中和对胶体的实验研究中能够对准二维系统整体的结构变化进行定量研究[6,7]; 但很难从局部上对“相”变化过程中的缺陷进行定量分析.
湿颗粒体系中粒子之间存在具有相互吸引的内聚力[8], 颗粒的空间结构变化易于测量, 因此可用湿颗粒体系研究固-液融化“相”转变过程[9].Bossler[10]研究了湿颗粒毛细管力对颗粒簇结构的影响, 证明了湿颗粒的结构受毛细管力影响.Ramming等[11]研究了垂直振动下单层湿颗粒融化的集体行为, 通过键向序参数来表征湿颗粒的结构形式, 得到了聚类、重组和融化的状态, 为研究湿颗粒在聚集后的结构与缺陷提供了一定的理论基础.Kosterlitz和Thouless[12]建立了固-液相变理论, Nelson等[13]对该相变理论进行了扩展, 统称为KTHNY理论.KTHNY理论认为固-液相位的转变在二维体系中不是一阶转变, 而是颗粒系统发生缺陷并进行结构重组的过程[14,15].在湿颗粒系统中, 颗粒之间所特有的内聚力使得其内部结构和运动状态均比较复杂[16].离散元仿真研究表明: 在准二维约束条件下, 非弹性颗粒在融化过程中其颗粒“相”图呈现出有序固相和流态无序相的状态[17].胶体颗粒固-液融化的研究是通过对径向分布函数和方向关联函数的分析以确定相位演化的规律和结构作用的范围[18,19].目前, 对固-液融化过程中结构缺陷的研究主要是通过取向关联函数的指数衰减来表征系统整体的结构变化, 但如果想进一步研究“相”转变过程中的结构缺陷, 就有必要探索一种更有效的量化指标对固液融化过程中的局部的结构与缺陷进行定量分析.
本文对垂直振动下单层湿颗粒在固-液融化过程中颗粒的结构与缺陷进行了定量研究.首先, 介绍了实验系统及颗粒接触的物理模型, 利用图像处理技术对每个颗粒的位置进行精准定位; 然后, 以颗粒为点建立Voronoi图研究了颗粒系统的结构变化; 最后, 通过引入局部体积分数, 对缺陷的临界状态进行了定量研究.
实验装置如图1(a)所示, 垂直振动台采用电动式振动系统.该系统由扫频发生器(SA-SG030)产生的正弦信号驱动, 利用功率放大器(SA-PA080)来调整激振器(SA-JZ050)振动参数.采用高速相机对图像进行采集, 然后使用图像处理技术对采集的颗粒图像进行重构, 如图1(b)所示.正弦信号振动频率f、振幅r, 无量纲加速度Γ=4π2f2r/g, 其中,g为重力加速度.
实验中使用黑色玻璃珠, 颗粒参数如表1所示.将黑色玻璃珠与少量的超纯水 (LaboStar_TWF_7)混合后添加到一个以玻璃钢为材料的圆柱形容器中, 圆盘高度H=(2.5±0.05)mm , 内径D=13.00cm, 其中H的选择需要比d稍大一点,以确保颗粒系统保持单层.颗粒数N=2355 , 全局面积分数φ=Nd2/D2≈55.7%.颗粒的含水量被定义为W≡Vw/Vg[11], 其中Vw表示液体的体积,Vg表示颗粒的总体积.本实验中加入液体的含量大约为3%, 使相邻颗粒之间形成液桥[20].在实验过程中, 为了使液体的含量基本保持不变, 在圆盘一周圈添加密封圈, 减少水分的蒸发.
表1 颗粒参数Table 1.Particle parameters.
图1 实验装置示意图Fig.1.Schematic of experimental setup.
为了确保实验的可重复性, 在每次实验前, 首先将振动系统的加速度调整到Γ≈7 运行大约3 min, 使颗粒间的液体均匀分布, 提供一个稳定的初始聚集状态条件.为了防止颗粒与边壁的粘连, 颗粒初始状态尽可能地向中间聚集.随后, 逐渐增加Γ并对颗粒系统进行图像采集.在图像采集前, 调节振动台振动频率, 使相机采样与振动同步, 从保证连续两帧图像的颗粒像素大小不受振动的影响, 所采集的原始图像如图2(a)所示.
图2 (a)原始图像; (b)颗粒的定位Fig.2.(a) The original image; (b) particle localization.
图像处理方面, 为了实现对图片中白色光点准确的识别, 首先对相机采集到的图片进行去噪、二值化、连通域特征提取处理.然后, 以图像左下角为原点建立坐标系, 获得每个颗粒的位置信息, 如图2(b)所示.为了排除水面反光导致的定位错误,本研究采用相互位置查找算法来识别颗粒之间的距离[21], 进而判断该白色光点是否为颗粒.
在本研究中, 由于液体的存在使得相邻粒子之间形成的毛细桥, 从而发生近距离的相互吸引[22,23].为了测量邻近颗粒的结构, 需要对每个粒子进行键向序(bond-orientational order, BOO)参数Qn识别[22,24].图3为邻近颗粒的BOO参数识别示意图.
键向序参数的公式计算在湿颗粒系统中被广泛使用[11].其计算方法主要包括以下3个步骤.首先, 确定粒子i与其所有邻居之间的键.如:rij表示键从颗粒i指向j, 当颗粒i满足|rij|≤rc时(rc=1.05d), 确定颗粒j为颗粒i的邻居颗粒, 其中|rij|表 示两颗粒的键长,rc的选择通过颗粒系统的Voronoi图的邻域来确定.从方位角φij和极位夹角θij(θij=0 , 对于2D)的颗粒开始, 计算
其中,Ynm(θij,φij) 表示颗粒i与j的球面调谐[21].然后, 计算局部体积元上所有液桥的平均相关参数
其中,Ni为粒子i的键数.最后, 得到颗粒i在颗粒系统中的B00参数Qn,
通过BOO参数识别可以得到每个粒子的局部结构排列的相位图, 如图4所示.为了便于研究系统中颗粒结构的演变, 不同的结构标识为不同的颜色, 其中Q0表示该颗粒为自由颗粒, 用蓝色表征, 如图4(a)所示;Q2表示该颗粒与两个颗粒相接触, 用绿色表征, 如图4(b)所示;Q4表示该颗粒与4个颗粒相接触, 用黄色表征, 如图4(c)所示.Q6表示该颗粒与6个颗粒相接触, 结构呈六角相, 用红色表征.其中六角相在“相”位的变化中比较明显.与其他局部结构的测量方法(如配位数或局部面积分数)相比, 使用BOO参数识别的优点是最小化周围颗粒对颗粒簇边缘的影响, 这对于分析小簇的结构至关重要.
图3 颗粒的局部结构识别示意图Fig.3.Schematic diagram of local structure identification of particles.
图4 接触模型的结构形式图, 蓝色表示自由颗粒, 在颗粒系统中没有接触颗粒; 红色表示具有六个接触颗粒的六角相颗粒.Fig.4.Diagram of the structure of the contact model.The blue color represents free particles, and there are no particles in contact with the particle system.Hexagonal phase particles with six contact particles in red.
根据接触模型, 通过BOO参数识别来表征不同加速度下颗粒系统的结构变化.
在不同加速度下, 通过BOO参数识别所得到的单层湿颗粒在团聚后颗粒系统的结构变化情况,见图5.使用不同结构颗粒的数量占全局颗粒的数量的百分比ξ来衡量颗粒系统结构的变化, 红色和蓝色曲线分别表示六角相(hexagonal phase)和自由颗粒所占全局颗粒百分比随无量纲加速度Γ的变化.本文误差棒是以多次实验所测量的局部体积分数ξ的标准差为中心, 线段长度的一半表示不确定度.在较低的加速度下, 颗粒系统始终保持初始聚集状态, 颗粒的结构表现为长程的六角固相.在Γ=7时, 颗粒簇会在圆盘内无规律移动, 单个颗粒由于受到容器底壁垂直振动带来的能量的影响,颗粒在垂直面上有一个小的位移, 但是液桥力的作用使得相邻颗粒之间并不会发生相对移动.随着振动强度的增加, 颗粒获得的能量也增加, 大的颗粒簇逐渐解体重组为小的团簇, 导致六角相迅速降低.这与在晶体中观察到的晶体缺陷相似[25].由于颗粒所受到底壁的能量注入越来越大, 颗粒变得越来越活跃.当两个颗粒的相对分离距离大于毛细桥的破裂距离时, 单个颗粒所获得的能量则足以解除液桥力的束缚融化为离散的粒子.这种颗粒系统的集体行为在凝聚态物理学中被认为是一种从固态到液态的转变过程, 在转变的过程中颗粒系统会发生结构的变化[26,27].在平衡胶体颗粒体系[28]中所观察的融化行为与这种非平衡系统相类似, 这也表明了非平衡系统的融化转变与平衡系统之间必然存在某种联系.
在目前的研究中, 探索颗粒结构变化的方法大多都是通过分析颗粒的径向分布函数G(r) 和定向分数函数g6(r) 的衰减形式来确定“相”位的变化,但是要想研究融化过程中结构的缺陷还需要找到一个合适的变量来定量的对缺陷的演变进行分析.在实验中, 融化是通过振动幅值的增加使颗粒的平均动能增大所驱动, 这也暗示了颗粒温度(δv2)和热力学温度之间具有类似的特性.然而, 在平衡硬球体系的研究中发现颗粒温度并不是一个相关的变量, “相”转变完全由密度变化所驱动.其实, 在图5所示的融化过程中有一个类似的机制在起作用.随着振幅的增大, 单层的颗粒开始在有限的圆盘空间内探索所有可用的体积, 从而导致密度发生变化.因此, 我们可以寻找一个合适的变量来定量地研究这种密度的变化.
图5 不同加速度下单层湿颗粒在团聚后颗粒系统的结构变化图.颗粒的结构是采用BOO参数识别并颜色表征后的结果, 曲线图表示了随着加速度的增加, 六角相和自由颗粒所占全局颗粒的百分比 ξ 的变化.Fig.5.The structure change diagram of single layer wet particle system after agglomeration under different accelerations.The particle structure is the result of BOO parameter identification and color characterization.The curve shows the change of the proportion of hexagonal phase and free particles in the system.
局部体积分数是一种改进的Voronoi cell[29]法, 用来计算颗粒的局部密度.根据Voronoi图中胞体的几何形状与面积来定量地表征结构的变化,与其他方法相比, 该方法更接近于颗粒i在局部内可移动的空间.根据Voronoi cell可得到颗粒i的局部体积分数的定义[30]:
其中,Vi表示第i个颗粒的体积,表示该颗粒在维诺图中多面体的体积.实验中, 每个颗粒的局部体积分数可以直接根据平面图像中颗粒的像素和与所对应多边形的面积之比来定义.当振动加速度Γ=15 时, 每个颗粒在Voronoi cell中的分布见图6(a), 为了避免出现区域块边缘切割粒子所导致局部体积分数不准的问题, 在图像处理的过程中, 利用空间位置查找法排除区域边缘的颗粒.红色七边形表示颗粒的结构为7-fold缺陷, 蓝色五边形表示颗粒的结构为5-fold缺陷.可以发现7-fold颗粒和5-fold一起构成一个缺陷, 缺陷也可以相互结合形成一个5-7-5-7结构的缺陷对.
图6 (a), (b), (c)分别为在 Γ =14,14.4,15 时颗粒在Voronoi cell中的结构形式; (d), (e), (f) 的柱状图分别展示了图6(a), (b),(c)加速度下颗粒的局部体积分数, 红色柱子表示结构为7-fold颗粒的局部体积分数, 蓝色柱子表示结构为5-fold颗粒的局部体积分数, 黑色柱子为发生缺陷前结构为六角相颗粒的局部体积分数Fig.6.(a), (b), (c), respectively in when Γ = 14, 14.4, 15 particles in the Voronoi cell structure of the form; the histograms of (d),(e), (f) respectively show the local volume fraction of particles under the acceleration of Figs.6 (a), (b) and (c).The red column represents the local volume fraction of 7-fold particles, the blue column represents the local volume fraction of 5-fold particles, and the black column represents the local volume fraction of hexagonal phase particles before the occurrence of defects.
目前, 对这些结构缺陷的形式[18,19]研究的比较多, 但很少有对这种局部结构转变的临界状态进行定量的研究.图6(a)为颗粒系统发生局部结构缺陷前的Voronoi图, 其中颗粒的局部体积分数φ对应于图6(d)(6个柱子是随机选择的6个颗粒的局部体积分数); 此加速度下颗粒系统表现为类晶体状态, 在Voronoi图中的相图全部为六角相.图6(b)为颗粒系统开始发生局部结构缺陷时的Voronoi图, 其中红色Cell和蓝色Cell分别表示颗粒产生局部的7-fold缺陷和5-fold缺陷, 并且5-7缺陷对总是同时出现.图6 (e)是对2个缺陷对的颗粒局部体积分数进行测量得到的结果, 红色与蓝色柱子分别表示7-fold和5-fold颗粒的局部体积分数, 黑色柱子表示在缺陷发生前六角相的局部体积分数.有趣的是, 由图6(e)可以很明显看出, 发生结构缺陷后7-fold缺陷的颗粒局部体积分数显著变小, 并明显小于5-fold缺陷和‘六角相’颗粒的局部体积分数.根据Olafse等[18]的研究可知, 7-fold缺陷的出现会导致系统产生结构缺陷[18], 5-7缺陷对的解绑导致颗粒系统开始发生固-液转变[19].因此, 可以预测: 通过测量颗粒系统中最小局部体积分数(7-fold缺陷的颗粒局部体积分数)的变化来量化缺陷的演变.为了验证颗粒在缺陷阶段7-fold缺陷颗粒的局部体积分数始终保持最小, 又测量 在Γ=15 时的颗 粒 的 局 部体积分数柱状 图, 如图6(f)所示, 发现在更高的加速度下7-fold缺陷颗粒的局部体积分数依旧保持最小.所以, 可以通过测量颗粒系统中最小局部体积分数的变化来量化缺陷的演变.
局部体积分数这一变量的引入, 成功的把这种耗散的非平衡系统的集体行为进行了量化处理.局部体积分数不仅适用于准二维系统, 而且在三维体系中更能充分地对局部可移动的空间进行表征,也有望对这种耗散的系统进行非平衡统计力学的分析.
在不同的振动加速度下, 颗粒系统中最小的局部体积分数和平均局部体积分数的变化如图7所示.当Γ≤8 时, 由于振动强度较小, 颗粒系统几乎处于静止状态, 局部体积分数保持不变.当8<Γ<14, 颗粒系统中处于团簇边缘的颗粒开始发生微弱的移动, 但由于受到液桥力的影响, 颗粒系统始终保持稳定状态, 此时颗粒的最小局部体积分数逐渐减小, 但平均局部体积分数几乎不变.在Γ<14这个阶段颗粒系统的“相”图始终处于以长程定向顺序为特征的结晶相.
逐渐增加Γ, 从Voronoi cell中可以看到, 在Voronoi cell右边缘位置(团簇边缘)颗粒的结构开始发生5-7 fold缺陷, 表明在Γ=14 时颗粒系统开始发生结构缺陷.此时颗粒系统的最小的局部体积分数φ=0.6652 , 并且颗粒系统的平均局部体积分数急剧下降.这是因为团簇边缘的颗粒受到的液桥力较小, 颗粒由于与底壁的碰撞所带来的动能使颗粒变得活跃起来, 开始探索可用的空间.这种现象与在旋转盘[2]中所得到的湿颗粒物质从表面开始发生融化的结果相一致.在 1 4<Γ<36 这个阶段, 颗粒的局部体积分数先显著降低, 然后达到稳定的状态, 衰减趋势呈指数型.此阶段颗粒系统的缺陷从团簇边缘的位置逐渐向中间蔓延, 有趣的是缺陷更多呈现链状的5-7-5-7结构的链状缺陷对.究其原因, 是由于液桥力而产生的力链所导致的,使得颗粒更趋向于成排一起移动.但是, 从Voronoi cell中并没有观察到缺陷解绑形成独立的缺陷, 整个系统还是处于比较稳定的“类晶体”结构.当Γ=36时, 颗粒系统的缺陷急剧增加, 并出现7-fold或5-fold孤立的缺陷, 这说明颗粒系统开始发生从固体到液体转化.在更高的加速度下, 整个颗粒系统结构杂乱无章, 系统逐渐被液化, 这也与图5中所分析的结果一致.
图7 不同振动强度下颗粒系统中局部体积分数的变化曲线.其中红色曲线表示系统中最小的局部体积分数变化, 蓝色表示系统中所有颗粒的平均局部体积分数变化.图中Voronoi cell表示在该加速度下颗粒系统的“相”位图, 其中Voronoi cell右边缘处于颗粒系统中团簇的边缘位置Fig.7.The change curve of local volume fraction in particle system under different vibration intensity.The red curve represents the minimum local volume fraction change in the system, while the blue curve represents the average local volume fraction change of all particles in the system.Voronoi cell in the figure represents the phase diagram of the particle system under this acceleration, where the right edge of Voronoi cell is located at the edge of the cluster in the particle system.
为了验证临界值是否具有迟滞效应, 分别对加速度递增和加速度递减情况下最小局部体积分数φ的增减变化趋势进行了研究, 如图8所示.在加速度Γ从 0至44 逐渐增加时, 最小局部体积分数的衰减趋势与图7 相吻合; 但是在加速度Γ从 44至0 逐渐减小时, 最小局部体积分数的增长曲线在图8两虚线之间的这个区间具有很大的不确定性, 与加速度递增时最小局部体积分数的增减趋势并不相同, 所以此临界值不具有迟滞效应.
图8 最小局部体积分数 φ 的增减变化趋势, 红色曲线表示加速度 Γ 从0依次增加到44时颗粒系统中最小局部体积分数的变化, 黑色曲线表示加速度 Γ 从44依次减小到0时颗粒系统中最小局部体积分数的变化Fig.8.The trend of the increase and decrease of the minimum local volume fraction φ.The red curve shows the change of the minimum local volume fraction in the particle system as the acceleration gradually increases from 0 to 44, and the black curve shows the change of the minimum local volume fraction in the particle system as the acceleration decreases from 44 to 0.
从以上分析可知, 颗粒系统的结构与缺陷是从团簇边缘位置开始发生且结构变化存在一个临界值, 即 当φ≤0.6652 时 颗 粒 发 生 结 构 的 缺 陷, 当φ≤0.4872时颗粒系统发生从固体到液体的转变.该临界值的误差是通过误差棒来量化的.误差棒的测量是在相同的初始状态下重复30次得到的实验数据, 然后以局部体积分数的标准差为中心, 其中误差棒一半的长度表示不确定度, 只要临界值的波动在误差棒的范围内, 该临界值就可以接受.在准二维湿颗粒体系中, 颗粒可以被看作分子, 用相邻颗粒之间的液桥代替分子键, 通过施加的外力来控制系统的融化状态.外部驱动力的增加使“晶体”克服摩擦而相对移动, 由于“晶体”表面粒子间的相互作用较弱, 导致“晶体”表面颗粒变得稀疏而发生结构的缺陷, 并向“晶体”内部结构松散区域扩散.在缺陷阶段, 由于外部能量的注入与液体桥的破裂能量达到平衡, 因而分子解除键的束缚, 最终导致“晶体”产生向液体融化的突变[23].这与图7中φ=0.4872时, 从缺陷到开始发生固-液转变的临界状态一致.然而, 对于在缺陷扩散时所出现的链状和环状的结构的物理机制尚不清楚.
综上所述, 在今后对二维融化系统的研究中,可以打破仅使用序参量(如取向关联函数)的情况,利用局部体积分数来量化结构的演变.目前对二维融化的集体行为和结构演变的量化较为成熟, 但是对于结构缺陷的物理机制的研究仍具有挑战性, 所以今后的研究重点应放在对相转变过程中的dislocationvs.disorientation进行研究, 定量研究结构缺陷的发生与转化的物理机制.
本文对垂直振动下单层湿颗粒在聚集后发生的结构与缺陷进行了定量研究.根据实验及理论分析构建了湿颗粒体系的接触模型, 量化了准二维湿颗粒体系融化过程中颗粒的结构.实验结果表明,缺陷是导致固-液“相”转变的重要因素.缺陷在团簇的边缘开始发生, 并呈现链状的5-7-5-7结构的缺陷对向中间蔓延.并且, 颗粒发生缺陷时7-fold缺陷颗粒的局部体积分数会显著的减小, 并远远小于其他类型结构的颗粒.对局部体积分数的分析表明: 颗粒系统在最小局部体积分数φ≤0.6652时发生结构的缺陷, 当φ≤0.4872 时颗粒系统发生从固体到液体的转变.
感谢昆山杜克大学黄锴老师的细心指导与帮助, 感谢上海理工大学光电信息与计算机工程学院杨晖老师提供的实验设备与支持.