颗粒柱塌落中的尺寸效应和瞬态流变性研究

满 腾， 葛 转， Herbert E.Huppert, Sergio Andres Galindo Torres*

(1.西湖大学 工学院 浙江省海岸带环境与资源研究重点实验室，杭州 310024；2.剑桥大学 国王学院，英国剑桥 CB2 1ST)

1 引 言

对于颗粒材料的瞬态流变性研究，颗粒材料柱体塌落行为常作为理想的研究对象，其运动规律虽体现瞬态动力学特性却较为简单，同时颗粒柱的塌落行为也可与更为复杂的滑坡和泥石流等自然颗粒流建立联系，方便进一步将小尺度的柱体塌落问题转化为大尺度的自然地球物理流的研究[12]。文献[13-15]通过不同颗粒材料(沙粒、大米颗粒和小米颗粒等)塌落试验研究以及量纲分析总结出颗粒材料柱体塌落扩散程度随颗粒柱初始高径比的变化趋势，认为当柱体的初始高径比α小于某个临界值αc时，无量纲扩散比R*=(R∞-Ri)/Ri与初始高径比α=Hi/Ri成线性关系;而当α>αc时，R*与初始高径比α0.5成线性关系，其中R∞为柱体的最终塌落范围的半径，Ri为颗粒柱截面的平均初始半径。此相关关系只针对二维颗粒柱塌落以及圆形截面柱体塌落的情况。

Staron等[16,17]在此基础上研究了颗粒间接触特性对颗粒柱整体塌落行为的影响，认为颗粒之间的摩擦系数对柱体的塌落扩散程度有较大的影响。前期研究通过改变颗粒材料的粒间摩擦系数、颗粒-底面摩擦系数、柱体系统相对尺寸以及柱体截面形状等研究了相关参数对颗粒柱塌落的影响。首先，本文认为颗粒柱的塌落扩散比可以通过等效柱体高径比αeff进行描述[18]，

(1)

式中μw为颗粒与底面的摩擦系数，μp为颗粒之间的摩擦系数，β为常数,用以表征两种摩擦之间的能耗比。其次，颗粒柱塌落的尺寸效应可以用有限尺寸标度模型(2)表述[19]，

R*=(Ri/d)- β1/νFr[(αeff-αc ∞)(Ri/d)1/ν]

(2)

式中Fr[·]为尺度标度函数，d为颗粒平均粒径，β1=0.28±0.04和ν=1.39±0.14为常数，且与二维随机网络逾渗问题的有限尺寸标度模型参数相似。

本文在上述研究的基础上总结颗粒材料塌落扩散距离的计算方法，并通过对颗粒材料柱体塌落过程的应力和应变率的分析，研究颗粒材料在瞬态过程的流变特性。

2 离散单元法及模型的建立

2.1 离散单元法和接触模型

本文采用扩展多面体离散单元法(DEM)。使用DEM可以方便获得颗粒材料柱体塌落过程中的整体动力学响应，特别是能够获得颗粒之间的相互作用力，便于进行细观尺度上的力学分析，也有利于从颗粒尺度对问题进行深入研究。多面体单元由若干平面组成，其几何构成主要为角点、棱边和平面，在多面体单元的基础上采用球体和任意多面体的Minkowski Sum来得到扩展多面体[20,21]，这样颗粒之间相互作用的求解转变为两个多面体的扩展球体之间的相互作用。

扩展多面体DEM的颗粒接触有角点与角点的接触、棱边与棱边的接触以及角点与平面的接触三种形式。在计算过程中，首先根据颗粒的位置以及其几何形状确定两个颗粒的接触变形δn，利用线性粘弹性模型求得颗粒之间的接触面法向相互作用力。而颗粒之间的切向相互作用则需要在每一时间步计算两个颗粒接触点位置的切向位移，并确保利用线性模型计算颗粒之间的切向相互作用并且考虑颗粒之间的切向作用不得大于颗粒之间的滑动摩擦力。根据颗粒之间的相互作用力，可以实时利用牛顿第二定律计算每一个颗粒的加速度、速度和位移，获取颗粒之间的相互作用力信息。本文DEM计算参数为正向刚度系数Kn=1×104N/cm,切向刚度系数Kt=5×103N/cm,颗粒碰撞回弹比e=0.1。本文采用多物理场仿真程序库MechSys的DEM模块进行计算和研究，该数值计算方法已在文献[18-20]中得到多次验证，可以进行颗粒材料柱体塌落的研究。

2.2 模型的建立

本文对两种不同截面形状(圆形截面和长方形截面)的颗粒材料柱进行研究，利用MechSys在设计区域内利用三维Voronoi图进行分网，建立不规则扩展多面体颗粒模型，此方法生成的Voronoi颗粒的尺寸遵循高斯分布[20]，在此基础上随机选取20%的颗粒删除，使得初始的堆积柱体具有80%的固体颗粒体积分数，且孔隙分布较为均匀，此方法可减少颗粒流动过程中的剪胀效应，同时增加系统的随机性。此后，移除限制颗粒材料运动的边界条件(如颗粒柱截面为圆形，则此边界条件为圆管)使颗粒材料在自重条件下自然塌落。待颗粒停止运动，测量其最终扩散距离。对于圆形截面颗粒柱的塌落，本文测量四个方向的扩散距离并求平均，对于其他非圆形截面颗粒柱，则选取特定方向的扩散距离予以测量，如对于长方形截面颗粒柱，仅测量对角线方向的颗粒扩散距离以及垂直于截面棱边中点方向的平均扩散距离。长方形截面长边 6 cm，短边3 cm，其初始高度从1 cm增加至 50 cm。为研究尺寸效应，圆形截面柱的直径选取Di=1.0 cm,3.0 cm和7.0 cm。同时，颗粒之间的摩擦系数为0.2,0.4和0.6，而颗粒与底面的摩擦系数保持为0.4不变。

图1 不同截面形状颗粒柱的塌落情况以及塌落后的颗粒数量沿不同方向的分布情况

3 结果与讨论

3.1 颗粒材料柱体塌落及其尺寸效应

首先考察颗粒之间的摩擦系数对颗粒柱塌落的影响。文献[18]给出了不同摩擦系数条件下塌落扩散比R*与初始高径比α之间的关系，增加颗粒间摩擦系数可以显著降低其塌落扩散的距离。对于同一粒间摩擦系数的情况，虽然塌落扩散比R*随初始高径比的增加而单调递增，但是其在双对数坐标系下的斜率在某一临界高径比αc前后发生改变，一般认为这一临界点表征着颗粒材料具备不同的动力学特性。

当颗粒材料柱的初始高径比小于临界高径比时，颗粒材料处于拟静态;当大于临界高径比时,颗粒材料处于惯性流状态或近似液态。更为重要的是，改变颗粒之间的摩擦系数不仅能够改变其塌落扩散距离，也能改变临界高径比的大小。根据前期研究，当横坐标为有效高径比αeff时(如式(2))，本文得到的R*塌缩到一条曲线，从而在计算塌落扩散比时得以考虑摩擦系数的影响。本文进一步考察了两种具备不同摩擦系数的颗粒混合在一起时的颗粒柱塌落扩散行为，此时颗粒与底面的摩擦系数为μw=0.4，颗粒1#之间的摩擦系数为μp 1=0.1，颗粒2#之间的摩擦系数为μp 2=0.6，颗粒1#和颗粒2#之间的摩擦系数通过μ12=2μp 1μp 2/(μp 1+μp 2)计算。图2给出了两种不同配比情况的颗粒柱塌落扩散比与初始高径比的关系。当颗粒2#所占比例从10%增加到50%时，同样的初始堆积柱体获得的塌落扩散比显著增加，因而，即使是颗粒混合的情况，增加高摩擦系数的颗粒也可以显著提高颗粒材料整体的摩擦效应。

图2 两种不同颗粒配比条件下颗粒柱塌落扩散比R*与初始高径比α的关系，其中颗粒2#所占比例分别为10%和50%

本文同时考虑系统尺寸的大小对颗粒材料柱塌落扩散距离的影响。考虑三种不同的圆形截面颗粒柱尺寸，其初始直径分别为1 cm,3 cm和 7 cm，颗粒尺寸约为0.2 cm，所以相对系统尺寸Ri/d=2.5,7.5和17.5。当Ri/d=2.5时，临界有效高径比出现在αeff≈7.5；当Ri/d=7.5时，临界有效高径比出现在αeff≈3.0；当Ri/d=17.5时，临界有效高径比出现在αeff≈1.7。因而可以认为改变系统的相对尺寸，其临界高径比也随之变化。如图3所示，增加系统尺寸可以增加颗粒材料柱的相对塌落扩散距离，然而根据Cabrera等[22]的研究，这种扩散距离随系统尺寸增加而增加的趋势随着相对尺寸Ri/d>50而逐渐消失，前期研究指出，扩散距离随系统尺寸增加趋势的消失具备幂律衰减的特性。根据颗粒柱塌落中表现出的幂律衰减规律，本文在前期研究中利用有限尺度标度理论去分析模拟数据，从而得到如式(3)所示的尺寸效应规律。其中，ν=1.39±0.14为临界指数，αc ∞=0.2为系统尺寸无穷大时的临界有效初始高径比。而纵坐标为R*(Ri/d)β / ν，其中，β=0.28±0.04为临界指数且与常见材料的实验测量临界指数相似。在考虑有限尺度标度理论的情况下，所有计算结果均分布在一条规律的曲线上，从而在考虑颗粒间摩擦系数的基础上将系统相对尺寸也考虑在内。

图3 不同相对系统尺寸对R*和αeff关系的影响

3.2 柱体截面形状的影响

上述研究均以圆形截面颗粒柱的塌落为切入点，然而，非圆形截面颗粒柱往往表现出更为复杂的动力学特性。本文选取截面尺寸为6 cm ×3 cm的长方形截面柱作为研究对象。重点考察垂直于两种棱边方向以及沿对角线方向的颗粒柱塌落扩散情况。对于长方形截面柱，垂直于长棱边方向称为短棱向，垂直于短棱边方向称为长棱向，沿对角线方向仍称为对角向。

图4(a)为长方形截面柱DEM计算结果，其中，实心圆代表短棱向，空心圆代表长棱向，而×代表对角向计算结果。可见短棱向相对于其他方向具有更强的塌落扩散能力。需要明确的是，正方形截面的棱向以及长方形截面的短棱向均具有最小的初始半径Ri，因而该结果显示，较小的初始半径对塌落扩散有利。

然而，较小初始半径Ri对塌落扩散结果有利的结论与3.1节中尺寸效应的结果相反，因而重新定义无量纲的塌落扩散比为

(3)

图4 长方形截面颗粒柱不同测量方向对和αeff关系的影响

在此基础上，进一步考虑利用现有尺寸效应工具，即有限尺度标度理论对截面形状影响进行分析，相关结果如图5所示。可以看出，尽管截面形状改变和测量方向不同，考虑有限尺度标度理论后，DEM计算结果均与式(3)的结果吻合，其尺寸效应的临界指数也与式(3)给出的参数一致。

图5 考虑尺寸效应后的不同截面形状对等效塌落扩散比的影响

3.3 颗粒柱塌落与地球物理流

塌落扩散距离这一指标虽然可以表征颗粒材料在自重条件下的输运能力，但只展现了塌落堆积结果的其中一个维度，难以用于全面描述其堆积形态，而堆积形态对研究和预测地球物理流等复杂流体的运动和堆积更为重要。在前期研究中发现，当颗粒柱的初始高径比较大时，颗粒材料的塌落具有拟液态的特性，其塌落更多地体现为上覆颗粒推动下层颗粒塌落及扩散。塌落后，在初始高径比较小的情况下其堆积形态呈现为圆台或圆锥体，但当初始高径比较大时，颗粒更多地堆积在中心以及靠近外围的区域，而使0.5R∞处出现颗粒较少的情况，其总体堆积形态可总结为如图6(a)所示的中心开花，环状分布的情况。

图6 正方形截面颗粒柱塌落后形成的高度云图及网络中获取的法属波利尼西亚Bora-Bora岛卫星图

图6(a)为大初始高径比正方形截面颗粒材料柱塌落堆积云图，该算例初始高径比为20，颗粒间摩擦系数为0.4，颗粒与地面摩擦系数为0.4，可以看出，大部分颗粒分布在两处，一部分在中心处呈圆锥体分布，另一部分呈圆环状分布。该分布与常见的火山喷发形成的岛礁的形态有一定的相似性。如图6(b)所示，法属波利尼西亚的Bora-Bora岛的形态与颗粒柱塌落堆积形态类似，该岛屿中心为一死火山，外围有高出水面的环状岛礁，与文献[19]指出的颗粒柱塌落与火山碎屑流等地球物理流之间的联系相似。

3.4 塌落过程中应力和应变率的计算

(i,j=1,2,3) (4)

3.5 塌落过程中的瞬态流变性

图7 颗粒柱塌落过程的流变特性

(5)

图8 颗粒塌落后过程中颗粒速度分布情况以及代表单元分网情况 (红色边界区域的大于0.6(近底板固相区域))

同时，本文研究了不同相对系统尺寸(Ri/d=15和30)的颗粒柱塌落过程的流变特性，两组颗粒柱的初始直径分别为6 cm和12 cm，初始高度分别为12 cm和24 cm，其初始高径比相同，颗粒平均粒径均为0.2 cm。由于系统的整体尺寸较大，颗粒柱较高，因此在颗粒柱塌落过程中产生了两部分固相区域，如图9所示，一部分为近底板部分流动速度较小的区域，该区域形成圆锥体状的低速区域，另一区域为靠近颗粒柱顶端的具有高运动速度部分，此部分虽然颗粒流速较大，但颗粒之间的相对速度较小，而在这两部分中间形成高速剪切带。大系统的这种特性为流变性研究提供了丰富的信息。本部分研究集中于相对系统尺寸对流动起来的颗粒特性的影响，因此选取远离底板区域进行统计分析以减少底板边界条件的影响。

图9 颗粒塌落后过程中颗粒速度分布情况以及代表单元固相系数分布(黑色虚线以上为不受底板影响塌落区域)

4 结 论

本文系统地分析了颗粒材料柱在不同条件下的塌落堆积情况，特别是分析了颗粒柱的相对系统尺寸和截面形状对塌落扩散距离的影响。研究表明，减少颗粒之间的摩擦系数以及颗粒与底板直接的摩擦系数可以有效增加颗粒材料柱塌落扩散的相对距离，借助于αeff可以得到较为普适的描述不同摩擦系数条件下颗粒材料柱塌落扩散距离的计算方法。同时，发现增加颗粒材料柱的相对系统尺寸可以增加其塌落堆积的相对扩散距离，因而具备显著的尺寸效应，对尺寸效应的研究和理解有助于将实验室尺度的颗粒材料柱塌落的研究推广至更为宏观尺度的工程应用。针对尺寸效应，本文利用统计物理中常用的有限尺度标度理论进行分析，不同系统尺寸条件下，因考虑有限尺度标度理论，所有计算结果均塌缩到同一曲线，具有较好的拟合效果。

另外，对不同截面形状的颗粒材料柱的塌落堆积形态也做了相应的总结，指出垂直于截面棱边方向颗粒更易塌落扩散至较远的距离。同时考虑等效塌落扩散比，不同方向上颗粒的塌落堆积亦符合有限尺度标度理论。前期研究表明，颗粒材料柱的塌落往往可以与真实地球物理流(泥石流和火山碎屑流等)相关联，大初始高径比条件下颗粒材料柱的堆积形态与某些火山喷发后物质堆积形态具有一定程度的相似性，该领域仍需进一步的研究。

本文也从细观力学的角度探究颗粒材料在瞬态条件下的流变性能。以往对颗粒材料流变性的研究往往集中于稳态条件下，颗粒柱的塌落过程则与稳态相去甚远，本文求得颗粒材料柱塌落过程中不同时刻不同位置的计算单元的局部应力-应变率关系，发现塌落过程中时刻的选取对于应力应变关系的影响并不明显，然而改变颗粒材料中的摩擦系数以及系统的尺寸则对于其瞬态流变性有较大的影响。下一步将重点研究颗粒材料瞬态流变性量化分析，以期得到可以应用的瞬态本构关系方程。