致密硬球堆积是研究无序固体、颗粒物质以及无热介质中相变与阻塞(jamming)行为的基础模型体系之一^[1–6]. 硬球(hard sphere, HS)模型能够在最小相互作用模型中描述无序堆积、玻璃化以及结晶等典型结构现象, 因此长期被视为研究无序固体与相变问题的重要参照体系^[7–9]. 对于由宏观颗粒组成的无热体系, 堆积结构不仅决定力学稳定性与流动响应, 还在外部驱动(振动、剪切等)作用下呈现出丰富的结构和动力学演化过程^[10–12]. 然而, 在无热且具有摩擦的颗粒体系中, 无序球形颗粒在致密化过程中实现结构重排与对称性选择、并最终形成晶态有序结构的机制尚不清楚; 尤其是微观结构表征与统计力学描述之间的对应关系仍缺乏明确认识^[13,14], 理清这些问题是理解剪切驱动颗粒结晶的关键^[7,15,16].

围绕无序-有序转变的结构表征方面, 已有研究提出了多种可用于识别晶态生长起始点与增长的结构判据. 以键取向序(bond-orientational order, BOO)为代表的局域取向分析能够区分类液态与类晶态构型, 并追踪fcc与hcp类团簇的形成与演化, 揭示在体积分数跨越随机密堆积密度附近时无序与晶态区域的共存特征^[7,17]. 此外, 基于Voronoi划分的局域体积分布可揭示无序态、共存态与晶态区间在几何结构上的差异: 无序态中体积分布呈现宽广的类高斯分布形式, 而在共存区间出现双峰结构, 并在高密度晶态中收缩为窄峰分布, 反映出局域体积统计在转变过程中的系统变化^[13]. 另一方面, 结合高分辨率CT成像技术与配位数分析的研究表明, 在无序区、共存区以及晶态区间内, 几何配位数与力学配位数可呈现不同分支; 即使在近乎完美的晶体结构中, 接触网络仍可能保持显著的空间非均匀性^[12]. 这些结果表明: 无序到有序的转变不仅体现在序参量的变化, 也对应着体积统计与接触拓扑结构的协同重构.

除了几何结构表征之外, 近年来研究者们也尝试在Edwards体积系综的框架下, 为致密颗粒堆积建立统计力学描述. 在该理论中, 体系体积在形式上对应能量, 而其共轭的“温度”变量为可压缩率(compactivity)^[5,18,19]. 借助Voronoi体积涨落等可观测量, 可以在阻塞态颗粒体系中构建类似于自由能的热力学量, 从而为不同制备路径下的堆积状态提供统一表征. 已有针对振动与剪切驱动堆积的实验研究表明, 该方法能够在一定程度上一致地描述无序阻塞态所呈现出的临界特征^[20]. 然而, 对于三维硬球颗粒系统, Edwards体积系综能否合理刻画剪切诱导结晶过程中由无序到有序转变, 尤其是可压缩率与Edwards描述是否能够同时适用于阻塞态与晶态区间, 目前仍缺乏系统检验^[19].

基于上述背景, 本文聚焦于循环剪切诱导结晶的过程中, 三维颗粒体系的结构演化与统计性质. 通过键取向序分析方法与高分辨率CT成像技术相结合, 系统表征无序-有序转变过程中的结构演化. 结果显示, 当体系体积分数接近随机密堆积附近时, 结构序参量出现显著变化, 对应晶态区域的形成与增长; 在Edwards 统计框架下, 以体积涨落为基础得到的统计量能够对堆积的状态演化进行有效描述. 上述结果为Edwards体积系综在无热体系中的适用性提供了实验支持, 并为理解剪切驱动下无序到有序转变的结构与统计机制提供进一步依据, 也为相关离散元模拟(DEM)与理论研究提供直接的实验参照.

本文采用医用牙科CT(UEG Medical Group Ltd., 空间分辨率为0.2 mm)对循环剪切作用下颗粒体系的结构演化进行实验表征. 实验颗粒为约1万颗光滑 ABS(丙烯腈–丁二烯–苯乙烯)塑料球, 颗粒直径d ≈ 6 mm, 其多分散效应约为0.2%(约为0.012 mm). 颗粒装填于内部尺寸为$ 20 d\times 20 d\times 20 d $的有机玻璃剪切盒中. 为减弱重力导致的密度梯度, 剪切盒顶部放置一块可自由竖直移动的不锈钢盖板, 其质量约为颗粒总质量的1.5倍. 剪切加载通过底板步进电机驱动的直线位移台实现.

实验开始时将颗粒直接倒入剪切盒中形成可重复的疏松堆积. 随后施加循环剪切直至体系达到稳态, 其中前120次循环用于预处理. 循环剪切应变幅度$ \gamma =15{\text{%}} $, 剪切速率$ \dot \gamma = 0.15 {\mkern 1 mu} {{\text{ s}}^{{{- 1}}}} $, 对应惯性数$ I=\dot{\gamma }d/\sqrt{p/\rho } $= 5.79×10^–4, 其中p为顶部不锈钢盖板在重力作用下施加于颗粒堆积的平均法向压强, ρ为ABS颗粒的材料密度. 表明体系处于准静态剪切区间^[21]. 正式实验共施加4000次剪切循环保证结构充分演化并达到稳态. 为获得结构演化信息, 在剪切过程中按固定循环间隔采集三维CT图像, 并统一在$ \gamma =0 $(回到零应变)状态下进行扫描, 从而获得对应的阻塞构型.

此外, 为获得无序堆积的生成, 在完成上述预处理后我们对体系进行等步长分步剪切: 每一步将剪切应变增大$ \Delta \gamma =1.25{\text{%}} $, 直至累计剪切应变达到$ {\gamma }_{\text{max}}​=33.75{\text{%}} $. 我们在每完成一次应变加载后均进行一次CT扫描, 以记录颗粒构型随剪切应变的演化.

CT数据的颗粒识别与三维重构在Matlab中完成^[22]. 具体而言, 采用基于标记的分水岭分割算法提取颗粒的位置与尺寸信息, 其识别不确定度小于粒径的0.3%. 实验共进行3次独立重复. 为降低边界效应对统计量的影响, 我们在后续分析中剔除距容器壁小于3倍粒径范围内的颗粒, 仅对内部颗粒进行分析. 关于本研究选取3d (d为颗粒直径)作为体相(bulk)截取标准的考量如下: 首先, 这是一个关于减弱边界影响与保证统计样本量之间的权衡. 虽然增大剔除厚度(如5d)能更彻底地消除壁面约束, 但由于本研究样本尺寸有限, 过大的剔除厚度会导致内部统计样本量急剧减少, 从而引入较大的统计涨落. 其次, 我们进行了系统的敏感性分析(见图1). 结果表明, 当剔除厚度大于2d, 体积分数$ \phi $的演化趋势已趋于收敛. 选取3d不仅能有效规避明显的壁面有序化效应, 同时也保留了足够的颗粒数量以确保统计分析的鲁棒性. 体系总体体积分数$ \phi $定义为$ \phi = {\displaystyle\sum\nolimits_{i}{v}_{i, g}}\bigg/\bigg({\displaystyle\sum \nolimits_{i}v_{i}^{\text{vor}}}\bigg)$, 其中$ {v}_{i, g} $与$ v_{i}^{\text{vor}} $分别表示第$ i $个颗粒体积及其对应的Voronoi元胞体积.

图2(a)为循环剪切实验装置示意图. 图2(b)给出了体积分数$ \phi $随剪切循环次数t的演化. 体系初始体积分数约为$ \phi \approx 0.636 $, 随后随t增大逐渐升高并趋于平台. 不同实验在稳态平台值上存在轻微差异, 可能与初始微结构、动力学路径、边界以及有限尺寸效应有关; 总体而言, 稳态体积分数集中在$ \phi \approx 0.70 $附近. 为提高统计可靠性, 将3次独立实验的构型合并为统计样本, 并划分多个小体积分数区间对相关物理量进行平均, 给出其均值与标准差.

为定量识别剪切诱导结晶过程中由无序到有序的结构转变, 采用键取向序(bond-orientational order, BOO)方法对局域结构有序性进行表征. 对于每个中心颗粒$ i $, 选取其12个最近邻颗粒进行结构表征. BOO的思想是将中心颗粒与其邻居连线的方向分布映射到球谐基底上, 从而刻画局域排布相对于理想晶体参考结构的角向相似性与对称性特征^[16,23]. 相应地, 定义复数序参量:

其中$ {\hat{r}}_{ij} $为中心颗粒$ i $指向其近邻颗粒j的单位矢量, $ {\text{Y}}_{lm} $为球谐函数, $ {N}_{\text{b}}(i) $为颗粒$ i $的近邻数^[24]. 为消除坐标旋转对结果的影响, 进一步引入旋转不变量:

在三维硬球堆积中, 主要计算$ {q}_{4} $, $ {q}_{6}​ $以及三阶旋转不变量$ {w}_{6} $等 BOO指标, 用于描述结构随剪切演化的变化. 图3(a)给出了平均序参量$ \left\langle {q}_{6}\right\rangle $随总体体积分数$ \phi $的演化趋势. 结果表明, 体系可划分为两个具有不同结构演化的区间: 当$ \phi \lt 0.636 $时, $ \left\langle {q}_{6}\right\rangle $随密度缓慢增大, 说明无序堆积的局域有序性主要通过几何填充与局部重排逐步增强; 当$ \phi \gt 0.636 $时, $ \left\langle {q}_{6}\right\rangle $的增长明显加快, 并与晶态区域的出现相对应, 表明体系开始发生更具协同性的结构重排, 有序域迅速扩展. $ {\phi }_{\text{RCP}} $附近的斜率变化可理解为: 在无序致密化阶段, 体系基于几何填充与排斥约束下系统进行结构重排; 而体积分数在接近与超过随机密堆积(RCP)后, 在剪切驱动下体系更易通过协同重排进入低缺陷的晶态几何构型, 使BOO有序性出现更快的增长. 这与DEM 研究中基本趋势是一致的, 但实验中观察到的$ \left\langle {q}_{6}\right\rangle $转变点位置, 并不是在RCP的位置, 这个可能与外部驱动和表面摩擦相关.

为提高对结构转变测量的灵敏度, 进一步引入三阶旋转不变量$ {w}_{l}​ $(本文取$ l=6 $)作为有序性的补充表征^[17]:

其中$ (\cdot ) $是Wigner 3j符号. 相较于二阶不变量, $ {w}_{6} $对局域结构的相位关联更敏感, 因而能在$ \left\langle {q}_{6}\right\rangle $尚未显著变化时揭示对称性相关的微弱结构重排.

我们进一步定义归一化累积分布$ P({w}_{6}) $, 并将满足$ P(w_{6}^{\text{hh}})=1/2 $的半高位置$ w_{6}^{\text{hh}} $作为表征量. 图3(b)表明, 随着体积分数增大, $ w_{6}^{\text{hh}} $在$ \phi \approx 0.636 $附近出现明显跃升, 这与晶态区域的出现相一致, 表明体系在该体积分数附近发生结构重组, 并由无序向有序结构转变, 随后进入晶核快速生长阶段. 同样发现这与DEM 研究中基本趋势是一致的, 但实验中晶核出现的位置小于RCP的位置.

为识别类晶(crystal-like/solid-like)团簇, 采用ten Wolde-Frenkel提出的晶态键判据^[17]. 具体做法是: 先基于BOO向量定义相邻颗粒$ i $与j的取向相关性:

并将满足$ {S}_{ij} \gt 0.5 $的相邻对定义为一条晶态键^[17]. 随后, 对每个颗粒统计其晶态键数$ {n}_{\text{b}} $; 当$ {n}_{\text{b}}\geqslant b $(本文取b = 6)时, 将该颗粒判定为类固态(solid-like). 阈值的选取遵循经典成核研究中的“相态分布分离”原则^[17]: 通过比较纯无序参考态与晶态参考结构中$ {n}_{\text{b}} $的概率分布(图4(a)), 选择能够最大程度区分两相、同时降低无序相中误判的阈值, 从而提高类晶识别的稳健性.

图4(b)给出了类晶体颗粒分数随体积分数$ \phi $的演化. 结果表明体系经历了由无序向晶态逐步转变的过程: 在初始阶段($ \phi \lt 0.636 $), 类晶分数接近零, 体系主要保持无序几何构型; 随着体积分数增大, 体系进入无序与晶体共存的转变区, 类晶分数持续上升, 表明晶态团簇逐渐形成并扩展. 当$ \phi \gtrsim 0.69 $时, 类晶分数进入平台, 此时体系整体已经以晶态结构为主, 而进一步致密化主要表现为晶态内部缺陷的缓慢消除而非类晶体积分数的显著增长. 平台行为可理解为剪切诱导有序化与动力学约束之间的竞争平衡: 部分区域由于局域几何不匹配或重排受限, 难以完全满足晶态对称性要求, 从而在稳态中保留一定比例的无序成分.

在Edwards及其合作者提出的统计框架中, 无热颗粒体系的“阻塞态”(jammed microstates)可按体积进行计数^[18]. 将体积为$ V $的阻塞构型数记为$ {\varOmega }(V) $, 并定义Edwards熵为

在该框架下, 体积在形式上对应平衡统计中的能量, 而与体积共轭的“温度型”状态变量为可压缩率(compactivity) $ \chi $^[25], 其定义为

这一思想为堵塞态颗粒体系提供了一种以体积为核心的统计描述途径: 不同制备方法若能在$ \chi $上实现一致表征, 则有望建立类似热力学的状态变量体系. 实验上, $ \chi $通常无法通过直接计数$ {\varOmega }(V) $得到, 而需要借助体积(或体积分数)涨落来间接测量. 对于基于Voronoi 剖分得到每个颗粒对应的元胞体积$ {V}_{i} $, 定义其涨落为

在大量实验与模拟中, Voronoi(或radical Voronoi)体积涨落已被证明是表征堵塞态结构统计的重要观测量, 并可与Edwards体积系综的构造相衔接.

类比平衡统计中能量涨落与比热的关系, Edwards体积系综下可得到可压缩率与体积涨落之间的积分关系(涨落法):

该方法在颗粒体系的振动和剪切实验中被广泛用于从体积涨落反推出可压缩率$ \chi $. 其中, ϕ^r ≈ 0.60 (随机疏松堆积)为参考态的体积分数, χ = χ(ϕ^r)为该参考态对应的可压缩率. 在选取的随机疏松堆积参考态下, 取$\chi^{\rm r} \to \infty $, 即无穷温.

图5(a)展示了$ \sigma _{\text{V}}^{2} $随体积分数$ \phi $的演化. 对无序分支而言, $ \sigma _{\text{V}}^{2} $随$ \phi $增大总体呈下降趋势, 这反映了致密化过程中自由体积逐步收缩、构型空间被压缩的几何效应. 进入结晶共存区后, 体系中无序与晶态局域结构在空间上并存, 使得局域体积分布表现出明显的“混合”特征并被显著拉宽, 从而导致表观涨落$ \sigma _{\text{V}}^{2} $增强并在某一体积分数(实验观察在$ \phi \approx 0.68 $)附近出现峰值; 当晶态团簇进一步占主导且缺陷逐步减少时, 体积分布重新收敛, $ \sigma _{\text{V}}^{2} $再次下降并趋于稳定.

需要强调的是: 共存区涨落的增强很大一部分来源于“两类结构的混合导致的分布变宽”, 它未必等同于单一相内在的热力学涨落. 为避免将“共存混合”误当作单相本征涨落, 本文在计算$ \sigma _{2}^{2}(\phi ) $与后续$ \chi (\phi ) $时, 分别对无序颗粒与类晶颗粒(按前文BOO与类晶判据分类)进行条件统计, 并在各自分支内进行拟合, 从而得到两条分支对应的$ \sigma _{2}^{2}(\phi ) $.

针对无序分支, $ \sigma _{1}^{2}(\phi ) $可用低阶多项式拟合以获得平滑函数$ \sigma _{1}^{2}(\phi )=-28.656{\phi }^{3}+53.956{\phi }^{2}- 33.963\phi +7.152 $. 对于$ \phi \gt 0.68 $的区间, Jin和Makse^[19]指出当$ \phi \gt {\phi }_{m} $时共存态在热力学意义上并不稳定, 因此本文仅选取靠近晶态稳态($ \phi \approx 0.70 $)的高密数据拟合晶态分支, 从而得到晶态分支对应的$ \sigma _{2}^{2}(\phi ) $: $ \sigma _{2}^{2}(\phi )=20.547{\phi }^{3}-42.189{\phi }^{2}+28.709\phi - 6.467 $.

由(8)式可见为确定$ \chi $, 必须指定一个参考态 ($ \chi ({\phi }_{0}) $). 对无序分支, 通常取随机疏松堆积(random loose packing, RLP)作为参考态^[26]. RLP常被视为给定体系可实现的最疏松力学稳定堆积, 并在Edwards观点下对应“无穷可压缩率”极限$ {\chi }_{\text{RLP}}\rightarrow \mathrm{\infty } $, 即认为在该极限下阻塞构型等概率^[5]. 近期研究表明RLP并非普适的状态, 而与体系制备过程和颗粒摩擦直接相关^[26,27]. 在本文剪切实验中, 取实验中可稳定获得的最松无序稳态$ {\phi }_{0}=0.60 $作为参考 RLP 态, 并据此得到无序分支的$ \chi (\phi ) $(图5(b)左支), 其结果与我们过去的工作一致^[20].

晶体态分支对应参考态的选择更为微妙. 原则上, 完美晶体可作为Edwards 熵为零的参考态; 但实验上难以获得任意接近完美晶体的状态, 并且若直接以“完美晶体”为$ \chi $的参考, 积分在实际数据区间往往对积分上下限假设高度敏感. 在Jin和Makse^[19]的工作中, 他们认为无序分支中的凝固点与有序分支中的熔化点具有相同的可压缩率状态, 并且其自由能相等的关系. 其凝固点对应于随机致密堆积($ \phi \approx 0.64 $), 熔化点对应$ \phi \approx 0.68 $^[19].

我们之前的研究中建立了Edwards系综理论的热力学第零定律^[26], 在此基础上可认为在一个恒定的剪切幅度下, 此时无序支(亚)稳态与晶体稳态两者的可压缩率是相同的. 从而在可压缩率的两条分支中对应的稳态位置用同一$ \chi $标定. 在本文实验中, 晶态分支稳态位于$ \phi \approx 0.70 $. 无序分支的稳态难以直接达到, 因为当$ \phi $接近RCP即开始结晶. 考虑到近期双分散球体的循环剪切实验表明稳态体积分数对$ \gamma $的依赖较弱^[28], 本文将“理想”的无序稳态近似取为$ \phi =0.636 $(RCP附近). 一旦固定该无序稳态点与晶态稳态点满足$ \chi $相等, 即可确定晶态分支的$ \chi (\phi ) $, 结果如图5(b)所示.

根据可压缩率定义, 可将熵写为对体积分数(或体积)的积分形式:

其中$ {S}_{\text{r}} $为参考态熵. 在实际处理时, 通常通过对$ \chi (\phi ) $的拟合获得平滑积分. 对晶态分支, 取$ {\phi }_{\text{FCC}}= \text{π}/\left(3\sqrt{2}\right)\approx 0.74 $的“近乎完美晶体”作为参考态, 并令$ S=0 $(即$ S({\phi }_{\text{FCC}})=0 $), 结合拟合形式即可得到$ S(\phi ) $(图5(c)). 对无序分支, RCP态的熵长期存在争议^[29]. 在Edwards框架下, Jin和Makse^[19]提出无序与晶态阻塞态在“类熔化点与类凝固点”处满足平衡条件$ {\chi }_{\text{freezing}}={\chi }_{\text{melting}} $且$ {f}_{\text{freezing}}={f}_{\text{melting}} $, 其中f为自由能密度. 本文以临界体积分数$ {\phi }_{\text{c}} $作为两分支的连接点, 用以确定无序与晶态分支之间的潜热及相应自由能差.

在Jin和Makse^[19]工作中, $ {\phi }_{\text{m}}=0.68 $被提出作为类熔化点, 而类凝固点对应RCP^[9]; 但其物理依据并非完全确定. 本文实验中, 结晶在$ {\phi }_{\text{RCP}} $以下已经发生, 说明有序化起始并不必然锁定在$ {\phi }_{\text{RCP}} $. 同时, 也缺乏直接证据表明类熔化点在不同摩擦条件下固定为$ {\phi }_{\text{m}}=0.68 $. 事实上, 近期研究提示熵与自由能可能依赖摩擦^[30], 因此相变位置随摩擦变化在物理上是自然的.

结合结构分析结果, $ \phi \approx 0.68 $附近出现清晰的弯折与涨落峰值(见图5(a)、图4), 该特征既可能源于动力学效应, 也可能指示类熔化转变. 鉴于其与既有建议值接近, 本文取$ {\phi }_{\text{c}}=0.68 $(及对应$ {\chi }_{\text{c}} $)作为连接点, 用于确定无序与晶态分支之间的潜热(图5(c))以及由此得到的自由能密度(图5(d)). 潜热定义如下:

为便于表征, 定义Edwards自由能密度$ f= V(\phi )-\chi (\phi )S(\phi ) $, 其中$ V(\phi ) $为体积(或体积分数的等价量)对应的密度表述, $ S(\phi ) $为熵. 由此可对无序、晶态两分支在不同$ \phi $区间的自由能进行比较, 并与结构序参量的分段演化对应起来(图5(d)). 自由能曲线在$ {\phi }_{\text{c}} $附近出现清晰的相交与分支分离: 在较高$ \chi $区域内, 无序分支自由能更低, 因此体系不会发生结晶; 随着$ \chi $逐渐降低, 两条分支会在自由能相等的位置相交; 在这之后的更低的$ \chi $位置下, 此时体系会不可逆的进入晶态. 这一趋势与结构分析中晶态团簇持续生长、最终达到平台的现象相一致. 与此同时, $ {\phi }_{\text{c}} $附近的自由能差变化与体积涨落峰值相伴随, 指示转变区间的主要贡献来自晶态与非晶共存导致的体积分布展宽与结构异质性增强.

需要指出的是, 上述自由能的“稳定性比较”应理解为对外驱阻塞稳态的一种有效统计描述: 其定量结果不可避免地依赖参考态选取、分支拟合方式以及等可压缩率连接假设. 尽管如此, $ f(\phi ) $的分支结构仍提供了一个可检验的判据: 若等可压缩率假设成立, 则在相同剪切协议下, 不同初态或不同重复实验所达到的稳态应在$ \chi $与$ f $上表现出一致性; 反之, 自由能分支的偏移可用于量化协议依赖与非平衡效应的强弱. 该思路为后续通过调控剪切幅度、摩擦条件的实验来检验 Edwards框架的适用性提供了明确路径.

本研究利用X射线断层成像对准静态循环剪切下三维光滑球形颗粒体系的致密化与结晶过程进行了实验表征, 并将结构序参量与Edwards体积系综的统计变量相结合形成统一的表征框架. 结果表明: 体系体积分数随循环次数逐步升高并在高密区达到稳态平台; 基于$ {q}_{6} $与$ w_{6}^{\text{hh}} $的分析在随机密堆积附近识别出显著结构拐折, 且晶态键统计给出的类固态颗粒分数随$ \phi $连续增长并在$ \phi \gtrsim 0.69 $进入饱和, 反映从无序态到共存态再到晶态主导态的分段演化. 体积涨落在无序分支随$ \phi $单调减小, 而在共存区间显著增大并出现峰值, 指示结构差异性增强; 由此得到的可压缩率、熵与自由能密度呈现分支分离, 说明循环剪切稳态在Edwards意义下对有序构型具有更高统计稳定性. 综上所述, 本研究为理解循环剪切下无热颗粒体系由无序致密化向剪切诱导结晶的转变, 并为剪切驱动堵塞态的统计力学描述提供实验依据.