颗粒物质是由大量宏观颗粒组成的离散体系^[1], 广泛存在于自然界与工程领域, 如土壤、砂石、谷物及药片等. 这类体系表现出独特的气、液、固三相行为及相应的相变特性^[2–5], 其力学响应显著区别于连续介质流体, 因此无法直接采用经典流体力学模型来描述物体在其中运动时所受阻力. 研究物体在颗粒介质中的受力机制, 对材料科学、土木工程和生物医学工程等多个领域具有重要的理论与应用价值^[1,6–9].

在早期研究中, 研究者们对颗粒介质中运动物体的阻力开展了一系列宏观力学实验, 发现物体低速运动时, 阻力与速度基本无关, 而主要受物体埋深、尺寸及材料等因素影响^[10]. 此外, 阻力随时间呈现周期性波动, 研究表明该现象源于介质内部堵塞状态周期性失效, 而非由运动物体与颗粒的界面特性所决定^[10–14]. 这些研究多集中于物体在颗粒介质中水平运动或从外部侵入过程^[15,16], 而对物体从颗粒介质中提升的非对称过程关注较少. 同时, 相关宏观实验所获结论多基于流体动力学模型或经验推测, 尚缺乏从微观尺度揭示其物理机制的系统分析^[17].

此后, Ding等^[18]围绕颗粒介质中物体所受阻力的定量计算展开研究, 提出了适用于物体水平运动的阻力模型, 该模型基于对颗粒流场的分析确定了剪切带位置, 并对剪切带所包围的未变形区域进行受力分析, 从而推导出作用于物体上的阻力表达式. 然而, 该模型仅限于描述水平运动情况. 为拓展至更复杂的运动形式, Li等^[19]进一步提出了新的阻力预测模型. 该模型的核心假设是: 物体所受合力与其埋深呈线性关系, 且复杂形状物体的总阻力可通过对其表面微元受力的积分获得. 研究结果表明, 该模型能较准确地预测物体在浅层颗粒介质中运动时的阻力; 但在深层条件下, 阻力与埋深之间并非线性关系, 导致模型预测精度显著下降. 总体而言, 现有研究主要关注水平运动或浅层条件下的阻力特性^[18–20], 对物体从深层颗粒介质中向上运动时所受阻力的机制仍缺乏系统探讨.

近年来, Athani和Rognon^[21]针对物体从颗粒介质内部垂直提升的过程开展了系统的离散元(DEM)模拟研究, 重点探究了不同因素对提升阻力的影响. 他们首先考察了提升加速度的作用, 通过模拟无黏性颗粒介质中水平板的动态提升过程发现: 在低速条件下($ v\leqslant \sqrt{gd}/10 $)阻力与速度无关, 主要由板上方锥形区域内颗粒的自重控制; 而在高速条件下($ v \gt \sqrt{gd}/10 $)峰值阻力可视为准静态阻力与惯性阻力的叠加, 其中惯性阻力与最终提升速度呈线性关系. 研究进一步提出了弹塑性惯性阻力模型, 揭示了由弹性波传播主导的颗粒逐步动员机制, 并指出较短的加速时间可能引发颗粒流化, 从而导致峰值阻力后的急剧下降^[21]. 随后, 他们进一步研究了颗粒间黏附力的影响, 通过模拟黏性颗粒介质中水平板的提升过程发现: 提升阈值力与颗粒间黏附力呈线性关系. 在浅埋条件下, 板与下方颗粒发生黏附分离, 需克服拉伸黏附力; 在深埋条件下, 板上方颗粒锥形区域发生垂直分裂, 黏附力对阻力的贡献与埋深相关. 整体阻力可表示为无黏性颗粒情形下的锥台重量与黏附引起的额外阻力之和, 并据此建立了量化两种机制的近似模型^[22]. 这些工作虽然深入地揭示了峰值阻力的形成机制及若干关键因素的影响规律, 但仍主要聚焦于特定条件或单一阶段的阻力分析, 对物体在整个提升过程中阻力演变的全阶段行为缺乏系统性研究.

本文通过实验与数值模拟相结合的方法, 研究了颗粒介质中圆柱体垂直提升过程的阻力特性, 获得了宏观阻力演化曲线及对应的颗粒尺度细观信息. 结果表明, 圆柱提升过程中所受阻力曲线呈现明显的三阶段特征: 初始阶段阻力迅速增大至峰值; 随后较快下降至某一水平; 之后进入缓慢衰减阶段. 基于颗粒流场分析及$ \mu $值(平均法向压力与等效剪切应力之比)分布的计算, 本文建立了相应的阻力计算模型, 模型计算结果与实验阻力曲线吻合良好. 进一步分析表明, 该三阶段行为与颗粒体系的细观结构演变密切相关: 初始阶段颗粒发生大规模重排, 随后逐渐形成稳定的承载结构, 从而对应阻力由快速上升、急剧下降转为缓慢波动的宏观响应.

实验所用颗粒与圆柱均采用光敏树脂3D打印技术制备. 其中, 球形颗粒直径为3 mm, 圆柱体直径为28 mm、高度为30 mm. 实验装置示意图如图1所示, 主体为一个内径120 mm的亚克力圆柱容器. 实验开始前, 首先调节支架高度, 使圆柱底面距容器底部保持10 mm; 随后将颗粒自容器上方均匀倒入, 填充至初始高度90 mm, 对应的堆积分数约为0.60. 圆柱体通过一根直径0.8 mm的直径细小且材质轻质高刚的碳纤维杆与高精度测力计连接, 形成硬连接. 其对颗粒流场的扰动及对测力结果的影响均可忽略不计.

实验启动后, 由电机驱动滑台带动圆柱以0.5 mm/s的速度匀速垂直上升, 同时通过高精度测力计(IMADA ZTS-DPU-5N)实时采集拉升过程中圆柱所受阻力数据. 该测力计的最小分辨率为0.001 N, 量程为0—5 N, 采样频率为100 Hz.

离散单元法(discrete element method, DEM)是研究颗粒体系力学行为的核心数值手段^[23], 该方法基于离散介质视角, 将颗粒体系视为由独立刚性或可变形颗粒组成的集合, 通过直接追踪每个颗粒的运动(位移、速度、加速度)及颗粒间的相互作用(接触力、摩擦、碰撞等), 从微观尺度迭代演化得到系统的宏观力学响应. 与宏观实验或连续介质模型相比, DEM的优势在于能够精确地获取颗粒尺度的细观信息, 如单个颗粒受力、力链演化、剪切带的形成与发展等, 从而揭示物体在颗粒介质中运动阻力的微观机理, 这是传统宏观实验难以实现的突破.

针对颗粒介质与圆柱体相互作用的复杂特性, 本研究基于自主开发的DEM程序开展模拟^[24], 该程序以高精度还原颗粒系统的力学行为为目标, 不仅能够计算颗粒间接触力与颗粒运动轨迹, 还可同步输出力链分布、颗粒位移场与应力场等关键细观参数, 为建立宏观阻力曲线与微观物理机制之间的关联提供了系统的数据基础.

颗粒间的接触力学模型是决定DEM模拟精度的关键. 本研究基于赫兹-明德林理论描述颗粒法向弹性接触行为, 并在此基础上引入库仑摩擦定律刻画切向滑动效应, 从而构建了法向-切向耦合的完整接触模型.

1)赫兹-明德林接触理论的适用性与计算逻辑

赫兹-明德林理论是弹性颗粒接触力学的经典模型, 适用于描述两个球形颗粒间的非黏性弹性接触, 其核心优势在于同时考虑了法向接触的弹性变形和切向接触的剪切变形, 且能够准确地反映接触力与变形量之间的非线性关系, 与颗粒介质的实际接触特性高度吻合^[23,25]. 在法向接触力计算中, 该理论基于弹性力学推导, 认为法向接触力$ {F}_{\text{n}} $与两颗粒的法向重叠量$ {\delta }_{\text{n}} $满足非线性关系:

其中$ {E}^{*} $为等效弹性模量, $ {R}^{*} $为等效半径, 其计算公式分别为

其中, $ {E}_{1},\; {E}_{2} $为两颗粒的弹性模量, $ {v}_{1},\; {v}_{2} $为泊松比, $ {R}_{1},\; {R}_{2} $为颗粒半径. 本研究中所有颗粒与圆柱体的材料特性一致, 弹性模量$ E=1.5\;\text{GPa} $, 摩擦系数$ \mu =0.4 $, 因此等效弹性模量和等效半径可进一步简化计算, 确保接触力传递的一致性. 在切向接触力计算中, 明德林在赫兹理论基础上进行了拓展, 考虑了颗粒接触点的剪切变形$ {\delta }_{\text{t}} $, 切向接触力$ {F}_{\text{t}} $与剪切变形满足线性关系$ {F}_{\text{t}}={k}_{\text{t}}{\delta }_{\text{t}} $ ($ {k}_{\text{t}} $为切向接触刚度), 同时引入黏弹性效应, 以模拟颗粒接触过程中的能量耗散, 避免因刚性碰撞导致的计算不稳定.

2)库仑摩擦定律的约束机制

为了准确地描述颗粒间的相对滑动行为, 本研究在切向接触力计算中引入库仑摩擦定律作为约束条件^[26,27]: 即颗粒间的切向接触力不能超过最大静摩擦力, 当$ \left| {F}_{\text{t}}\right| \gt \mu {F}_{\text{n}} $时, 颗粒间发生相对滑动, 此时切向接触力保持为$ {F}_{\text{t}}=\mu {F}_{\text{n}} $, 方向与相对滑动方向相反; 当$ \left| {F}_{\text{t}}\right| \leqslant \mu {F}_{\text{n}} $时, 颗粒间无滑动, 切向接触力随剪切变形线性增长. 这一约束机制能够真实还原颗粒介质在受力时先弹性变形、后滑动摩擦的力学响应, 是模拟剪切带形成和阻力演化的关键. 此外, 考虑到颗粒的实际运动特性, 本研究在DEM代码中还引入了滚动阻力模型, 以抑制颗粒过度滚动导致的模拟失真, 确保颗粒堆积状态和运动模式与实际物理场景一致. 滚动阻力矩的大小与颗粒的滚动角速度成正比, 比例系数根据颗粒材料的表面特性确定, 最终通过数值调试验证, 确保其对整体力学行为的影响符合实际规律.

1)数值求解方法

本研究采用显式时间积分方法求解颗粒的运动方程. 对于每个颗粒, 根据牛顿第二定律$ {m}_{i}{\ddot{x}}_{i}=\displaystyle\sum {F}_{ij}+{G}_{i} $($ {m}_{i} $为颗粒质量, $ {\ddot{x}}_{i} $为加速度, $ {F}_{ij} $为颗粒$ i $受到的其他颗粒及边界的接触力, $ {G}_{i} $为重力), 通过时间步长迭代计算颗粒的速度和位移. 为保证计算稳定性, 时间步长的选取遵循瑞利时间步准则, 即取系统瑞利时间步的$ 1/5\hbox{---}1/10 $, 本研究中最终设定时间步长为$ 1\times {10}^{-6}\;\text{s} $, 既确保数值稳定性, 又兼顾了计算效率.

2)系统边界与初始状态制备

模拟系统采用圆柱形刚性边界, 直径为120 mm, 边界材质与颗粒一致, 颗粒与边界的接触作用同样遵循赫兹-明德林接触理论和库仑摩擦定律, 以避免边界条件引入的模拟误差. 颗粒为均匀球形, 直径$ 3\;\text{mm} $, 圆柱体尺寸为直径$ 28\;\text{mm} $、高度$ 30\;\text{mm} $, 二者材料特性完全一致, 摩擦系数0.4, 弹性模量$ 1.5\;\text{GPa} $.

初始状态制备分为两个阶段: 首先将圆柱体固定于距离边界底部$ 10\;\text{mm} $处, 然后通过重力沉积将颗粒从系统顶部均匀撒落, 颗粒在重力作用下自由堆积; 随后进入弛豫阶段, 保持圆柱体固定, 让颗粒系统在重力和接触力作用下自由运动$ 10\;\text{s} $, 直至系统达到力学平衡, 判断标准为: 所有颗粒的平均速度低于$ 1\times {10}^{-6}\;\text{m/s} $, 且系统总动能占总势能的比例小于$ 1{\text{%}} $. 弛豫完成后, 圆柱体上表面距离颗粒介质上表面的距离为$ 50\;\text{mm} $, 此时颗粒体系的体积分数$ \phi $稳定在$ 0.59\hbox{---}0.61 $之间, 与实际颗粒堆积的致密状态相符.

3)加载与数据采集

弛豫完成后, 对圆柱体施加恒定的垂直向上提升速度($ 0.5\;\text{mm/s} $), 直至圆柱体上表面到达颗粒介质上表面. 在加载过程中, DEM代码实时采集关键数据, 包括: 圆柱体受到的总阻力(通过计算圆柱体表面所有接触颗粒的作用力得到)、每个颗粒的位置、速度、位移及所受接触力、颗粒间的接触关系(用于构建力链网络)等. 数据采集频率为每迭代1000步记录一次, 最终得到阻力曲线、颗粒流场、力链演化等数据, 为后续宏观力学分析和微观机制探究提供了全面的基础数据.

实验与模拟获得的阻力曲线对比如图2所示, 为了验证实验的可重复性, 我们在相同条件下进行了三次实验, 三次实验的结果几乎完全重合, 红色曲线为三组实验得到的平均阻力曲线. 蓝色曲线为基于下文所提阻力模型的计算结果, 我们将在3.3节详细介绍. 由于模拟中忽略了部分实际因素的影响, 导致模拟曲线与实验曲线在定量上存在一定偏差, 但两者在定性上呈现出高度一致的三阶段演化特征: 第一阶段阻力快速累积, 在上升约7 mm时达到峰值; 随后进入第二阶段, 阻力迅速下降, 到达某一点后下降速率明显减缓, 并过渡到第三阶段的缓慢变化. 为了确定第二阶段与第三阶段的分界点, 首先对峰值之后的阻力曲线进行Savitzky-Golay平滑处理, 进而计算平滑后曲线的斜率k随提升高度的变化(如图2插图所示). 结果表明, 斜率发生明显转折的位置对应上升高度约为22 mm处, 因此将此点定义为阻力曲线第二阶段与第三阶段的分界点.

首先对圆柱上升过程中的颗粒流场进行分析. 考虑到系统的中心对称性, 为清晰地展示提拉过程中颗粒的运动状态, 沿圆柱轴线方向截取了20个厚度为1.2 mm的垂直切片. 计算每一提拉步后各切片内颗粒的竖向(沿提升方向)位移, 经插值处理后对20个切片取平均值^[28], 得到不同提升阶段的颗粒流场分布, 如图3所示. 流场结果显示, 在圆柱上升过程中, 受扰动的颗粒流场呈放射状分布, 位移量由圆柱表面向外围逐渐衰减. 不同上升高度对应的流场特征存在明显的差异: 在初始上升阶段(加载阶段), 仅圆柱正上方局部区域的颗粒位移较大, 其余区域位移较小; 上升至约10 mm后, 圆柱上方较大范围内的颗粒均开始出现显著位移; 而当上升距离达到约22 mm后, 仅有靠近圆柱的较小范围内的颗粒继续跟随上升, 位移区域明显收窄. 进一步对力链演化进行分析(图4). 结果表明, 力链的断裂与重组主要发生在颗粒屈服面上. 圆柱正上方高位移区域内的力链密度最高, 且力链方向大多垂直于圆柱顶面, 直接传递向上的提拉力; 而在外围低位移的锥形区域内, 力链密度较低, 力链呈放射状分布, 主要承担剪切传递作用. 此外, 为从微观力学角度更精确地识别剪切带位置, 我们计算了参数$ {\mu }^{*}=\tau /P $. 其中$ \tau $为等效剪切应力, $ P $为平均法向压力, 该参数反映了剪切作用相对于法向压力的相对强度, 当$ {\mu }^{*} $超过屈服摩擦系数时, 颗粒将发生相对滑动. 由图5所示结果可知, 系统的屈服摩擦系数约为0.5. 在圆柱正上方及锥形区域边界处, $ {\mu }^{*} $值呈现由低至高的过渡, 该过渡区对应着颗粒由稳定状态向流动状态的转变. 综合以上分析, 可对阻力曲线达到峰值后先急剧下降再趋于平稳的现象给出如下解释: 阻力达到峰值时, 系统尚未形成稳定的剪切带, 颗粒发生大规模重排, 导致阻力剧烈波动; 当圆柱提升至约22 mm后, 其上方颗粒形成相对稳定的承载结构, 运动颗粒主要集中在圆柱上方的小锥形区域内, 并出现稳定的剪切带, 此时阻力的变化主要由剪切带处的摩擦行为所主导, 因而呈现平缓波动特征.

接下来, 我们对圆柱所受阻力展开定量建模. 我们认为, 圆柱在提升过程中所受阻力主要由两部分组成: 一是圆柱侧壁与颗粒间的摩擦阻力, 二是圆柱正上方颗粒对其施加的压力. 对于圆柱侧壁的摩擦阻力, 本文借鉴流体力学中经典的侧压力假设, 认为颗粒对圆柱侧壁的压力等于该处颗粒上方小圆柱区域内颗粒的自重. 摩擦阻力$ {F}_{\text{f}} $的计算公式如下:

其中, $ n $为与圆柱侧面接触的颗粒数; $ \mu $为摩擦系数; $ l $表示从圆柱中心至颗粒上边界的垂直距离; $ \phi $为颗粒体系的体积分数; $ r $为颗粒半径; $ \rho $为颗粒的密度. 圆柱上方颗粒所施加的压力$ {F}_{\text{G}} $由下式描述:

其中$ N $为系统中所有颗粒总数; $ {K}_{i} $为第$ i $个颗粒对压力的贡献系数. 为反映不同位置颗粒对阻力贡献的差异, 我们根据其竖向位移设定$ {K}_{i} $的取值范围: 位移大于区间上限$ {d}_{\max } $(取为 0.01 cm)的颗粒, 其$ {K}_{i} $取值为1; 位移小于区间下限$ {d}_{\min } $(取为0.001 cm)的颗粒, 其$ {K}_{i} $取值为0; 位移介于两者之间的颗粒, 其$ {K}_{i} $的取值则按位移大小线性映射至$ \left[0, 1\right] $. 该阈值$ {d}_{\min } $的设定旨在排除因误差引起的微小位移干扰. 此方法物理上对应于: 圆柱正上方区域的颗粒将其全部自重作用于圆柱, 而外围锥形区域内颗粒的贡献随其与中心距离增大而逐渐衰减. 将上述两部分阻力相加, 即得到圆柱所受总阻力:

基于该模型计算得到的阻力曲线如图2蓝色曲线所示, 模型结果与模拟获得的阻力曲线吻合良好, 验证了该建模方法的有效性.

本研究通过实验与数值模拟相结合的方法, 系统地探究了圆柱体从颗粒介质中垂直提升过程中的阻力演变机制, 并建立了相应的定量计算模型, 克服了以往研究多依赖二维体系、对阻力全程演化缺乏系统定量分析的不足. 宏观上, 圆柱所受阻力呈现典型的三阶段演化特征: 在初始约7 mm提升阶段, 阻力迅速累积至峰值; 随后急剧下降并逐渐趋于稳定, 进入小幅波动阶段. 微观上, 综合颗粒流场、力链网络及$ {\mu }^{*} $参数分布的分析, 揭示了峰值阻力后颗粒体系由大规模重排向稳定剪切带形成的转变过程, 从而阐明了阻力先骤降后缓变的细观机理. 本研究为颗粒介质中运动物体的阻力估算提供了简明、定量的理论框架, 在土木工程与地质灾害防治等领域具有一定的应用潜力. 后续工作可在本模型基础上, 进一步考虑颗粒材料特性、非圆柱几何、动态提升速度等多种因素的影响, 以拓展模型的适用性与工程实用性.