摘 要
驱动粘性岩土材料普遍表面侵蚀的基本机制可以从微观力学角度结合著名的数值技术进行分析。在这里，结合离散元和点阵玻尔兹曼方法的耦合模型(DEM-LBM)提供了在垂直流体射流撞击颗粒组件表面引起的瞬态侵蚀过程中固体-流体相互作用的深入见解。本文采用单参数屈服曲面的简单键合模型来描述固体颗粒间的黏聚性。地表侵蚀的初始拓扑结构倾向于模拟直接在土壤表面上方的流体速度剖面。我们发现，侵蚀速率和侵蚀质量的大小直接取决于单个固体键的微观力学强度。
背 景
土壤表面可蚀性的地质力学评估通常是根据经验进行的，通常是通过标准的实验试验，如喷射侵蚀试验和孔侵蚀试验。这些测试手段的对象一般是有粘性的岩土材料，而在其表面流动的水起着腐蚀剂的作用。
可蚀性通常是根据土壤表面相对于起始侵蚀的阈值(临界剪应力)的平均水动力剪应力来量化的。然而，在实践中，侵蚀过程要复杂得多，并涉及局部瞬态流动，如带有局部水动力应力间歇峰值的湍流射流。另一方面，由于与粘性土的颗粒间力及其复杂关系有关的泥沙性质种类繁多，导致临界应力的可靠经验表达式的推导受到阻碍。所有这些导致了对不同测试程序的稳健性以及分析方法的普遍缺乏共识，特别是在评估临界剪切应力的假设方面。
从理论方面来看，迄今为止用于研究岩土材料侵蚀的大多数数值方法的出发点仅仅是宏观力学，将水和土壤描述为连续介质。因此，人们提出了不同的技术来跟踪它们界面的瞬态演化，其中相之间的质量交换(即实际侵蚀)发生。一方面，纯欧拉方法允许使用固定网格，但不允许在界面上精确定义流体力学变量。类似地，混合欧拉-拉格朗日方法可能在界面处显示出良好的精度，但它们通常会导致严重和耗时的重新网格问题。
本研究的目的是介绍一种基于离散元模拟的方法，即采用格点玻尔兹曼方法模拟流体射流对粘性颗粒床的影响。这是一种微观力学方法，它允许研究驱动侵蚀基本机制的不同因素，例如土壤表面颗粒团聚体的分离(表面侵蚀)或自由表面以下物质块的移动(质量侵蚀)。
下面我们将介绍我们的数值方法及其在颗粒间黏结不断增加的颗粒床上的JET测试中的应用。将考虑粒子间凝聚力的简单模型，并将讨论其对床层可蚀性的后果。
数值方法
目前已有多种耦合方法用于颗粒尺度流固相互作用的数值模拟。本文采用的方法结合了用于模拟流体的格子玻尔兹曼方法(LBM)的计算效率和用于模拟固体颗粒的离散元方法(DEM)的相对简单性，这种方法已经应用于几个地质力学问题。
固 相
本文采用基于分子动力学方法的光滑离散元方法对颗粒土进行了描述。我们使用圆形粒子的二维组合，其轨迹是通过牛顿运动方程的简单积分来计算的。粒子间的光滑接触和摩擦规律由内聚相互作用补充。这里，对于接触关系我们考虑图1中描述的一个简单的内聚定律。
图1：左：固体接触流变模型；右图：作用力空间内粘性键的屈服面
因此，两个粒子之间的相互作用可以用施加在接触点上的相互作用力F和相互作用力矩M来表示。将作用力分解为其法向分量Fn和剪切分量Fs，分别对应于接触处的法向方向和切向方向。法向力Fn用局部互渗δn描述，用经典的Kelvin-Voigt粘弹性关系定义法向刚度和阻尼系数kn和ηn。另一方面，在摩擦(非粘性)接触处产生的剪切力在这里由粘性正则化库仑定律定义，其特征是静摩擦系数μ和剪切正则化粘性系数ks。这样，剪力被计算为滑动速度 s的函数，而滑动速度在实践中从未真正消失。由于本文考虑的时间步长较小(10~6s)，在静摩擦力作用下，颗粒通过残余微滑移达到准静力平衡。作用于每个颗粒中心的力矩由以颗粒半径为杠杆臂的剪切力和控制颗粒滚动运动的滚动摩擦分量定义(相对旋转而不滑动)。后者类似于滑动摩擦，用相对旋转速度代替滑动速度，用滚动摩擦系数μω和正则化系数kω代替滑动摩擦系数及其正则化系数。
到目前为止所介绍的成分描述了纯摩擦接触。为了引入允许固体颗粒之间的拉伸法向力的内聚力，在接触力空间中采用了一个具有抛物线屈服面的简单弹塑性模型(图1)。通过这种方式，为所有最初接触的颗粒定义了固体键，其特征是线弹性流变，其特征是法向键和剪切键刚度kn,b和ks,b。当接触力到达屈服面时，粘结断裂，接触变成纯摩擦。为方便起见，假定正态屈服阈值和剪切屈服阈值Cn和Cs仅依赖于一个参数C=Cn= 2·Cs，该参数表示固体粘结的强度或胶结程度。多分散组合中键的相对强度可以用一个无量纲数来表征，该数被定义为键内聚力与颗粒自身浮力重量的比率。这个内聚键数可以表示为：
其中Δρ是固体颗粒的淹没表观密度，g是重力加速度，V是颗粒的体积。一般来说，可以注意到B≥0和B < 1的固体键是不稳定和短暂的，因为在其自身重量下，任何轻微的组合重排都有可能导致键断裂。我们的数值试验表明，当B≥3时，样品中的所有颗粒在重力作用下保持完全粘结。
液 相
本文用格Boltzman方法的经典D2Q9格式描述了流体相的瞬态流动，该格式是连续Boltzmann方程的显式有限差分格式。LBM基于空间中的笛卡尔网格和每个网格点上流体粒子的概率密度函数和速度的离散集。然后，流体动力学计算只需两个基本步骤，即流体粒子的碰撞和平流，遵循一套确保质量、动量和能量守恒的规则(即不可压缩的纳维-斯托克斯行为)。然而，这只适用于低马赫数，即当流体速度远低于晶格速度。为了克服经典Bhatnaggar-GrossKrook方案(单松弛时间)的众所周知的缺点，我们使用了具有多重松弛时间(MRT)的d 'Humières的广义公式。因此，流体相可以通过指定晶格网格大小Δx和晶格速度cs，对角弛豫矩阵S的分量Sα(不同弛豫时间的倒数)以及密度ρf和运动粘度ν的流体材料参数来完全定义。
固液耦合
流体与固体相之间的耦合以及对颗粒的水动力计算是通过Bouzidi等提出的简单动量交换方法来实现的，该方法依赖于LBM中任意几何形状的广义反弹和线性插值策略。Feng等提出的子循环时间积分技术，利用流体计算所需的时间步长通常比粒子动力学大得多的事实。然而，为了保持固体颗粒上计算的水动力的准确性，每个LBM步长的DEM子循环数被限制为2。为了解释流体在二维颗粒密集组合中的流动，使用了略小于颗粒半径的“液压”半径。
结论
采用DEM-LBM耦合模型研究了流体射流作用下粘性岩土材料表面侵蚀的细观力学特性。我们表明，对于给定的射流，存在一个发生表面侵蚀的粒子间凝聚的临界值。我们目前正在研究更复杂的内聚力定律的影响以及微观力学和宏观屈服阈值之间的关系。
,
