1、第40 卷第2 期2023年6 月长江工程职业技术学院学报Journal of Changjiang Institute of TechnologyVol.40 No.2Jun.2 0 2 3SPH-DEM流固耦合方法模拟金沙江白格滑坡应用研究唐岳灏1，姜清辉?（1长江工程职业技术学院，武汉430 2 12；2 武汉大学，武汉430 0 7 2）摘要：从力学本质来讲，滑坡属于滑坡动力学和水动力学复杂流固耦合问题。采用离散元方法（DEM)模拟滑坡体大变形、不连续问题，以光滑粒子流方法（SPH)模拟水体的流体动力学特性，充分利用不同数值计算方法的特点，构建SPH-DEM流固耦合算法，并应用于2 0
2、 18 年10 月10 日西藏白格滑坡堵江的灾害模拟中。结果显示计算模型能较为准确地捕捉滑坡碎屑体入江激起水砂射流的规模以及范围，再现主堆积区向次堆积区产生二次滑移的动力过程，完整刻画堆积体的厚度、几何形态和断面特征。仿真结果有助于进一步深化对高速高位滑坡产生堵江过程和机理的认识，为后续进行抢险施工组织设计，制定防灾减灾策略提供参考。关键词：白格滑坡；堵江；流固耦合；离散元法（DEM)；光滑粒子流体动力学（SPH)中图分类号：P642.22DOI:10.14079/42-1745/tv.2023.02.0011 引 言我国西南山区受到青藏高原隆升影响，活动断裂带发育密度大，新构造活动强烈，深切
3、沟谷发育，地势高差大。活动断裂带附近斜坡岩体结构复杂，斜坡完整性差，断裂带岩体破碎，附近斜坡松散层厚度较大，极易发生高位滑坡链式灾害。高位滑坡具有高势能，下滑过程中伴随滑坡的解体，产生高速的碎屑流，冲击坡下江面后激起水砂射流，并堵塞河谷形成堰塞湖，产生二次灾害风险。2 0 18 年10 月10日（以下简称“10.10 滑坡”)和11月3日，在西藏自治区江达县波罗乡白格村与四川省白玉县绒盖乡则巴村交界处金沙江西藏岸先后两次发生大规模高位滑坡，堵塞金沙江，形成堰塞湖。滑坡堰塞湖淹没了上游的村庄和各种生产生活设施，同时溃坝洪水对滑坡下游的村庄、农田和公路、桥梁等基础设施的冲毁也十分严重。因此高位、高
4、速远程滑坡的启动过程、失稳机理以及滑坡碎屑流灾害链成灾模式研究成为当前国内外地质灾害防治领域重点研究内容。国内很多学者对白格滑坡产生的地质原因 1-2、诱发机制 3.4、动力过程 4-6、发展趋势 7.8、次生灾害 9及抢险救灾 10 进行了广泛的研究，取得了丰硕的成果，但仍然有一些未解的问题。例如在白格滑坡中，滑坡体撞击金沙江河水后产生了强烈收稿日期：2 0 2 3-0 2-17基金项目：湖北省教育厅科学技术研究计划项目“尾矿库溃决精细化模拟与环境安全技术研究”（项目编号：B2022598）作者简介：唐岳灏（197 7-），男，湖南湘潭人，副教授，博士，研究方向：水动力学、水资源优化运行、离
5、散单元模型开发和应用。文献标识码：B文章编号：16 7 3-0 49 6(2 0 2 3)0 2-0 0 0 1-0 8的水砂射流，在四川岸形成了一大片水沙射流区。由于射流的速度大，范围广，使得对白格滑坡产生的射流进行精细化分析成为一个巨大的挑战。这里面主要困难在于滑坡与水体间复杂的“流固”耦合作用以及射流传播过程中表现出高度分散的自由表面流。一些研究者为了计算简化，将滑坡体和水分别视为两种不同的流体，从而将其转化为求解多相流问题 5.8。但大多数滑坡的运动特征与流体有很大的差别,因此将滑坡简化为流体虽然简化了流固耦合计算分析的难度，但也势必会给结果带来一些无法解释的现象。其次，一些研究者受到
6、计算能力的制约，将三维的滑坡过程模拟简化为二维计算 5-7。也有研究者采用国外成熟三维软件进行仿真 8 ，但受限于软件二次开发的可扩展性，结果的展示显得较为单一。就目前所报道的数值模拟结果来看，大部分都没有直观反映堰塞体堆积高度、厚度、纵横剖面形态等基本几何参数，也无法准确判断口位置，更无法为堰塞坝应急抢险中导流槽的快速优化设计和施工提供参考依据。根据拉格朗日体系的力学理论，不连续介质力学计算方法，如DEM、D D A 等，可以很好地模拟大变形的运动过程，而流体力学方法如SPH、M P M等，则可以较好地模拟具有自由表面的流体运动过程。例如,徐文杰探讨了拉格朗日算法在滑坡涌浪研究应用中的可行性
7、 11；Zhao等人使用DEM与OpenFOAM软件二次开发的方法进行流一固耦合研究 12，Vacondio等人则使用三维SPH-DEM模拟了意大利瓦依昂大坝的滑坡涌浪事故 13。本文基于DEM和 SPH理论,充分利用不同的计算方法一1一唐岳灏，姜清辉SPH-DEM流固耦合方法模拟金沙江白格滑坡应用研究各自的优势，开发了基于并行加速的SPH-DEM耦合分析程序，更好地描述滑坡（固体）与水体（流体）间的相互作用。同时，由于本程序具有完全自主知识产权，可以无缝开发更为完善的后处理程序，实现滑坡灾害动力学过程模拟，这对于现阶段开发基于数字李生的流域模拟器，全景真实的展现、评估“滑坡堰塞湖溃决”灾害链
8、的全过程具有重要的理论和现实意义。2SPH-DEM 耦合分析方法2.1光滑粒子流体动力学(SPH)SPH是一种近似方法，其本质是将场函数及其导数的积分表示为其支持域中所有相邻粒子相应数值的求和。在SPH方法中，函数f()可以用积分表达式表示为：f(x)=-f(x)(x-x)dxf(x)W(x-x,h)dx2f(x,)W(x-x,h)2 j(1)=其中是Dirac delta函数，满足（x=0）=1和（x 0）=0。W 是平滑核函数。表示支持域,h是平滑长度。对于一个三维问题，点x代表一个粒子，2 表示在影响域中粒子所代表一个微元体积。由于体积可以表示为,=m;/pj，因此公式（1)也可被写成：
9、f(X)=25(x)W,本文中使用的核函数W是由以下公式给出的五次函数：W(r)=0(1-)(2g+1)00=i可以证明，公式（11)是一种精确的投影格式，即向量场在投影中的散度为零 14。另外可以看出，计算公式（11)需要其包含i，j,k 影响域内所有的粒子，这对提高数值精度和稳定性是有利的，因为核函数的二阶导数对粒子无序性非常敏感。最后，将公式（9），（10)和（11)带到流体控制方程（12），（13）中，即可求解流体粒子的密度场，速度场和压力场等物理参数。1d+V u=0dtd=-1vp+vVu+gdt其中，u为流体粒子速度矢量，为流体颗粒的一2 一其中，在三维问题中是2 1/（16 元
10、h），r 是两个粒子之间的距离，q=r/h。根据公式（3），粒子i的密度可以通过其影响区域内所有粒子的加权质量之和来估计。0;=20,W对公式（2)两边取导数：V.f(x)=2/(x)W.同样地，使用公式(2)或（4)来计算密度的梯度，可以得到：Vp:=2m,VW,SPH中的散度算子可以通过将放在公式(2)和(5)中来构造：V f(x)=(pf(x)-f(x)Vp(7)梯度算子可以用同样的方法得到：(8)将公式(4)和(6)中的和替换到公式(7)和(8)中得到:V.f(x)=-2m,f(x)-(x)VW,0;j=1(9)Vf(xi)=piZm;(fx)+f)vW,(10)(2)公式(9)、（1
11、0)是优于公式（6)和（7)的，因为场函数f()以成对粒子的形式出现，因此能保证线动量守恒性质。根据拉普拉斯算子表达式=（V f），将（10)代人公式(9)中，得到如下形式。2（3)密度，p为压力，g为外力，是运动粘滞系数。方程（12)描述质量守恒，即任意粒子所代表流体微元内流速的散度和其密度的变化满足平衡关系。动量方程（13)表明，粒子的运动由压力梯度p，黏性力vu以及外力g造成。本文采用投影法对（12)和（13)进行离散。投影法是一种算子分割算法，其理论基础是任何矢量场都可以分解为无散分量和无旋分量的梯度。由于压力是一个非旋标量场，即p=0,因此速度场可以被分解为流动粘度和外力(12)引起
12、的中间速度和压力引起的压力速度。最终压力(13)场可以通过求解半隐格式的泊松方程获得。具体过程可参考相关文献 15 mi=2m,W=1(4)(5)(6)(11)2023年6 月2.2描述颗粒摩擦特性的离散单元模型(DEM)DEM方法用牛顿第二定律的微分形式计算固体粒子的位置、线速度、角速度。假设固体粒子都是质量均匀的标准球体，粒子i沿v方向运动并围绕w,轴旋转，粒子j沿v,方向并围绕w轴旋转（图1)。dvimig+ZF,+Fadtmidpi=VidtFlgXr;n,+FlaXr;nado:=idt其中,F;是粒子i对粒子i的力,Fa是边界对粒子i的力，g是重力加速度常数，m;是粒子的质量，V;
13、是粒子的速度，P;是粒子的位置矢量，;是粒子i的角速度,Fli是粒子i和粒子之间的切向力。r;是粒子的半径，I是粒子i的惯性矩,Fila是边界给粒子i的切向力。n.是与i粒子作用边界的外法线单位向量方向，n是从i粒子的位置指向粒子位置的单位向量。Vi图1DEM计算模型根据摩擦理论，粒子之间的动摩擦是由粒子间的相对速度和压力决定的，即：Fl;=hl,n,-CiVliFali=-min(Fl,l,CrVl)vl(18)其中C是切向粘壶的阻尼常数,CI是法向粘壶的阻尼常数,Fli是粒子i和粒子i之间的法向力，Fai是粒子i和i之间的切向动摩擦力，h是法向胡克系数，l;是球体i和j之间的重叠距离（l,
14、0)，V i 是粒子i和i之间的相对速度，Vl;是V;的速度分量，方向垂直于粒子中心线。类似地，V工；是平行于粒子中心线的速度分量。是动摩擦系数。F,=Flj+FliVli=(Vin)njVi=Vij-VliV,=(V;+o;Xr;nj)=(V,+o,rj nj;)(22)长江工程职业技术学院学报过大。Fug=min(|Fll，|F s )Frlj（2 4)Fslj=-h,Sj，/v li|U m a x(25)式中,F;粒子i和i之间的静摩擦力,h,是切(14)向虎克系数，s;是垂直于粒子i和i法向力的滑移距离，满足以下近似微分关系。Umax是产生静摩擦(15)的最大速度阈值，超过此值即产生
15、动摩擦。dsi=-(P:P)Si (V-V)(16)dt当所有的位置、线速度、角速度和滑移合并为一个矢量y时，整个计算可以写成一个微分公式的形式：y=f(t,y),y(to)=yo方程(2 7)可以通过欧拉法或龙格库塔法求解，这里不再赞述。2.3SPH-DEM耦合模型流-固耦合过程中，流体的运动将对固体产生驱动力，而固体作为流体的边界对流体的运动产生阻碍作用。由于SPH和DEM都属于拉格朗日体系的方法，可以将粒子间相互作用力看作是粒子间的“碰撞”。因此，基于SPH和DEM的流-固耦合计ni算，实质上是将流体作用于固体颗粒表面上的作用力进行积分，从而得到流体对固体的作用力；而DEM固体颗粒表面元
16、素包含在与其相邻流体颗粒P,的SPH方程内插中，从而计算固体对流体的作用。Pi在SPH的计算模型中，固体粒子被视为流体粒子也0参与到N-S方程的计算中。从动量方程（13)可知，固体粒子的运动主要受到周围流体粒子的压力和粘性力的作用。将压力带人方程（10)可以得到固体粒子i受到周围流体粒子f的作用力为：(17)Fi,f=m;mf(根据作用力和反作用力原理，流体粒子f受到固体粒子i的作用力为F,=-F,f，即流固作用是相互耦合和动量守恒的。粘性力II 可采用层流的计算方法，即固体粒子i受到周围流体粒子对其产生的粘性力为：IIi.f=m;mf固体运动的微分方程（14)则改写为：(19)m:dt4=m.g+2F.+F.+Z,+21,(30)(20)(21)第40 卷第2 期摩擦力被摩擦系数和法向力所限制，以防止其Fllj=Fsllj+Falli(23)(26)I p;-p;12(27)（+）VWi.f(p%由于SPH和DEM在时间步积分上都是显格式，对步长有要求，因此必须将两者进行协调。本文一3 一(28)(29)唐岳灏，姜清辉SPH-DEM流固耦合方法模拟金沙江白格滑坡应用研究采用的半隐格式