重力驱动下的颗粒物质流动广泛存在于多种地球物理过程与自然现象^[1–4]. 然而, 这类现象的触发条件与演化机制仍难以预测. 从物理角度看, 颗粒物质是典型的非平衡耗散体系, 几乎不受布朗运动支配, 且由于摩擦与非弹性碰撞持续耗散能量, 同时其流态转变受惯性、摩擦、材料性质与边界条件等多因素耦合影响, 呈现显著复杂性与普适性. 特别地, 地震与雪崩等过程中常伴随迟滞现象, 其出现往往与失稳触发密切相关, 因此理解迟滞机制有助于更深入地把握地震与雪崩的发生条件与动力学过程^[5–8]. 关于迟滞现象的本质机理, 传统观点多将惯性效应视为颗粒体系产生迟滞的主要来源; 但也有研究指出, 颗粒介质的滞后行为更可能源于颗粒接触界面的摩擦作用及其引起的能量耗散^[1–4]. 除此之外, 外部振动、惯性强度、颗粒摩擦特性以及颗粒材料的弹性模量等因素均会对迟滞效应产生显著影响. 在微观层面, Perrin等^[9], Forterre和Pouliquen^[10], DeGiuli和Wyart^[11]进一步提出: 与弹性模量相关的接触形变同样是调控迟滞的重要机制之一. 具体而言, 当颗粒材料刚度较大时, 体系的响应更容易表现出明显的非单调性, 迟滞特征更加突出; 而当材料刚度较小时, 这种非单调性会逐渐减弱, 甚至趋于消失而进入稳流状态^[6,11].

基于旋转鼓实验的观察结果^[12], Henein和Watkinson使用弗劳德数$ Fr $这一无量纲参数表征了不同的流态, 其定义为$ Fr=\omega \sqrt{D/g} $^[10,13–18], $ \omega $为旋转速度, $ D $为鼓的直径, g为重力加速度. 弗劳德数$ Fr $表示转鼓系统中离心力与重力的比值^[19,20]. 已有研究表明, 在转鼓体系中, 动态休止角与等效摩擦系数会随弗劳德数的增大呈单调上升趋势^[21,22]. 为更定量地区分不同流态, Dai等^[23]进一步引入了基于堆积表面角度波动的形态学参数, 并构建了以弗劳德数与密度比为控制变量的流态相图. 此外随着弗劳德数的增大, 颗粒流动将发生一系列流态转变, 通常依次经历迟滞态即间歇性雪崩、稳流态以及离心态等典型流动形态. 具体而言, 在低弗劳德数条件下, 体系常表现出明显的滞后效应与间歇性雪崩; 在中等弗劳德数范围内, 颗粒进入连续且稳定的流动状态; 而当弗劳德数进一步增大至较高水平时, 颗粒会沿鼓壁贴附滑动, 逐渐进入离心状态^[21,22]. 迟滞态的出现表明颗粒体系的能量耗散具有间歇性特征, 其间歇性的能量积累与突发释放行为在物理机制上与地震过程中的成震事件具有一定相似性. 相较之下, 连续稳态流动对应于能量耗散的时间连续性, 体系中不存在显著的能量积累与突发释放过程. 然而, 尽管迟滞态与稳态流在物理机制和应用意义上具有明确区分, 现有研究尚未给出能够统一具有不同弹性模量的颗粒材料体系的迟滞–稳态流动转变判据^[16,24–31]. 为此, 本文探究弹性模量对颗粒物质流动状态的影响进行系统性研究, 从而统一具有不同弹性模量的颗粒的流态转变行为.

本研究通过实验与数值模拟系统研究了旋转鼓中颗粒体系从迟滞态向稳定连续流动转变的动力学机制. 实验设计涵盖多组参数组合, 包括不同颗粒半径、不同材料弹性模量以及不同旋转速度. 借助高速成像与图像处理技术, 我们分别在迟滞区与稳流区获得了一系列稳定、可重复的宏观可测量, 如自由表面倾角及其波动特征等, 并据此刻画流态转变的临界条件. 进一步地, 我们提出并验证了一个包含弹性模量效应的修正弗劳德数: 不同弹性模量的颗粒体系在相同的临界值$ F{r}_{\text{c}} ^{**} $发生迟滞-连续流动转变, 实现了对材料差异的有效“归一化”描述. 因此, 颗粒材料的弹性模量能够通过本研究提出的修正弗劳德数统一表征其流态转变阈值, 将为雪崩、滑坡与泥石流等自然灾害过程的触发预测提供更强的物理依据与判据, 从而提升风险评估与预警的可操作性并降低潜在损失. 同时, 这一方向的进展也有助于深化我们对颗粒介质本构关系与非平衡流变行为的理解, 为建立更具普适性的颗粒流动模型提供实验与理论支撑.

本研究中所使用的实验装置如图1(a)所示, 尺寸为内径D = 184 mm, 宽度W = 95 mm的圆柱形转鼓安装在旋转电机上. 转鼓以不同恒定角速度$ {\varOmega }=0.4—56\;(^\circ)/{\rm s} $进行驱动. 为避免颗粒在内壁滑动, 鼓壁通过粘贴直径为5—6 mm的ABS半球体进行粗糙化处理. 转鼓内填充约一半体积的单分散颗粒, 颗粒的直径d = 3—6 mm, 总体积约555 mL. 颗粒的运动由一台垂直于转鼓端面安装的高速摄像机记录. 旋转过程中可实时记录颗粒堆自由表面的演化情况. 如图1(b)所示, 实验选取了两种不同材料的颗粒, 分别是玻璃珠和ABS颗粒. 每种材料采用三种粒径d = 3, 4, 6 mm. 表1是实验颗粒的各种属性.

首先, 利用高速摄像机获得实验过程的原始灰度图像, 如图2(a)所示. 随后, 采用基于深度学习的图像处理技术对转鼓内的颗粒进行识别, 并以不同颜色对颗粒进行标记, 如图2(b)所示. 在完成颗粒识别后, 通过算法函数提取颗粒的空间坐标信息, 如图2(c)所示. 另外, 通过采用全局阈值法对图像进行二值化处理, 实现了颗粒与背景的有效分离; 从图像中提取所有连通区域的边界, 并通过比较每个边界的点数, 选择出点数最多(即包围区域最大的)的边界, 通过距离图像中心的阈值来筛选边界点, 保留有效的边界, 如图2(d)所示. 最终, 对识别得到的表面坐标点进行拟合, 并计算每一帧图像对应的堆积表面倾角.

本研究同时使用离散单元法(discrete element method, DEM)来研究颗粒材料的流态转变行为, 以便能够获取系统中详细的颗粒尺度信息. 在离散元数值模拟中, 颗粒间的接触力采用Hertz-Mindlin接触理论进行计算, 而阻尼成分则基于Taberlet等^[32]所概述的推导得出. 在此模型中, 切向接触力也遵循库仑摩擦定律. 质量为$ {m}_{i} $, 惯性矩为$ {I}_{i} $的颗粒$ i $的运动方程为

式中, $ {v}_{i} $和$ {\omega }_{i} $分别为颗粒$ i $的平移速度和角速度; $ {F}_{ij}^{\rm n} $, $ {F}_{ij}^{\rm t} $和$ {M}_{ij} $分别为颗粒$ j $或鼓壁作用在颗粒$ i $上的法向接触力、切向接触力和滚动摩擦力矩; $ {F}_{i}^{\text{g}} $为颗粒$ i $所受重力, $ {r}_{ij} $为颗粒$ i $指向颗粒$ j $的位置向量. 法向接触力$ {F}_{ij}^{\rm n} $和切向接触力$ {F}_{ij}^{\rm t} $, 可表示为

式中, $ k $和$ \eta $分别为弹簧刚度和阻尼系数; $ \delta $为颗粒形变量; $ {v}_{ij} $为颗粒$ i $和颗粒$ j $的相对速度; $ {\mu }_{\text{s}} $为滑动摩擦系数, n和t (上标或下标)分别为法向分量和切向分量. 相关参数具体形式如下:

式中, $ E $, $ \nu $, $ R $, $ e $, $ G $, $ {\mu }_{r} $和$ \overline{{\omega }_{i}} $分别为颗粒弹性模量、泊松比、半径、质量、恢复系数、剪切模量、滚动摩擦系数和单位角速度矢量. 另外, 具体模拟参数在表2列出.

堆积角度计算方法: 筛选已知空间位置的边缘目标颗粒, 调用自定义函数拟合计算堆积表面倾斜角. 通过滑动窗口提取颗粒堆积上表面点, 利用主成分分析(principal component analysis, PCA)表面点的主延伸方向, 计算得到堆积表面的倾斜角度. 图3展示转鼓内颗粒体系在模拟过程中形成的瞬时流动状态. 图中颗粒按照速度大小进行着色, 颜色由蓝色逐渐过渡到红色, 表示颗粒速度由低到高. 可以看出, 体系内部大部分颗粒处于相对低速或近似静止状态, 而靠近自由表面的颗粒具有较高速度, 形成了一条明显的表面流动层. 这说明在转鼓运动过程中, 颗粒的流动行为主要集中于自由表面区域, 即类流体状态主要出现在表层. 该模拟结果直接反映了颗粒体系中静止区和流动区的空间分布, 为后续分析流态转变、临界转速以及颗粒尺寸和材料性质对流动行为的影响提供了基础.

随着转鼓角速度$ {\varOmega } $的增加, 颗粒堆积的流动形态由低速下的间歇性崩塌状态逐步过渡到连续流动状态. 如图4(a)所示, 给出了同一种材料在较低转速条件下, 堆积表面倾角$ \theta $随时间$ t $的演化规律. 可以看到, $ \theta \left(t\right) $呈现出典型的锯齿状波动特征: 表面倾角在随时间缓慢增大后发生一次突降, 对应一次局部崩塌事件. 基于这一特征, 可在每一次崩塌过程中清晰地定义崩塌前的最大倾角$ {\theta }_{\text{m}} $以及崩塌后的回复倾角$ {\theta }_{\text{r}} $, 并引入角度差$ \Delta \theta =\overline{{\theta }_{\text{m}}}-\overline{{\theta }_{\text{r}}} $, 在物理意义上, $ \Delta \theta $可作为表征间歇流动“迟滞强度”的一个序参量: $ \Delta \theta $越大, 表明体系在发生崩塌前需要积累更大的倾角超调, 因而越不易触发崩塌事件; 相反, 如图4(b)$ \Delta \theta $越小则意味着体系逐渐接近连续流动的极限情形. 特别地, 对于理想无摩擦颗粒体系, $ \Delta \theta \approx 0 $, 对应几乎不表现出迟滞效应, 这和Perrin等^[9], Forterre和Pouliquen^[10], DeGiuli和Wyart^[11]的结论一致.

图5(a), (b)给出了不同材料和不同颗粒半径角度差$ \Delta \theta $随转鼓角速度$ {\varOmega } $的变化关系. 实验结果表明, 随着角速度的增加, $ \Delta \theta $在低转速区间快速减小, 随后逐渐趋于稳定的平台值. 当$ {\varOmega } $增大至某一临界角速度$ {{\varOmega }}_{0} $附近时, $ \Delta \theta $基本达到稳定值$ \Delta {\theta }_{0} $. 值得注意的是, 在实验中可清晰观察到如下流态差异: 当$ {\varOmega } \lt {{\varOmega }}_{0} $时, 体系主要表现为间歇性崩塌对应的迟滞流态, 而当$ {\varOmega }>{{\varOmega }}_{0} $时, 颗粒堆积进入持续的稳定流动状态. 基于这一现象, 本文将$ {{\varOmega }}_{0} $定义为迟滞态向稳定流动状态的转变点. 此外, 从同种材料的数据对比可以发现, 随着颗粒半径增大, 转变点$ {{\varOmega }}_{0} $整体向更高角速度移动, 即较大颗粒需要更强的驱动才能进入稳定流动. 换言之, 颗粒半径越小, 越容易形成稳定连续流动. 与此同时, 平台值$ \Delta {\theta }_{0} $也随颗粒半径增大而增大, 表明在稳定流动区间内, 大颗粒体系仍保持更显著的角度波动幅度.

为进一步验证上述$ {{\varOmega }}_{0} $的确可作为流态转变点, 对倾角$ \theta $的时间序列进行了概率分布统计. 图6给出了不同转速$ {\varOmega } $条件下$ \theta $的概率分布函数(probability distribution function, PDF). 实验结果表明: 在低转速区间, 概率分布呈现较矮且较宽的双峰分布; 而在较高转速下, 概率分布逐渐演化为较高且较窄的单峰分布. 该分布形态的变化反映了体系动力学状态的改变: 低转速时体系在两种倾角水平之间交替, 对应崩塌前后的两类状态, 从而产生双峰特征, 高转速时体系围绕某一稳定平均倾角进行涨落, 表现为单峰分布.

从相变理论的类比角度看, 当系统处于单一相态时, 某些可观测量的涨落往往近似呈现单峰高斯分布. 而在相共存或两态竞争的情况下, 分布可能表现为双峰特征. 基于这一类比, 图6中的分布演化进一步明确了在$ {{\varOmega }}_{0} $发生的相态转变: 当$ {\varOmega } \gt {{\varOmega }}_{0} $时, 体系进入稳定连续流动状态. 当$ {\varOmega } \lt {{\varOmega }}_{0} $时, 体系在两种亚稳态之间交替, 表现为间歇性崩塌的迟滞状态. 如图6(a)所示, 3 mm颗粒在转鼓转速大于$ {\varOmega }=4\;(^\circ)/\text{s} $时, 概率分布由双峰转变为稳定的单峰分布, 该转速可与图5(a)中得到的3 mm体系转变点$ {{\varOmega }}_{0} $相对应. 类似地, 图6(b)中4 mm颗粒在$ {\varOmega }=5.2\;(^\circ)/\text{s} $附近出现由双峰向单峰的转变, 与图5(a)中4 mm体系的转变点$ {{\varOmega }}_{0} $一致.

对于迟滞态向稳流态的转变过程, 其本质是转鼓旋转带中颗粒的向心力和重力的竞争, 因此, Henein等^[33]引入了无量纲弗劳德数$ Fr $对这一相态转变过程进行描述. 其表达形式为: 由于转鼓导致的向心力$ F $与重力$ G $的比值, $ Fr=\varOmega \sqrt{D/g} $. 在此基础上, Balmforth和McElwaine^[34]进一步对$ Fr $作了尺寸修正, 提出修正弗劳德数$ F{r}^{*}{=}Fr(D/d) $. 其中$ D $为转鼓直径, d为颗粒直径. 修正之后二者趋势虽然完全相同, 但如果你要比较: 不同鼓直径D, 不同颗粒尺寸d, 那么基础$ Fr $就不够, 因为它无法区分“同样$ Fr $但颗粒离散程度完全不同”的系统. 这时$ F{r}^{*} $的意义就出来了: 它试图把动力学强度和颗粒尺度效应统一到一个参量中; 更有希望让不同条件下的颗粒统一塌缩; 更适合做“小鼓实验”与“大尺度颗粒流/滑坡”之间的类比. 如图7为归一化系数$ \Delta \theta /\Delta {\theta }_{0} $随$ F{r}^{*} $在不同颗粒粒径条件下的关系, 转变点坍缩到统一的$ \Delta \theta /\Delta {\theta }_{0} {=} f(F{r}^{*}) $曲线位置. 对于玻璃颗粒, 不同半径样品均在$ F{r}^{*}= 0.6 $附近发生一致的流态转变, 即由迟滞状态进入稳定流动状态. 上述结果表明, 在引入$ D/d $的尺寸修正后, $ F{r}^{*} $能更有效地表征不同粒径体系的临界转变条件, 为流态转变的统一描述提供了更简洁的无量纲框架. 然而对于不同材料的颗粒并不具有普适性, 其转变点发生在不同的$ F{r}^{*} $.

如图8(a)所示对于玻璃颗粒, 不同半径样品均在$ F{r}^{*}=0.6 $附近发生一致的流态转变, 即由迟滞状态进入稳定流动状态. 而对于ABS颗粒, 该统一转变位置对应$ F{r}^{*}=0.72 $, 结果表明不同种类的材料在不同$ F{r}^{*} $发生统一转变.

对于同种材料的不同粒径体系, 上文已表明可以通过修正弗劳德数$ F{r}^{*} $对转变点进行较好的统一. 然而, 当材料弹性模量发生变化时, 各体系往往在不同的$ F{r}^{*} $处发生转变. 需要指出的是, 在实验条件下实现仅改变弹性模量而保持其他参数严格不变非常困难. 因此, 针对弹性模量效应的研究, 我们采用数值模拟方法开展系统对比. 为验证模拟结果的合理性, 基于数值模拟进行了ABS颗粒的实验复现, 图9(a)是实验和模拟结果得到的角度差$ \Delta \theta $随转鼓角速度Ω的变化关系. 可以看到, 无论实验还是数值模拟, $ \Delta \theta $均随Ω增大表现为显著下降并最终趋于稳定平台值$ \Delta {\theta }_{0} $. 这一总体趋势具有一致性. 尽管实验与模拟曲线在数值上并不完全重合, 但二者的演化形态相同. 造成数值偏差的主要原因在于实验系统中存在若干难以精确获取或在模拟中完全复现的参数, 如颗粒层面的表面粗糙度、实验用颗粒的形状的不均匀性, 以及转鼓内壁的粗糙/摩擦特性与碰撞耗散等. 进一步地, 为削弱参数不确定性带来的绝对量差异, 我们对每组数据进行归一化处理: 将纵坐标归一化为$ \Delta \theta /\Delta {\theta }_{0} $, 横坐标归一化为$ {\varOmega }/{{\varOmega }}_{0} $(其中$ {{\varOmega }}_{0} $为对应体系的转变点的转速, $ \Delta {\theta }_{0} $为稳态流的角度差). 如图9(b)所示, 经归一化后, 实验和模拟的结果塌缩到$ \Delta \theta /\Delta {\theta }_{0}= f\left({\varOmega }/{{\varOmega }}_{0}\right) $曲线位置. 这一结果表明实验与模拟在描述迟滞到稳流转变的主导机制上具有一致性, 从而在一定程度上验证了数值模拟方案的合理性, 并为后续引入弹性模量修正、讨论跨材料的统一转变提供了可靠基础.

因此, 基于经过验证的离散元数值仿真模型对不同弹性模量的颗粒进行了模拟, 保持颗粒半径和密度不变, 系统性改变模量和转鼓转速. 图10(a)展示了不同弹性模量材料的角度差$ \Delta \theta $随转鼓角速度$ {\varOmega } $的变化关系. 对于$ E=2.3\times {10}^{9} $的材料, $ \Delta \theta $随$ {\varOmega } $增加先快速下降并在$ {\varOmega }=14{\varOmega }\;(^\circ)/\text{s} $附近进入平台区, 表明体系由迟滞流态转入连续稳定流动流态. 其迟滞-稳态流转变点为$ {\varOmega }=14{\varOmega }\;(^\circ)/\text{s} $. 对于$ E=2.3\times {10}^{8} $的材料, 其转变点位于$ {\varOmega }=15{\varOmega }\;(^\circ)/\text{s} $. 实验结果显示, 随着材料弹性模量$ E $的降低, 转变点$ {{\varOmega }}_{0} $整体向更高角速度移动, 即体系需要更强的驱动才能进入稳定流动状态. 换言之, 颗粒材料越“软”, 要实现连续稳定流动越困难, 这表明材料刚度对迟滞到稳流转变具有显著调制作用.

弗劳德数$ Fr={\varOmega }\sqrt{D/g} $表征转鼓体系中惯性效应相对于重力效应的相对强度. 在本研究中, 由于不同模量会导致颗粒间的形变量不一致. 以此我们在Fr的基础上引入含模量的修正项, 提出修正弗劳德数:

其中$ \delta $是形变量, r是特征长度颗粒半径, $ \alpha $为待定指数. 该修正的目的在于不同模量材料之间实现转变点的统一表征. 在赫兹模型中接触应力F_n为

其中, $ E $为弹性模量. (13)式应理解为统计平均意义下的标度关系, 而不是对每一个接触形变的精确描述. 不同接触之间的形变差异会导致一定的分布宽度, 但在无量纲分析中, 这些差异通常只会影响数值系数, 而不会改变主要的幂律标度. 从这个角度看, 所使用的$ \delta /r $可以理解为接触网络中的有效平均形变量.

在实际颗粒体系中, 一个颗粒通常同时参与多个接触, 这些接触共同承担颗粒所受的法向载荷. 因此, 颗粒的重力并不是由某一个接触力单独平衡, 而是由若干接触力的合力共同承担. 在统计平均意义下, 接触网络中的典型接触力尺度仍然由宏观载荷尺度(例如颗粒自重或流动层中的典型法向应力)决定. 力尺度$ {F}_{n} $应理解为接触网络中典型接触力的数量级估计, 而不是单一接触的精确力值. 更严格地说, 颗粒重力应由若干接触力分担, 因此单个接触的典型力规模应为$ {F}_{n}\sim {mg}/{Z} $, 其中$ Z $是平均接触数. 由于$ Z $在同一体系中通常只在有限范围内变化, 而不会改变主要标度关系, 因此其接触力$ {F}_{n}\propto mg $,

代入(13)式整理可得

整理(15)式:

基于上述讨论, 最后整理得修正弗劳德数:

其中$ \beta =\dfrac{2}{3}\alpha $, $ \rho $为粒子密度, $ g $重力加速度.

如图10(b)所示, 通过对数据分析, 取$ \alpha = -0.052 $, $ \beta =-0.0344 $将不同弹性模量材料的数据按修正弗劳德数$ F{r}^{**} $重新标度, 并得到角度差$ \Delta \theta /\Delta {\theta }_{0} $随$ F{r}^{**} $的变化关系. 随着$ F{r}^{**} $的增大, $ \Delta \theta /\Delta {\theta }_{0} $先出现显著下降, 随后逐渐进入平台区并趋于稳定, 表明体系由间歇性崩塌主导的迟滞流态逐步过渡到连续稳定流动流态. 更为关键的是, 不同弹性模量材料的迟滞-稳态流呈现统一的转变点$ F{r}^{**}=F{r}_{\mathrm{c}} ^{**}$, 各体系基本在$ F{r}^{**}=8.8\times {10}^{-3} $附近完成由迟滞态向稳流态的转变. $ F{\text{r}}^{**} $的物理意义在于, 它不是单纯的经验拟合参数, 而是一个把“宏观驱动”和“微观接触力学”联系起来的半经验判据. 对于同种材料、不同粒径的体系, $ D/{d} $反映的是几何离散尺度对转变的影响; 对于固定几何、不同模量的体系, $ {\left(\delta /r\right)}^{\alpha } $反映的是接触形变对转变的影响. 两者共同说明, 迟滞-稳流转变并不只由惯性与重力的竞争决定, 还取决于颗粒网络在加载下积累、释放和连续耗散能量的能力. 最后, 该结果说明, 引入模量修正项后, 修正弗劳德数$ F{r}^{**} $能统一不同刚度材料的迟滞-稳态临界转变, 为跨材料体系的流态转变表征提供了统一判据.

本研究围绕转鼓中颗粒堆积由迟滞态向连续稳定流动态的转变开展了系统实验研究, 揭示了颗粒尺度与材料属性对转变点的影响规律. 实验结果表明, 不同半径的同种材料颗粒, 其转变点可通过修正弗劳德数实现较好的统一: 体系在近似相同的修正弗劳德数取值处由迟滞态转入稳流态. 进一步地, 通过对不同转速下倾角时间序列的概率分布统计, 我们观察到分布由低转速区间的双峰特征向高转速区间的单峰特征演化, 从统计意义上为转变点的存在提供了独立证据. 此外, 数值模拟结果显示, 对于不同刚度的材料体系, 颗粒堆积同样可在某一修正弗劳德数附近由迟滞状态过渡到稳定流动状态, 表明该无量纲表征在一定程度上具有跨材料的适用性. 上述结果为建立更简洁的转鼓流态判据与临界条件提供了实验与数值依据, 并有望为雪崩、滑坡等相关自然现象的阈值预测与风险评估提供参考. 同时, 本工作的发现也为进一步理解颗粒介质的宏观本构行为与其微观机制之间的联系提供了新的思路.