陈闻潇 石 崇 单治钢 张一平 李汪洋

(①河海大学， 岩石力学与堤坝工程教育部重点实验室， 南京 210098， 中国) (②河海大学， 岩土工程科学研究所， 南京 210098， 中国) (③中国电建集团， 华东勘测设计研究院有限公司， 杭州 310014， 中国)

0 引 言

水下滑坡作为在近陆范围常见的一种地质灾害，近年来已造成严重的人员伤亡或财产损失(Alves et al.，2015; Terry et al.，2017;Sassa et al.，2019)，仅在1616～1886年期间，世界上由于地震、火山喷发等引发的水下滑坡失稳事件达333次(Milne，1897)。此外水下滑坡对海洋设施如海洋油气生产平台、水下管道、水下电缆等也具有较大威胁(朱超祁等， 2016)，因此研究水下滑坡的形成过程与机理具有重要的工程意义。

由于水下滑坡的特殊环境，数值模拟是研究水下滑坡的重要手段。杨林青(2012)基于ABAQUS有限元强度折减法分析了水下斜坡稳定性，霍沿东等(2019)基于极限分析上限方法分析了水下斜坡在地震下的稳定性。但这些方法采用的连续力学分析方法，无法模拟分析滑坡造成的大变形破坏、滑坡运动过程、堆积形态与灾害影响范围等(陈闻潇等， 2020)。

颗粒离散元(PFC)方法具有模拟大变形破坏的优势而被广泛运用到滑坡过程分析中(陈晓等， 2020)，但目前基于PFC方法的滑坡研究大多为陆上滑坡(曹文等， 2017； 周礼等， 2019)，对水下滑坡研究不足。且在PFC中无法直接考虑水的作用，常见的达西渗流法(毕庆杰， 2018)、管道域法(Zhang et al.，2019;Zhang et al.，2020)分别适用于多孔介质渗流和水力劈裂，但无法模拟水下大尺度流体环境，故需要考虑引入流体力学计算方法，才能实现流固耦合计算。

OpenFOAM是一款计算流体力学的开源软件，相比商业软件，其具有良好的扩展性、丰富的模型库和求解高效的优势(Shojaeefard et al.，2017)。已经用来研究波浪、风场等方面的流体作用(Jacobsen et al.，2015; Liu et al.，2017)，但其在土木和地质方面鲜有研究，有待进一步开发。

本文通过PFC内置TCP套接字模块(Itasca， 2014)，实现基于OpenFOAM与PFC耦合计算方法，并针对PFC流体分析提出适用于该耦合方法的相似性方法。然后通过典型案例，验证该耦合方法进行水下滑坡分析的可行性，进一步研究水下滑坡的动力学特性及堆积形态，并与单向耦合水下滑坡和陆上滑坡结果进行对比分析。

1 耦合分析实现方法

在考虑滑坡与水动力学耦合作用时，流体动力学需要考虑连续性方程和动量守恒方程，此处采用考虑颗粒作用修正后的连续性方程和Navier-Stokes方程(Jackson， 2000)：

(1)

(2)

式中：V*为流体速度矢量；n为流体单元内的孔隙率;ρf为流体密度；p为水压力；t为时间；τ为黏性力张量；g为重力加速度；fint为单个流体网格内包含颗粒所受的拖曳合力。上述方程的计算使用OpenFOAM中的demIcoFoam求解器，demIcoFoam通过修改icoFoam求解器来模拟固体颗粒的作用，同时添加颗粒与流体之间的相互作用力和基于粗网格法(Wang et al.，2013)得到的单元孔隙率。

单个体积为Vi的网格内的拖拽合力fint通过PFC中的CFD模块计算，计算公式为：

(3)

式中：fdrag为单个颗粒所受流体对其的拖曳力，其通过下式计算：

(4)

式中：r为颗粒半径；v为流体速度；u为颗粒速度。孔隙率n通过基于粗网格法的多面体法计算，n-χ项是考虑局部孔隙度的经验系数。这个修正项使拖曳力同时适用于高孔隙和低孔隙度系统，流体雷诺数也可在大范围内取值。拖拽力系数Cd被定义为：

(5)

式中：当流体的动力黏滞系数为μf时，颗粒的雷诺数Rep通过下式计算：

(6)

上述耦合方法，基于PFC中提供的TCP套接字(Itasca， 2014)连接OpenFOAM与PFC，实现耦合计算中信息传递与更新。每一个流程步中，OpenFOAM将流体速度v、压强p和压力梯度p传递给PFC，而后在PFC中进行式(4)～式(6)的计算以得出每个颗粒所受拖拽力fdrag，继而通过力-位移运动方程得出土体颗粒的运动状态，并通过式(3)～式(6)更新单元格内拖拽合力fint，将其与单元格孔隙率n一起传递给OpenFOAM实现式(1)～式(2)的计算，这样便完成了一个耦合循环步的计算。耦合流程如图 1 所示：

图 1 耦合交互方法流程图Fig. 1 Flowchart of coupling interaction method

用颗粒水下自由沉降来验证该方法的可行性(Jing et al.，2019)，将半径为5 mm的土体颗粒置于流体环境中，计算颗粒在重力下自由运动0.5 s后的运动速度。图 2为运动速度随时间的变化曲线，其理论解由式(7)计算，可得颗粒速度的理论解与模拟值一致，可验证该流固耦合方法正确性。

(7)

图 2 不同时刻的颗粒下落速度Fig. 2 Falling velocities of particles at different times

在流固耦合中由于模型尺寸或颗粒尺寸相比真实值变化，雷诺数、流速、压力差、启动力等均会变化，需要考虑相似性理论(孙其诚等， 2008； 蒋明镜等， 2014)。现有的相似性方法主要有两种：粒径缩放-流体参数调整法和粒径恒定-流体参数调整法(倪小东等， 2009)。在离散元计算中，通常以大的离散颗粒来替代真实试验中的土体颗粒(Liu et al.，2015)，故相当于增大了颗粒尺寸同时模型尺寸不变，因此现有的相似性方法无法很好地模拟这种情况。通过调节流体参数的方法均无法较好地模拟，其会造成流场结果失真，故需要提出一种基于离散元PFC的相似性方法。

当颗粒尺寸为真实值的N倍时，本文采取程序二次开发的方式，通过修改内部函数使得模型满足相似性原理。首先通过二次开发，将雷诺数Rep赋以系数1/N，用以保证颗粒半径扩大N倍的同时雷诺数(式6)不变，即满足雷诺准则。同时将拖曳力fdrag赋以系数N，用以保证拖拽力与重力同步增长，满足启动力相似(倪小东等， 2009)，同时整体满足弗洛德(Froude)准则v2/gl=c1和欧拉准则Δp/ρv2=c2。其中拖曳力fdrag的扩大通过PFC中FISH函数二次开发，在CFD_AFTER_UPDATE注册点修改函数值实现。此外保持其他流体参数不变。

表 1显示了不同半径的颗粒，在通过该相似性方法模拟后，颗粒在重力下自由运动0.5 s后的平均速度和平均位移。为保证放大倍数变量的单一性，防止孔隙率变化对结果的干扰，故取颗粒总体积一致。由表 1 可见颗粒平均速度和平均位移基本一致，误差较小，可认为该相似性方法具有一定的可行性与准确性，并且通过放大颗粒半径，可以缩小模拟所需颗粒个数，极大地缩短计算所需时间。

表 1 相似性参数及验证结果Table 1 Similarity parameters and results

2 水下滑坡的耦合计算分析

基于典型水下滑坡案例，分析OpenFOAM与PFC耦合水下滑坡模拟方法的可行性，并进行水下滑坡的动力学特性及堆积形态研究。

为兼顾问题的三维条件和耦合方法的可拓展性，耦合方法建立于三维条件下，可以通过假三维模型来进行平面问题的模拟，如图 3 所示建立连续边坡模型，假三维厚度为1 m。基于随机种子法生成水下滑坡的离散元模型，边坡模型长120 m，高25 m，共生成颗粒6747个，颗粒粒径大小在0.25～0.5 m之间，呈线性分布。边坡右侧长80 m，使得滑坡体能够充分地运动堆积，同时避免在流体计算时边界效应的影响。对水下土体的模拟采用接触黏结模型，在颗粒破碎之后，接触模型由接触黏结模型，退化成线性接触模型，能够较好地模拟土体的力学响应和破坏特性(Shimizu， 2004)。本文模型参数如表 2 所示。

图 3 初始计算模型及网格划分Fig. 3 Initial calculation model and mesh generation

表 2 计算模型细观力学参数Table 2 The meso mechanical parameters of calculation model

采用固定粗糙网格法建立流体网格模型，网格尺寸大小为1 m 的正方形，单元数量为7200。同时假定研究区域内水下土体平均粒径为0.001～0.002 m之间，结合上文提出的放大系数为250，在OpenFOAM的demIcoFoam求解器中设置流体的密度ρf=1000 kg · m-3，动力黏滞系数μf=0.001 Pa · s。此外在耦合模块中，在CFD_AFTER_UPDATE注册点对雷诺数Rep和拖曳力fdrag赋以放大系数。

水下滑坡的耦合计算结果如图 4 所示，后文描述时间为滑坡实际时间。图 4a，图4c分别为模型计算至29.9 s和39.8 s时的土体位移图，可见滑坡呈现出较为明显的圆弧滑动面，随着滑坡的发展，滑坡土体沿滑动面运动，滑动区域下部的土体不断推进，边坡后缘土体逐渐破坏，土体最大位移分别为26.0 m和33.9 m，其堆积至坡脚的范围分别为14.3 m和20.1 m。值得注意的是，滑坡体前端厚度较大，且前端呈现椭圆状，同时在滑坡体的前端底部出现了滑水现象(修宗祥等， 2016)，其主要由于滑坡体快速运动而在底部形成了类似水楔的水层，这同时也和景路等(2019)研究中得到的现象相符。由图 4b 可见滑坡土体运动带动流体运动并形成了漩涡状，流场速度分布基本与土体所在区域重合，远离滑坡区域的流体速度较小。随着滑坡的发展，由图 4d可见边坡后缘部分由于土体滑落，原本由土体填充的区域由流体涌入，并在该区域内形成了杂乱的紊流。

图 4 耦合模型计算结果Fig. 4 The calculation results of the coupling modela. 29.9 s位移图(单位：m)； b. 29.9 s流体速度矢量图(单位：m · s-1)； c. 39.8 s位移图(单位：m)； d. 39.8 s流体速度矢量图(单位：m · s-1)

滑坡结束后堆积状态如图 5a 所示，滑坡体最终在坡脚堆积并停止进一步运动，堆积距离为27.6 m，堆积体前端变薄，不再呈现椭圆面。由图 5b 可见滑坡体孔隙率要明显高于未滑动土体，其孔隙率值处于0.4～0.5之间。这说明随着滑坡的发展，土体由紧密接触的整体，破碎成了松散的团体。

图 5 滑坡结束后堆积状态Fig. 5 The accumulation state after the landslidea. 滑坡体位移图(单位：m)； b. 网格孔隙率图

滑动土体的速度位移曲线如图 6 所示，可见在滑坡发展初期速度迅速增大，在15 s左右土体的平均速度达到峰值0.47 m · s-1，峰值的迅速出现是因为流体对颗粒的拖拽力fdrag与颗粒速度方向相反，大小与流固速度差相关，故会限制土体速度进一步发展，同时由于流固之间的相互作用，滑坡土体的能量一部分消耗于带动流体的运动。在20～80 s内土体速度不断减小，在50 s时土体平均速度下降到0.091 m · s-1左右，之后土体进入缓慢的堆积平衡阶段，此时流体对颗粒的拖拽力fdrag会减缓土体颗粒速度下降的趋势。如果该研究对象处于有海洋洋流的存在情况下，将会使得堆积体在流体的作用下运移出较远的距离，这将会导致滑坡影响范围更进一步地扩大。

图 6 滑坡体的平均速度位移曲线Fig. 6 Average speed and displacement curve of landslide body

3 耦合对比分析

为对水下滑坡的耦合计算方法进一步探究，并讨论水下滑坡与陆上滑坡的异同点，本节研究基于OpenFOAM-PFC的双向耦合水下滑坡、单向耦合水下滑坡与陆上滑坡(无流固耦合)3种不同的耦合方式。单向耦合是当前滑坡分析常用方法之一(丁欢欢等， 2013)，其在计算过程中信息只由流场传递至颗粒速度场，而颗粒运动不会对流体产生影响，即后续计算中流体速度v恒定不变(式4)。两个模型均采用上文的模型尺寸与土体参数，但由于陆上滑坡初始条件下不受水力作用，故需要通过强度折减法调整边坡的接触强度而导致其产生滑坡。

可见陆上滑坡( 图 7a) 和流固单向耦合滑坡( 图 7b) 与基于双向耦合的水下滑坡( 图 4) 滑坡状态并不一致，两者在滑坡过程中滑坡体前端厚度较小，均没有椭圆面出现，无滑水现象产生，这是因为单向耦合模拟中颗粒对流场作用无法体现，流固相互作用机理无法实现。陆上滑坡由于没有流体作用，滑坡时间更短，主要位移发生在0～10 s内，滑坡过程中的颗粒碰撞相对剧烈，大部分接触黏结已经破坏，滑坡体相对更加松散。而单向耦合中由于流体速度恒定不变，较大的流固流速差使得土体运动缓慢，滑坡过程持续达300 s，缓慢的运动过程使碰撞和剪切相对较弱，土体相对仍较密实。

图 7 滑坡位移曲线与运动状态图Fig. 7 Diagram of landslide displacement curve and motion statea. 陆上滑坡； b. 流固单向耦合

从图 7a 可以看出，不同耦合方式模拟的运动规律并不一致。双向耦合和陆上滑坡的平均速度峰值分别为0.47 m · s-1和1.95 m · s-1，出现的时间分别在14.3 s和3.0 s。而单向耦合中并不存在明显的速度峰值阶段，其滑坡速度在10～200 s内维持在0.15～0.2 m · s-1范围内，并随着时间逐渐降低。这是由于其假定流体速度v恒定不变，随着颗粒速度u增大，拖拽力fdrag不断增大进而当合外力趋于0时，颗粒速度不再明显增加，且只因颗粒接触力变化而小幅度变化，以上说明了单向耦合方法具有一定的局限性。受流体作用影响，水下滑坡的滑坡发展与运移要明显滞后于陆上滑坡。从图 8b 可以看出，单向耦合模拟的位移要小于双向耦合，这是由于单向耦合中更大的流固速度差，导致更多的能量被消耗在了流固相互作用中。而陆上滑坡由于运动剧烈，接触破坏更多，颗粒碰撞更加剧烈，导致能量损耗迅速，故平均位移要小于双向耦合结果。总体而言，陆上滑坡的位移要小于水下滑坡，且真实情况下由于水下存在水流作用，其可能会带动滑坡土体进一步运动，产生更大的破坏范围。

图 8 滑坡体速度和位移曲线Fig. 8 Velocity and displacement curves of the landslidea. 平均速度曲线； b. 平均位移曲线

4 结 论

本文基于PFC和OpenFOAM建立了流固耦合计算方法，提出了适用于该耦合方法的相似性赋值方法，利用案例分析了水下滑坡的动力学特性及堆积形态，并与基于单向耦合方法的水下滑坡和陆上滑坡结果进行了对比，得到主要结论如下：

(1)颗粒水下自由沉降模拟得到的计算值与理论计算值一致，表明了该耦合方法的正确性，同时对不同粒径的颗粒采用相似性方法考虑，验证了该方法的可行性。

(2)水下滑坡的运动形态表现为滑坡体前端厚度较大，呈椭圆面并有滑水现象的产生，由于流固相互作用，拖拽力会阻碍颗粒运动，滑坡土体运动速度较快达到峰值并缓慢下降，堆积土体孔隙率明显大于未滑动土体。结论表明该耦合方法能够较好地进行水下滑坡过程模拟。

(3)水下滑坡与陆上滑坡在运动特征、堆积形态和滑坡历程上具有较大的差异。单向耦合方法存在一定局限性，PFC-OpenFOAM双向耦合能更好模拟水下滑坡运动特征。