粒径分布特征对水下颗粒柱坍塌过程影响

Effect of Particle Size Distribution on Collapse of Immersed Polydisperse Granular Columns

摘要

基于典型的颗粒柱坍塌试验，采用计算流体力学与离散单元法耦合的数值方法，系统探讨了以分形维数表征的粒径分布特征在不同流态下对颗粒柱坍塌过程的流动能力和能量演化的影响机制.结果表明：随着环境流体由空气逐步过渡至低粘度流体和高粘度流体，颗粒柱的运动距离相比干燥条件降低约40%，流动能力显著减弱；在自由下落态和惯性态中，不同分形维数系统的流动性差异仅约1%，而在粘性态中，较高分形维数的系统表现出明显的运动延迟且运动能力降低幅度达11%；并且这种差异被认为是细颗粒含量更多的系统更容易在高粘度流体中耗散能量引起的；利用渗透性测试方法发现，水下颗粒系统的运动能力主要受初始渗透性影响，渗透性越差运动能力越弱.

Abstract

In this study, a typical granular column collapse scenario is investigated using a numerical approach that couples computational fluid dynamics (CFD) with the discrete element method (DEM). Granular columns with varying fractal dimensions—used to characterize PSD—are simulated in different fluids to examine the influence of PSD on flow behavior and energy evolution across distinct flow regimes. The results show that as the ambient fluid gradually changes from air to low-viscosity fluid and then to high-viscosity fluid, the runout distance of the granular systems decreases by about 40% compared to dry conditions, and the mobility is significantly weakened. In free-fall and inertial regimes, the mobility difference between systems with different fractal dimensions is only about 1%, whereas in viscous regimes, the system with a higher fractal dimension shows a significant movement delay and the mobility is reduced by 11%. This reduced mobility can be attributed to the increased presence of fine particles, which enhance energy dissipation under low Stokes number conditions. Permeability tests further reveal that mobility is primarily governed by the initial permeability of the immersed granular system-lower permeability corresponds to reduced mobility.

Graphical abstract

关键词

颗粒材料 / 分形维数 / 海底滑坡 / 颗粒柱坍塌 / 渗透性 / 工程地质学.

Key words

granular material / fractal dimension / submarine landslide / column collapse / permeability / engineering geology

0 引言

随着我国海岸带与海洋资源快速开发建设，海底滑坡等地质灾害风险日益凸显.海底滑坡不仅自身具有强烈的破坏力，更可能通过触发次生灾害如涌浪、浑水流等事件，演变为典型的链式灾害过程，从而对水下关键基础设施（管线、光缆、基础平台）安全构成重大威胁（Hsu et al.， 2008；年廷凯等，2021；彭铭等，2025）.如渤海浅海层发生的大规模海床滑塌事件，导致临近的多处输油管道被切断，引发严重的原油泄漏、环境污染等问题（吴时国等，2019）.因此，深入揭示水下滑坡运动演化机制及其与复杂海洋环境（如沉积物特性、水动力条件）的相互作用机理，对于精准预测滑坡灾害链、提高深远海地质灾害风险防控能力和灾害韧性具有重要意义.

由于与自然条件下地质灾害过程的动力学行为高度相似，颗粒柱坍塌常被视为地质灾害过程的一种简化构型，是深入研究碎屑体颗粒流复杂瞬态流动演化过程的重要手段（Lube et al.， 2004； 崔溦等，2021；满腾等，2022）.在水下环境中，颗粒柱坍塌涉及颗粒-流体耦合、孔压动态响应及运动演化等多个物理过程，有助于理解库岸及海底滑坡等地质灾害的触发机制与运动特征（景路等，2019；Yang et al.， 2020）.Courrech du Pont et al.（2003）指出，水下颗粒柱的宏观运动行为不仅受颗粒属性控制，还强烈依赖于所处的流体环境.斯托克斯数（St）和流体-颗粒密度比（r）是描述颗粒-流体耦合关系的两个关键无量纲参数，其中St与颗粒和流体密度、颗粒粒径、流体粘度有关.利用这两个无量纲参数能有效划分水下颗粒柱的运动状态：（1）自由下落态（FF），在该状态下周围环境流体的阻力可以忽略不计，一般对应于颗粒与空气相互作用的干燥状态；（2）惯性态（I），此时颗粒在环境流体中主要受到流体惯性力的影响；（3）粘性态（V），此时颗粒在环境流体中主要受到流体粘性力的影响（Bougouin and Lacaze， 2018； Jing et al.， 2019）.

自然界中的颗粒物质普遍具有一定的粒径分布特征（PSD），表现出明显的多分散特性（Breard et al.， 2023）.在处于自由下落态的干燥条件下，颗粒材料的流动性在没有尺寸效应影响时，被证明与材料的多分散特性无关（Cabrera and Estrada， 2021； Polanía et al.， 2022）.而当颗粒材料处于惯性态时，存在一定粒径分布形式的多分散颗粒系统相较于粒径均匀的系统在动力学行为方面将产生显著差异，并且体积分数被认为是控制这种差异的关键参数（Di Felice， 1994； Trulsson et al.， 2012）.试验表明，体积分数通过改变颗粒系统内部孔隙水压力的演化过程，从而显著影响其流动行为（Iverson， 1997，2005）.因此，该机制导致在惯性态下，颗粒系统的坍塌速度和流动能力随多分散性及其引起的体积分数的升高而减弱（Pailha et al.， 2008； Wang et al.， 2017； Lee， 2021； Polanía et al.， 2022，2024）.并且相较于惯性态，多分散颗粒系统在粘性态下的孔隙结构对孔隙水压力的变化更为敏感，从而放大了体积分数对坍塌过程与流动性的影响（Rondon et al.， 2011； Yang et al.， 2020； Lacaze et al.， 2021）.

尽管关于粒径分布特征对浸没颗粒柱动力学行为的影响已开展若干研究，并取得了一定进展，但是仍存在以下关键问题亟待解决：首先，目前针对粒径分布特征的研究多采用不同粒径组分均匀分布的简化形式（Polanía et al.， 2022，2024），与实际碎屑物质的自然粒径分布形式差异较大，从而限制了研究结论的代表性及其在工程实践中的适用性.其次，目前关于多分散性影响的研究普遍难以有效区分体积分数变化所带来的额外作用，因而无法准确评估粒径分布所引起的多分散性对流动性的独立影响.这主要是由于粒径分布的变化不可避免地会引起颗粒系统体积分数的改变（Oquendo-Patiño and Estrada， 2022），而体积分数已被广泛认为是控制水下颗粒流动行为的关键参数之一.最后，目前关于多分散颗粒柱坍塌的研究主要关注惯性态，尚未将研究拓展至粘性态，因此有必要进一步探究多分散性在该流态下对颗粒系统流动性的具体影响机制.

本文采用数值模拟方法，系统研究了不同流体环境下多分散颗粒柱的坍塌过程.通过对自由下落态、惯性态与粘性态的对比分析，揭示了以分形维数表征的粒径分布特征对颗粒系统流动性的影响.进一步从能量耗散与渗透性的角度出发，阐明了不同粒径分布系统流动能力差异的物理机制，并提出了一个考虑颗粒系统渗透性影响的流动距离标度律，用于量化多分散颗粒系统在水下环境中的运动能力.本研究旨在深化对粒径分布特征主导下的水下碎屑颗粒流动力学机制的理解，为自然颗粒流建模与地质灾害流动性预测提供理论基础和数值支撑.

1 流固耦合计算模型及工况设置

为模拟多分散颗粒柱在不同环境流体中的坍塌过程，本文采用计算流体力学（CFD）方法求解基于局部平均形式的纳维-斯托克斯方程以描述流体相的运动，同时采用离散单元法（DEM）求解牛顿运动方程以模拟颗粒相的动力学行为.本文基于开源软件OpenFOAM与LIGGGHTS分别进行CFD与DEM求解，并通过CFDEM框架实现两相间的流固耦合，耦合力包括阻力、浮力与粘性力等流体作用项（Jing et al.， 2016）.

1.1　流固耦合方法

针对颗粒系统，基于DEM求解控制平动和转动的牛顿方程，计算并更新每个颗粒的速度与位置（Cundall and Strack，1979）：

\left\{\begin{array}{c}{m}_i\frac{\mathrm{d}{\mathit{U}}_i^p}{\mathrm{d}t}={\sum }_{j=\mathrm{1}}^{n_i^c}{\mathit{F}}_{ij}^{\mathrm{c}}+{\mathit{F}}_i^{\mathrm{f}}+{\mathit{F}}_i^{\mathrm{g}}\\ I_i\frac{\mathrm{d}{\omega }_i}{\mathrm{d}t}={\sum }_{j=\mathrm{1}}^{n_i^c}\left({\mathit{M}}_{t,ij}+{\mathit{M}}_{r,ij}\right)\end{array}\right.，(1)

式中：m_i和I_i分别是颗粒i的质量和转动惯量，{\mathit{U}}_i^p和{\omega }_i分别是颗粒i的平动速度和角速度，n_i^c是颗粒i的总接触数，{\mathit{F}}_{ij}^c是颗粒i的总接触力，{\mathit{F}}_i^{\mathrm{f}}是颗粒i上的颗粒-流体相互作用力，{\mathit{F}}_i^{\mathrm{g}}是颗粒i的重力，{\mathit{M}}_{t,ij}和{\mathit{M}}_{r,ij}分别是作用在颗粒i上的切向力和滚动摩擦力所产生的力矩.其中，颗粒接触力F_{ij}^c基于Hertz-Mindlin接触模型计算：

{\mathit{F}}_{ij}^c=\left(k^{\mathit{n}}{\delta }_{ij}^{\mathit{n}}-{\gamma }^n{v}_{ij}^{\mathit{n}}\right)+\left[\left({\mathit{F}}_{\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}}^{s\mathrm{0}}+k^s∆{\delta }_{ij}^{\mathit{t}}\right)-{\gamma }^s{v}_{ij}^{\mathit{t}}\right]，(2)

式中：k^{\mathit{n}}{\delta }_{ij}^{\mathit{n}}表示法向弹簧接触力，{\gamma }^n{v}_{ij}^{\mathit{n}}表示法向阻尼力，{\mathit{F}}_{\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}}^{s\mathrm{0}}+k^s∆{\delta }_{ij}^{\mathit{t}}表示切向弹簧接触力，{\gamma }^s{v}_{ij}^{\mathit{t}}表示切向阻尼力.{\mathit{F}}_{\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}}^{s\mathrm{0}}+k^s∆{\delta }_{ij}^{\mathit{t}}代表总切向力，并且在运动中持续增长到\mu {\mathit{F}}^n，其中\mu是摩擦系数，{\mathit{F}}^n是第一项法向弹簧接触力.切向弹簧接触力会在随后的运动中保持不变直到颗粒间失去接触，{\mathit{F}}_{\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}}^{s\mathrm{0}}表示上一个时间步的初始切向弹簧力.k^{\mathit{n}}和k^s分别是法向和切向刚度，{\delta }_{ij}^{\mathit{n}}是法向重叠距离，∆{\delta }_{ij}^{\mathit{t}}是切向位移增量，{\gamma }^n和{\gamma }^s分别是接触面法向和切向上的阻尼系数并在后续计算中作为恢复系数e的函数，v_{ij}^{\mathit{n}}和v_{ij}^{\mathit{t}}分别是重叠颗粒相对速度的法向和切向分量.

针对环境流体，CFD采用对流体域的离散化处理，并在每个网格内求解连续性方程与局部平均的纳维-斯托克斯方程，以模拟环境流体的流动过程：

\left\{\begin{array}{c}\frac{\partial \left({\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}\right)}{\partial t}+\nabla \bullet \left({\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{\mathit{U}}^{\mathrm{f}}\right)=\mathrm{0}\\ \frac{\partial \left({\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{\mathit{U}}^{\mathrm{f}}\right)}{\partial t}+\nabla \bullet \left({\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{\mathit{U}}^{\mathrm{f}}{\mathit{U}}^{\mathrm{f}}\right)=\\ -\nabla p-{\mathit{f}}^{\mathrm{p}}+{\epsilon }_{\mathrm{f}}\nabla \bullet \mathit{\tau }+{\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}\mathrm{g}\end{array}\right.，(3)

式中：{\epsilon }_{\mathrm{f}}是孔隙率，{\rho }_{\mathrm{f}}是流体平均密度，{\mathit{U}}^{\mathrm{f}}是网格内流体的平均速度，p是网格内流体压力，{\mathit{f}}^{\mathrm{p}}是网格内颗粒对流体产生的相互作用力，\mathit{\tau }是粘性应力.计算过程中假设每个流体网格内的流体性质不变.

颗粒与流体之间的耦合作用是通过公式（1）中的流体作用力{\mathit{F}}_i^{\mathrm{f}}与公式（3）中的源项{\mathit{f}}^{\mathrm{p}}之间的信息交换实现.单个颗粒上的流体作用力可表示为：

{\mathit{F}}_i^{\mathrm{f}}={\mathit{F}}^{\mathrm{b}}+{\mathit{F}}^{\mathrm{d}}+{\mathit{F}}^{\mathrm{v}}，(4)

式中：{\mathit{F}}^{\mathrm{b}}，{\mathit{F}}^{\mathrm{d}}，{\mathit{F}}^{\mathrm{v}}分别是浮力、拖曳力与粘性力.每个流体单元网格内的颗粒-流体相互作用力为：

{\mathit{f}}^{\mathrm{p}}=\frac{\mathrm{1}}{V_{\mathrm{c}}}{\sum }_{i=\mathrm{1}}^{n_c}{\mathit{F}}_i^{\mathrm{f}}，(5)

式中：V_{\mathrm{c}}是流体网格的体积，n_{\mathrm{c}}是流体网格内颗粒的数量.浮力的计算公式为：

{\mathit{F}}^{\mathrm{b}}=\frac{\mathrm{1}}{\mathrm{6}}\pi {\rho }_{\mathrm{f}}{d}_{\mathrm{p}}^{\mathrm{3}}\mathbf{g}，(6)

式中：d_{\mathrm{p}}是颗粒粒径.拖曳力的计算公式为：

\left\{\begin{array}{c}{\mathit{F}}^{\mathrm{d}}=\frac{\mathrm{1}}{\mathrm{8}}{C}_{\mathrm{d}}{\rho }_{\mathrm{f}}\pi d_{\mathrm{p}}^{\mathrm{2}}\left({\mathit{U}}^{\mathrm{f}}-{\mathit{U}}_i^{\mathrm{p}}\right)\left|{\mathit{U}}^{\mathrm{f}}-{\mathit{U}}_i^{\mathrm{p}}\right|{\epsilon }_f^{-\chi }\\ C_{\mathrm{d}}={\left(\mathrm{0.63}+\frac{\mathrm{4.8}}{\sqrt[]{R{e}_{\mathrm{p}}}}\right)}^{\mathrm{2}}\\ R{e}_{\mathrm{p}}=\frac{{\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{d}_{\mathrm{p}}\left|{\mathit{U}}^{\mathrm{f}}-{\mathit{U}}_i^{\mathrm{p}}\right|}{{\mu }_{\mathrm{f}}}\\ \chi =\mathrm{3.7}-\mathrm{0.65}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\mathrm{1.5}-\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{10}}^{Re}p\right)}^{\mathrm{2}}}{\mathrm{2}}\right]\end{array}\right.，(7)

式中：C_{\mathrm{d}}是与雷诺数R{e}_{\mathrm{p}}有关的阻力系数，{\epsilon }_{\mathrm{f}}^{-\chi }作为校正系数，用于考虑系统中其他颗粒对计算颗粒阻力的影响，主要与孔隙率{\epsilon }_{\mathrm{f}}和由雷诺数反映的流动状态\chi有关.由流体剪切引起的粘性力的计算公式为：

{\mathit{F}}^{\mathrm{v}}=-\left(\nabla \bullet \mathit{\tau }\right)V_i^{\mathrm{p}}，(8)

式中：V_i^{\mathrm{p}}是颗粒i的体积.

1.2　粒径分布特征构建

在颗粒流动过程中，由于颗粒之间频繁发生摩擦和碰撞，常伴随着颗粒的破碎与细化，进而改变其粒径分布特征（尹振宇等，2012；彭建兵等，2023）.现场调查发现，碎屑颗粒流的沉积物通常表现出明显的逆级配特征，即颗粒尺寸由沉积层上部到底部逐渐减小（李坤等，2022；Shi et al.， 2024）.这种由破碎过程驱动的粒径细化趋势，其粒径分布（PSD）往往具有自相似特征，并可通过幂律关系N\left(d\right)=N_{\mathrm{0}}{\left(d/d_{\mathrm{0}}\right)}^{-D}利用分形维数D对这种分布进行定量描述（Lai et al.， 2021； Breard et al.， 2023）.其中，N\left(d\right)是直径为d的颗粒数量，N_{\mathrm{0}}是参考粒径为d_{\mathrm{0}}的颗粒数量（Hyslip and Vallejo， 1997）.较大的分形维数对应于系统中细颗粒含量的增加，反映了更剧烈的破碎过程.

鉴于碎屑流普遍呈现出的粒径分布特征，本文引入分形维数D对颗粒材料的多分散性进行定量建模，以构建更贴近自然条件的颗粒系统.实地调查表明，不同类型的碎屑流（如岩崩、冰岩混合流和泥石流等）所对应的分形维数大致分布在1.30~3.54之间（Crosta et al.， 2017）.本文选取D=1.0， 3.5代表两类极端粒径分布特征，以涵盖自然界中常见的多分散性变化范围.Tang et al.（2025）指出，分形维数D在1.0~3.5范围内变化时，多分散颗粒柱的坍塌动力学过程呈现出近似线性变化趋势.因此，本文选取具有代表性的两个极端分形维数（D=1.0和D=3.5）作为研究对象以揭示粒径分布对系统动力学行为的影响.本文采用半径为1 mm、2 mm、3 mm、4 mm、5 mm、6mm的6种不同尺寸颗粒，并根据幂律关系计算各粒径对应的颗粒数量，以生成具有不同分形维数的粒径分布.颗粒柱的高宽比设为1，试样尺寸设定为对应分形维数下颗粒平均粒径的50倍，以避免有限尺寸效应对结果的影响（Cabrera and Estrada， 2021）.图1展示了D=1.0和D=3.5的颗粒柱试样，以及相应的粒径数量分布与累计分布曲线.

1.3　工况及模型设置

图2展示了水下多分散颗粒柱坍塌计算模型示意图.模拟计算域为充满流体的矩形区域，其中流体域顶部为开放边界，在x方向和底部施加无滑移固壁条件，颗粒则采用固定边界条件处理.在坍塌流动的展向y方向对流体和颗粒均采用周期性边界条件，并且厚度为10倍平均粒径，以有效避免侧向边界效应的影响.

在本研究中，分形维数D=1.0时，自然堆积状态下的体积分数约为0.624；而当D=3.5时，体积分数提升至0.688.由于体积分数是影响水下颗粒流动力学行为的关键因素，因此本文采用固定体积分数（φ=0.667）的方法生成颗粒柱，以避免在研究粒径分布对水下颗粒流动力学过程的影响时耦合体积分数的作用.具体方法是，在颗粒柱填充区域内根据目标体积分数（φ=0.667）生成对应数量的颗粒.对于D=1.0的颗粒柱，该体积分数对应较为致密的堆积状态.而对于D=3.5的颗粒柱，该体积分数则相对较松散.因此通过改变重力加速度控制颗粒沉降的方法调整不同分形维数的颗粒柱内的颗粒排布，使颗粒柱高度达到预设值H₀从而满足目标体积分数（Wu et al.， 2017）.为避免重力沉降对颗粒接触力的长期干扰，本文将成样时间设置为运动时间的两倍，使试样内部充分平衡以达到准静态系统，并且在试样内部沿流动方向均布4个监测断面，在每个断面内逐层统计平均孔隙率，确保不同分形维数下试样的初始状态保持一致性.

颗粒柱的挡板被垂直放置在距离左侧边界L₀的位置，并在试样堆积完成后移除挡板使颗粒柱坍塌.颗粒柱受重力作用向x方向扩散，并最终沉积稳定，形成一个长度为L_f的堆积体（虚线部分）.数值模型中的物理参数如表1所示.

本文流体域高度为60倍平均粒径，长度方向为250倍平均粒径，以保证计算结果不受边界效应影响.为了保证计算模型的准确性并提升计算效率，计算域内的流体网格尺寸在各个方向均保持在2倍平均粒径大小.CFD和DEM的时间步长分别为10^-5 s和10^-6 s，并且每10个DEM时间步后进行一次CFD-DEM耦合信息交换.

为了研究不同流态下颗粒柱的坍塌动力学行为，图3展示了所有模拟工况在（St， r）参数空间中的分布情况，并依据Courrech du Pont et al. （2003）提出的临界参数（St_c=10，r_c=4，Re_c=2.5）对相应流态进行划分.图中蓝色数据点代表D=1.0的试样，黄色数据点代表D=3.5的试样.本文流态的变化是通过改变流体的动力粘度以及颗粒大小的方式实现，这种方式已被证明是有效的（Jing et al.， 2019）.表2总结了各计算工况的平均粒径、流体动力黏度、St及其所对应的流态区间.

为了验证本文流固耦合计算模型的准确性，在利用DEM模型通过标定参数获得与Polanía et al.（2024）针对高宽比a=1的干燥颗粒柱坍塌试验达到相同的运动距离后，本文对Polanía et al.（2024）高宽比a=1的水下单分散颗粒柱坍塌试验进行了模拟.图4展示了3个典型时刻（T=t_f/3， 2t_f/3t_f； t_f为最终运动停止时刻）水下颗粒柱坍塌过程的沉积轮廓，并与Polanía et al.（2024）的试验结果在对应时刻的沉积形态进行对比.可以看出，在选取的3个关键时刻，通过数值模拟得到的沉积物轮廓高度以及前沿运动距离均与试验结果保持了较好的一致性，定量地再现了相关试验的观测结果.这种比较验证了本文的数值模型在捕捉浸没颗粒坍塌的瞬态流动过程中的准确性和有效性.

2 考虑粒径分布特征的颗粒柱坍塌模拟

2.1　多分散颗粒系统动力学过程

图5~8记录了多分散颗粒柱在4种不同流态下的几个典型时刻对应的运动形态，黑色虚线代表了初始颗粒柱的位置.可以看出，多分散颗粒柱的动力学过程在不同流态内表现出明显的差异.

在自由下落态中，当通过设置颗粒柱尺寸与平均粒径之比保持一致从而避免出现有限尺寸效应后，具有不同分形维数的颗粒柱在倒塌后，颗粒系统前端快速向前发展，颗粒流呈现明显的流态化特征，最终形成近似锥形的堆积体.并且两者运动过程和最终运动距离都没有明显的差别，与之前的研究结论相同（Cabrera and Estrada， 2021）.与自由下落态类似，多分散颗粒柱在惯性态中的动力学过程没有因为具有不同粒径特征而出现显著差异.这与在均匀粒径中得到的结论一致（Bougouin and Lacaze， 2018），证明在惯性态下流体对颗粒的影响较小，颗粒的运动过程主要受颗粒惯性力支配.并且由于流体阻力的存在，导致惯性态中颗粒的运动距离小于自由下落态（Tang et al.， 2024）.

在粘性-惯性过渡态中，分形维数越高的颗粒柱运动距离越短.表明随着流体粘性作用增强，流体-颗粒相互作用行为对运动演化过程的影响逐渐显现，粒径分布特征已成为影响颗粒柱运动能力的重要因素.进入粘性态后，坍塌的颗粒柱会因为受流体的明显阻力而移动缓慢，由于强烈的粘性耗散作用，颗粒柱的移动距离进一步降低.与其他流态相比，粘性态下不同分形维数的颗粒系统动力学差异更为显著，分形维数越高的系统流动性越差.

2.2　不同流态下流动性对比

图9比较了多分散颗粒柱在不同流态内流动前端随运动时间的演变曲线.由于颗粒柱高度不同，利用特征时间T=\sqrt[]{H_{\mathrm{0}}/\mathrm{g}}对颗粒柱的运动时间进行归一化处理，式中H₀是初始颗粒柱高度，重力加速度g根据干燥或浸没状态取相应值.从图中可以看出，在自由下落态和惯性态中，不同粒径分布的颗粒柱在运动过程和最终流出距离方面都展现出较高的一致性，此时粒径分布对颗粒系统的动力学过程几乎没有影响.在粘性-惯性过渡态中，不同粒径分布的颗粒系统在流动性方面开始展现出差异，高分形维数的颗粒柱运动能力整体偏低.而在粘性态中，这种差异更为显著，并且在运动开始，多分散颗粒柱能在一定时间内保持稳定.高分形维数的颗粒柱运动时间显著增长，并且运动过程也与低分形维数的颗粒柱出现差异，最终运动距离降低.

图10统计了不同工况下颗粒柱的无量纲运动距离与对应St的关系.从图中可以看出，随着St的降低，流态逐渐从自由下落态向粘性态转变，颗粒系统的流动性逐渐降低，并且这种关系呈现出良好的一致性.这一现象表明流体黏性对颗粒运动能力具有持续抑制作用.同时，在惯性-黏性过渡态及粘性态下（蓝色区域），粒径分布特征对最终流动距离的影响开始显著体现，分形维数较高（即细颗粒比例较大）的试样表现出更短的运动距离，说明多分散性在强流固耦合环境中可能通过改变系统能量耗散机制与渗透性从而显著影响了颗粒系统的宏观迁移能力.

3 粒径分布特征对流动特性的影响机制

3.1　细颗粒耗散机制

为了证明多分散性在流固耦合环境中可能通过改变能量耗散机制对颗粒系统的宏观迁移能力产生影响，图11a记录了2种分形维数的颗粒柱在不同流态下的势能转化（红线）、累计耗散能（黑线）、颗粒动能（黄线）以及流体动能（蓝线）随时间的变化情况.所有的能量均按照初始势能E₀进行归一化处理，图中实心点代表D=1.0的颗粒柱，空心点代表D=3.5的颗粒柱.对于自由下落的干燥状态，约有50%的初始势能（红线）参与到了向动能转换并最终耗散的过程，并且随着流态从自由下落态向粘性态转换，这种参与转换的能量占比逐渐降低.在转换的能量中，大部分的能量都被耗散（黑线），仅有小部分转换成了动能（黄线、蓝线）.并且在粘性态时，黑线和红线几乎重合，表明几乎所有的能量都在运动过程中被耗散，导致颗粒和流体获得的动能最少，运动能力最差.

比较自由下落态和惯性态的能量曲线可以发现，不同粒径分布特征的颗粒柱在能量转换形式上表现出一致性，这种在图中表现为实心点和空心点随时间变化的曲线几乎平行.因此两者的颗粒动能演化曲线也几乎重合，表现为在自由下落态和惯性态下，不同粒径分布的颗粒柱的运动能力没有显著差异.而在过渡态和粘性态时，不同粒径分布特征的颗粒柱对应的能量转换曲线和能量耗散曲线存在明显区别：细颗粒含量更多的高分形维数颗粒柱（空心）参与转换的能量占比更低，并且转换的能量几乎都转变成了能量耗散，证明在流体粘度更大的情况下，细颗粒的存在将显著增加颗粒系统的能量耗散，从而大大降低流动性.图11b统计了不同分形维数试样在各流态下的能量耗散比例.结果表明，在自由下落态和惯性态中，D=1.0和D=3.5的试样之间能量耗散差异较小，且整体耗散比例较低，系统能够获得更多动能，从而表现出较强的流动性.而在惯性-粘性过渡态和粘性态中，能量耗散占比显著增加，且不同分形维数试样之间的差异开始显现：D=3.5试样中由于细颗粒含量更高，导致耗散能量显著增加，流动能力相应减弱.这种现象与Shi et al.（2024）在干颗粒流中观测到的结果一致.

3.2　渗透性差异机制

上一节从能量角度分析了具有不同粒径分布特征的颗粒柱存在流动性差异的原因.但是能量演化曲线的不一致仅仅是流动性差异的另一种表现形式，不同分形维数引入了哪些物理机制造成了能量演化过程的区别从而导致了流动性出现明显差异值得去进一步探索.在研究体积分数对浸没颗粒系统动力学行为的影响时，渗透性被广泛认为是其流动性的关键因素（Bougouin and Lacaze， 2018）.因此，本节利用渗透性测试方法，对不同分形维数的颗粒系统进行渗透性评估，以进一步揭示粒径分布特征如何通过改变渗透性来影响浸没颗粒系统的流动行为.

流体在颗粒系统中运动时，流固耦合作用会产生流动阻力，导致沿流动方向出现压力降，压力降的程度直接反映系统的渗透性特性.近年来，Ergun方程被广泛用于预测颗粒系统中的压力降（Ergun， 1952； Gao et al.， 2024）.对于不可压缩流体通过均一粒径系统的情况，沿程压力降∆P/∆L可表示为：

\frac{∆P}{∆L}=C_{\mathrm{1}}\frac{{\mu }_{\mathrm{f}}}{d_{\mathrm{p}}^{\mathrm{2}}}\frac{{\left(\mathrm{1}-\epsilon \right)}^{\mathrm{2}}}{{ϵ}^{\mathrm{3}}}{\mathit{U}}_s+C_{\mathrm{2}}\frac{{\rho }_{\mathrm{f}}}{d_{\mathrm{p}}}\frac{\left(\mathrm{1}-\epsilon \right)}{{ϵ}^{\mathrm{3}}}{\mathit{U}}_{\mathrm{s}}^{\mathrm{2}}，(9)

式中：{\mu }_{\mathrm{f}}是流体动力粘度，d_{\mathrm{p}}是颗粒平均粒径，\epsilon是孔隙率，{\mathit{U}}_{\mathrm{s}}是流体速度，C₁=150，C₂=1.75是基于实测值给定的经验系数.然而，Ergun方程的压力降预测仅适用于均匀粒径系统.尽管已有研究针对双组分及多组分颗粒体系提出了修正的Ergun公式（Gao et al.， 2024），但对于分形维数控制的多分散颗粒系统，其适用性尚未得到充分验证.因此，仍需通过数值模拟方法，系统研究不同粒径分布特征对颗粒柱渗透性的影响.

本节通过固定初始试样内的颗粒位置，使水流沿x方向（颗粒柱坍塌的主要流动方向）通过不同粒径分布的颗粒柱，进而测定其渗透特性.图12展示了不同粒径分布特征的颗粒柱在各流态下压力降随流速的变化规律.结果表明，颗粒系统沿程压降在从惯性态向粘性态的过渡过程中而增大，反映出流体黏性在不同流态下对系统阻力的调控作用.在相同流体粘度条件下，低分形维数的颗粒系统始终表现出更高的压降增长率，且该差异在从惯性态向粘性态过渡的过程中愈加显著，表明粒径分布特征在高黏性环境中对流体-颗粒耦合阻力的影响更加敏感.

为进一步定量评估多分散颗粒系统的渗透性，通过计算图12所示沿程压降曲线的平均斜率，并将其定义为渗透性的表征参数ζ.图13比较了D=1.0和D=3.5的两种颗粒柱在不同流态下的渗透性.从图中可以看出，渗透性与分形维数之间具有良好的对应关系，分形维数越高的颗粒柱渗透性越差，最终的运动距离越短.分形维数越高其内部孔隙结构将更为复杂，而通过渗透性测试则可利用沿程压力降以及渗透性等指标定量地反映出这些孔隙结构的复杂程度.这种渗透性的降低将直接限制孔隙水从复杂孔隙结构中快速排出，进而加剧流体-颗粒之间的耦合阻力和能量耗散，从而削弱颗粒系统的运动能力.综上所述，在相同的体积分数和流态条件下，不同粒径特征的渗透系数呈现明显差异，并进一步引起宏观流动特性的明显区别.这一现象证明，与体积分数相比，渗透性能够更本质地表征水下颗粒系统与环境流体相互作用的动力学特征，是描述此类多相体系流体动力学行为的核心物理参量.

3.3　基于渗透性修正的流动能力统一描述方程

以上研究结果表明，渗透性作为关键物理参数，不仅建立了粒径分布特征与宏观动力学行为之间的关系，同时定量表征了由粒径分布差异所导致孔隙结构的复杂程度.因此，图14提出了一个考虑渗透性修正的St与无量纲流动距离的标度律，以反映粒径级配特征在不同流态下对流动性的影响.从图中可以看出，利用计算的渗透性指标ζ对St值进行修正后，所有数据与运动距离表现出良好的一致性.为了证明这种标度律的有效性，本文进一步考察了自然堆积条件下D=1.0和D=3.5的两种颗粒柱在不同流态中的运动距离.数据显示，自然堆积条件下的流动性数据依然在不同流态下表现出与修正后的St具有良好的一致性.这一结果表明，渗透系数作为反映多分散颗粒系统整体渗透能力的初始参数，能够在均匀流场条件下较为可靠地预测浸没颗粒系统的流动性.然而本文中渗透性参数的测定采用的是整体平均尺度，未对塌缩过程中颗粒柱内部可能出现的空间异质性特征进行定量分析.在实际坍塌过程中，颗粒系统内部的孔隙结构与渗透性可能存在显著的空间非均匀性，这可能进一步影响局部孔压耗散行为及流动模式.因此，在未来研究中，有必要结合局部渗透率场的提取与演化分析，进一步揭示空间异质性对系统整体流动性的调控机制，特别是在非均匀流场或复杂地形条件下的适用性验证仍需深入探讨.

4 结论

本文采用CFD-DEM方法对具有不同粒径分布特征的颗粒柱在不同流态下的坍塌过程进行模拟，系统展示了颗粒系统在全流态下的运动特征，并进一步分析了不同流态下多分散颗粒柱流动性出现差异的原因，具体结论如下：

（1）流态转变对颗粒柱坍塌过程具有显著影响：随着流态由自由下落态逐渐向粘性态演化，颗粒柱的坍塌距离逐渐减小，流动持续时间延长，运动模式也发生明显变化.进一步分析表明，在高粘度流体环境下，细颗粒含量较高的颗粒系统更容易发生能量耗散，从而显著削弱整体流动性.

（2）颗粒分布特征对流动性的影响具有明显的流态依赖性：在流体粘度较低的条件下（即自由下落态与惯性态），粒径分布对颗粒系统流动性的影响较弱；而在流体粘度较高的条件下（惯性-粘性过渡态及粘性态），分形维数越高的颗粒系统，其运动能力越显著降低.

（3）体积分数作为反映颗粒系统平均密实度的物理量，在描述水下多分散颗粒系统的动力学行为时存在一定局限性，特别是在面对具有复杂孔隙结构的系统时.相比之下，渗透性不仅能够反映颗粒系统的整体密实程度，还能表征其内部孔隙结构的复杂性，并可直接影响孔压耗散过程.因此，渗透性相较于体积分数更适合作为控制水下多分散颗粒系统流动性的关键参数.研究结果表明，渗透性越差，系统的运动能力越弱.

（4）基于渗透性与运动能力之间的耦合关系，本文提出了一个新的流动距离标度律，能够有效量化多分散颗粒系统在不同流态下的运动能力.该标度律不仅具有明确的物理意义，也可为地质灾害风险评估中的滑动距离预测提供理论指导与模型支持，进一步加深了对库岸滑坡、海底滑坡等链式地质灾害中高流动性本质的认知.

参考文献

[1] Bougouin, A., Lacaze, L., 2018. Granular Collapse in a Fluid: Different Flow Regimes for an Initially Dense-Packing. Physical Review Fluids, 3(6): 064305. https://doi.org/10.1103/physrevfluids.3.064305

[2] Breard, E. C. P., Dufek, J., Charbonnier, S., et al., 2023. The Fragmentation-Induced Fluidisation of Pyroclastic Density Currents. Nature Communications, 14: 2079. https://doi.org/10.1038/s41467-023-37867-1

[3] Cabrera, M., Estrada, N., 2021. Is the Grain Size Distribution a Key Parameter for Explaining the Long Runout of Granular Avalanches? Journal of Geophysical Research: Solid Earth, 126(9): e2021JB022589. https://doi.org/10.1029/2021JB022589

[4] Courrech du Pont, S., Gondret, P., Perrin, B., et al., 2003. Granular Avalanches in Fluids. Physical Review Letters, 90(4): 044301. https://doi.org/10.1103/physrevlett.90.044301

[5] Crosta, G. B., Hermanns, R. L., Dehls, J., et al., 2017. Rock Avalanches Clusters along the Northern Chile Coastal Scarp. Geomorphology, 289: 27-43. https://doi.org/10.1016/j.geomorph.2016.11.024

[6] Cui, W., Wei, J., Wang, C., et al., 2021. Discrete Element Simulation of Collapse Characteristics of Particle Column Considering Gradation and Shape. Chinese Journal of Geotechnical Engineering, 43(12): 2230-2239 (in Chinese with English abstract).

[7] Cundall, P. A., Strack, O. D. L., 1979. A Discrete Numerical Model for Granular Assemblies. Géotechnique, 29(1): 47-65. https://doi.org/10.1680/geot.1979.29.1.47

[8] Di Felice, R., 1994. The Voidage Function for Fluid-Particle Interaction Systems. International Journal of Multiphase Flow, 20(1): 153-159. https://doi.org/10.1016/0301-9322(94)90011-6

[9] Ergun, S., 1952. Fluid Flow through Packed Columns. Chemical Engineering Progress, 48 (2):89.

[10] Gao, S., Theuerkauf, J., Pakseresht, P., et al., 2024. A Modified Ergun Equation for Application in Packed Beds with Bidisperse and Polydisperse Spherical Particles. Powder Technology, 445: 120035. https://doi.org/10.1016/j.powtec.2024.120035

[11] Hsu, S. K., Kuo, J., Lo, C. L., et al., 2008. Turbidity Currents, Submarine Landslides and the 2006 Pingtung Earthquake off SW Taiwan. Terrestrial, Atmospheric and Oceanic Sciences, 19(6): 767. https://doi.org/10.3319/tao.2008.19.6.767(pt)

[12] Hyslip, J. P., Vallejo, L. E., 1997. Fractal Analysis of the Roughness and Size Distribution of Granular Materials. Engineering Geology, 48(3/4): 231-244. https://doi.org/10.1016/S0013-7952(97)00046-X

[13] Iverson, R. M., 2005. Regulation of Landslide Motion by Dilatancy and Pore Pressure Feedback. Journal of Geophysical Research: Earth Surface, 110(F2). https://doi.org/10.1029/2004JF000268

[14] Iverson, R.M., 1997. The Physics of Debris Flows. Reviews of Geophysics, 35(3): 245-296.

[15] Jing, L., Guo, S.Y., Zhao, T., 2019. Understanding Dynamics of Submarine Landslide with Coupled CFD-DEM. Rock and Soil Mechanics, 40(1): 388-394 (in Chinese with English abstract).

[16] Jing, L., Kwok, C. Y., Leung, Y. F., et al., 2016. Extended CFD-DEM for Free-Surface Flow with Multi-Size Granules. International Journal for Numerical and Analytical Methods in Geomechanics, 40(1): 62-79. https://doi.org/10.1002/nag.2387

[17] Jing, L., Yang, G. C., Kwok, C. Y., et al., 2019. Flow Regimes and Dynamic Similarity of Immersed Granular Collapse: A CFD-DEM Investigation. Powder Technology, 345: 532-543. https://doi.org/10.1016/j.powtec.2019.01.029

[18] Lacaze, L., Bouteloup, J., Fry, B., et al., 2021. Immersed Granular Collapse: From Viscous to Free-Fall Unsteady Granular Flows. Journal of Fluid Mechanics, 912: A15. https://doi.org/10.1017/jfm.2020.1088

[19] Lai, Z. Q., Chen, D., Jiang, E. H., et al., 2021. Effect of Fractal Particle Size Distribution on the Mobility of Dry Granular Flows. AIP Advances, 11(9): 095113. https://doi.org/10.1063/5.0065051

[20] Lee, C. H., 2021. Two-Phase Modelling of Submarine Granular Flows with Shear-Induced Volume Change and Pore-Pressure Feedback. Journal of Fluid Mechanics, 907: A31. https://doi.org/10.1017/jfm.2020.838

[21] Li, K., Cheng, Q. G., Lin, Q. W., et al., 2022. State of the Art on Rock Avalanche Dynamics from Granular Flow Mechanics. Earth Science, 47(3): 893-912 (in Chinese with English abstract).

[22] Lube, G., Huppert, H. E., Sparks, R. S. J., et al., 2004. Axisymmetric Collapses of Granular Columns. Journal of Fluid Mechanics, 508: 175-199. https://doi.org/10.1017/S0022112004009036

[23] Man, T., Ge, Z., Huppert, H.E., et al., 2022. Transient Rheology and Size Effect in Granular Column Collapses. Chinese Journal of Computational Mechanics, 39(3): 381-388 (in Chinese with English abstract).

[24] Nian, T. K., Shen, Y. Q., Zheng, D. F., et al., 2021. Research Advances on the Chain Disasters of Submarine Landslides. Journal of Engineering Geology, 29(6): 1657-1675 (in Chinese with English abstract).

[25] Oquendo-Patiño, W. F., Estrada, N., 2022. Finding the Grain Size Distribution That Produces the Densest Arrangement in Frictional Sphere Packings: Revisiting and Rediscovering the Century-Old Fuller and Thompson Distribution. Physical Review E, 105(6): 064901. https://doi.org/10.1103/PhysRevE.105.064901

[26] Pailha, M., Nicolas, M., Pouliquen, O., 2008. Initiation of Underwater Granular Avalanches: Influence of the Initial Volume Fraction. Physics of Fluids, 20(11): 111701. https://doi.org/10.1063/1.3013896

[27] Peng, J. B., Zhang, Y. S., Huang, D., et al., 2023. Interaction Disaster Effects of the Tectonic Deformation Sphere, Rock Mass Loosening Sphere, Surface Freeze-Thaw Sphere and Engineering Disturbance Sphere on the Tibetan Plateau. Earth Science, 48(8): 3099-3114 (in Chinese with English abstract).

[28] Peng, M., Zhao, Q. X., Li, S., et al., 2025. Two-Phase SPH Simulation of Granular Landslide-Tsunamis Processes Considering Dynamic Seepage. Earth Science: 1-13(in Chinese with English abstract).

[29] Polanía, O., Cabrera, M., Renouf, M., et al., 2022. Collapse of Dry and Immersed Polydisperse Granular Columns: A Unified Runout Description. Physical Review Fluids, 7(8): 084304. https://doi.org/10.1103/physrevfluids.7.084304

[30] Polanía, O., Estrada, N., Azéma, E., et al., 2024. Polydispersity Effect on Dry and Immersed Granular Collapses: An Experimental Study. Journal of Fluid Mechanics, 983: A40. https://doi.org/10.1017/jfm.2024.176

[31] Rondon, L., Pouliquen, O., Aussillous, P., 2011. Granular Collapse in a Fluid: Role of the Initial Volume Fraction. Physics of Fluids, 23(7): 073301. https://doi.org/10.1063/1.3594200

[32] Shi, A. N., Yang, G. C., Kwok, C. Y., et al., 2024. Enhanced Mobility of Granular Avalanches with Fractal Particle Size Distributions: Insights from Discrete Element Analyses. Earth and Planetary Science Letters, 642: 118835. https://doi.org/10.1016/j.epsl.2024.118835

[33] Tang, H., Lin, B. S., Wang, D. M., 2024. Granular Collapse in Fluids: Dynamics and Flow Regime Identification. Particuology, 92: 30-41. https://doi.org/10.1016/j.partic.2024.04.013

[34] Tang, H., Wu, Y. S., Wang, D. M., 2025. Dynamics and Flow Regime of Immersed Granular Collapse: Role of Fractal Particle Size Distribution. Computers and Geotechnics, 186: 107392. https://doi.org/10.1016/j.compgeo.2025.107392

[35] Trulsson, M., Andreotti, B., Claudin, P., 2012. Transition from the Viscous to Inertial Regime in Dense Suspensions. Physical Review Letters, 109(11): 118305. https://doi.org/10.1103/physrevlett.109.118305

[36] Wang, C., Wang, Y. Q., Peng, C., et al., 2017. Dilatancy and Compaction Effects on the Submerged Granular Column Collapse. Physics of Fluids, 29(10): 103307. https://doi.org/10.1063/1.4986502

[37] Wu, S. G., Sun, Y. B., Li, Q. P., 2019. Deepwater Geohazards in the South China Sea. Science Press, Beijing (in Chinese).

[38] Wu, Y. L., An, X. Z., Yu, A. B., 2017. DEM Simulation of Cubical Particle Packing under Mechanical Vibration. Powder Technology, 314: 89-101. https://doi.org/10.1016/j.powtec.2016.09.029

[39] Yang, G. C., Jing, L., Kwok, C. Y., et al., 2020. Pore-Scale Simulation of Immersed Granular Collapse: Implications to Submarine Landslides. Journal of Geophysical Research: Earth Surface, 125(1): e2019JF005044. https://doi.org/10.1029/2019JF005044

[40] Yin, Z. Y., Xu, Q., Hu, W., 2012. Constitutive Relations for Granular Materials Considering Particle Crushing: Review and Development. Chinese Journal of Geotechnical Engineering, 34(12): 2170-2180 (in Chinese with English abstract).

基金资助

国家自然科学基金-联合基金重点项目(U23A2044)

国家自然科学基金-联合基金重点项目(U24A20170)

广西重点研发计划项目(No.桂科AB25069121)

国家自然科学基金-面上项目(12472412)

广东省基础与应用基础研究基金项目(2022B1515130006)