基于PFC-CFD的土体管涌数值模拟研究

陈闻潇 石 崇 单治钢 陈 亮 倪小东 李汪洋

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

管涌是导致河道堤防、坝基等失稳破坏的重要因素[1-2].在渗透水流的作用下,填充在土体骨架中的细颗粒逐渐流失,管涌通道的直径将不断扩大,导致堤坝发生不均匀沉陷,最终导致崩塌[3-5].

近年来,越来越多的学者对管涌现象展开了研究.理论计算方面:施倩等[6]建立了管涌破坏的数学模型,介玉新等[7]研究了管涌破环的时间过程模拟.室内试验方面:曹洪等[8]研究了管涌通道边壁松散层对管涌发展的影响,Mohamed Elkholy等[9]研究了土壤组成对土堤管涌的影响.数值模拟方面:常利营等[10]研究了砂土渗流规律的适用性,Y.Guo等[11]研究了不同颗粒形状对管涌的影响.离散元方法如颗粒流(PFC)[12],因为能够较好地模拟大变形破坏,而逐渐被较多的学者采用.周建等[13]通过PFC模拟了砂土管涌的基料-滤层系统,从细观角度揭示管涌发展过程中颗粒的运动特性和滤层防治机理.倪小东等[14]基于PFC对室内砂槽试验进行了模拟.张刚等[15]通过PFC从细观角度研究了管涌的形成和发展机制.但学者们的模拟大多基于室内试验或小尺度模型,针对现场管涌案例建立模型的离散元分析还有待进一步研究.

本文基于PFC模拟管涌现象的发展,渗流流体的CFD计算通过python中的fipy偏微分方程求解器求解达西公式来实现,通过将PFC-CFD结合来实现土体-流体的相互作用计算.参考南京长江大堤管涌案例,建立管涌计算模型,分析随着管涌发展,渗透流速、水压力和土体孔隙率的变化规律和动力学特性,并研究堤坝位移的变化趋势.

1 基于PFC-CFD的多孔介质渗流计算方法

利用PFC-CFD模块计算三维条件下的多孔介质流动,多孔介质中的低雷诺数流动通常可以通过达西定律描述:

式中:v为流体的速度矢量;K为渗透矩阵;μ为流体黏度;ε为孔隙率矩阵;p为流体压力.通常假定流体的压缩性很小,可以忽略不计,即认为流体不可压缩:

稳态不可压缩渗流方程可以通过对方程(1)两边同时取散度导出:

将公式(3)代入(2),得泊松方程:

式(4)边界条件为入口处有

其中:vin为入口速度;出口处有p=0.这个方程通过隐式求解可以很容易得出流体的压力场.求解方案基于稳态流,即流入量与流出量相等,一旦压力已知,流体速度可以由方程(1)直接获得.

流体方程(1)和(2)是在粗流体网格单元[16]上求解.通过计算PFC3D颗粒与流体单元之间的重叠量确定孔隙率.考虑流动受颗粒运动的影响,渗透系数由PFC3D模型的孔隙率计算.

当re为PFC颗粒平均半径时,第i个粗网格的渗透系数Ki可以通过其孔隙率εi和Kozeny-Carman关系[17]来计算:

2 南京长江大堤管涌案例

2016年7月3 日晚,南京长江扬子江段江堤背水坡坡脚出现渗水现象,本文通过建立PFC-CFD计算模型,实现管涌发展的动态分析.如图1所示,出险段长江大堤堤顶高程12.6 m,堤顶宽8.0 m,背水坡堤脚高程约8.5 m,长江下关最高潮位为9.67 m.如果渗水不断发展,渗透水流逐步掏走通道中的土体,通道不断扩大,大堤就会坍塌.因此南京市防办及相关部门及时采取应急处置,渗水情况得到明显控制,大堤安全总体可控.本文基于该案例进行PFC-CFD方法的管涌过程模拟分析.

图1 南京长江大堤管涌案例

3 管涌计算模拟

如图2所示建立边坡模型,在此基础上考虑管涌现象的发展及其对堤坝的稳定性影响.

图2 管涌计算模型

由于耦合方法是基于三维模型,为兼顾问题的三维性,采用假三维模型来进行平面问题的模拟,假三维厚度为3 m.模型尺寸如图2(a)所示,整体模型长50.9 m,高12.4 m,土体的孔隙率在0.36左右.由于计算能力受限,将颗粒粒径扩大10倍,管涌影响区域内细颗粒和粗颗粒半径分别为0.04～0.08 m和0.12～0.24 m.管涌影响区域之外均为0.12～0.24 m半径的颗粒,为提高计算效率,共生成颗粒63644个.考虑颗粒的尺寸效应[18],需要对流体的黏滞系数扩大10倍以满足相似性原则.

如图2(b)所示,基于python调用gmsh软件建立非均匀六面体网格,共生成3612个网格,并将网格信息导入到python中的fipy偏微分方程求解器中,通过fipy实现公式(1)～(6)的计算.通过PFC中的CFD模块为颗粒施加流体-颗粒相互作用力,在每一个时间步更新孔隙率和渗透系数,并将数据传递回python中的fipy求解器中,随之进行下一个时间步的计算.模型采用接触粘结模型,该模型能够较好地模拟土体的力学响应和破坏特性[19].堤防土体细观力学参数见表1.

表1 计算模型细观力学参数

管涌计算结果如图3所示,如果从细颗粒逐渐启动流失到管涌逐渐形成开始模拟,则需要模拟的时间过长难以接受.故从长江潮水位9.67 m开始模拟,此时已经具有一定的初始渗透流速.从图3(a)、3(b)可见,土体中各流体单元孔隙率的不同导致渗透系数不同,引起了各测点流场中压力及流速的不均匀.随着管涌的发展,各个测点的渗透性也在不断地增强,渗透流速不断增大.

图3 管涌模型计算结果

从图3(b)可见,管涌由背水侧的细小颗粒流失而不断发展,渗透流速不断变大,进而导致整个渗流通道内的土体开始流失.近水侧监测点1在15万～20万时间步时,孔隙率急剧上升,表明这部分的土体已经大部分流失而接近破坏.背水侧土体虽有流失,但同时水流也将近水侧流失颗粒携带而来,故孔隙率只在较小幅度内变化.从图3(c)可见前10万时步时,孔隙率的变化幅度在0.1左右,10万～40万时步,孔隙率的变化幅度达到了0.2～0.4.

从图3(d)可见,在渗流初期,水压力随着渗透路径呈现近似地均匀分布,随着管涌发展,靠近渗流入口处的土体逐渐流失,水压力上升明显.从图3(e)可见,初始各个位置的水压力梯度是近似均匀的,在2万～20万步内,沿渗透路径0～10 m处的水压力梯度逐步下降,是因为该处的土体不断流失(图3(c)),而10～20 m处压力梯度逐步提高.计算至40万步时,0～15 m处的压力梯度陡然下降,这是由于此处已达到土体的临界启动速度而发生土体流失(图3(c)).此时的压力梯度集中在了20～30 m处,随着管涌的进一步发展,这一部分土体也将会达到临界启动速度并流失.

模型对土体的流失量进行了测量,图3(f)为模拟过程中流失量随时间的变化曲线,20万时间步时流失量已经达到了2.5 m3,在40万时间步流失量已经迅速增加至9 m3,巨大的流失量表明土体中可能已经产生了空洞,这对于整个堤坝的安全是十分不利的.随着管涌的进一步发展,流失量将会继续增大,这可能会导致堤坝发生破坏.

图4 模型位移云图(单位:m)

由图4位移云图可以看出,渗流初始阶段堤坝位移主要发生在堤坝内部渗流通道附近,在2万、10万、20万、40万时步下的最大位移分别为0.042、0.068、0.093 6和0.148 m.可见位移是逐渐扩大的,主要发生在堤坝内部和堤内坡脚处,这对于堤坝土体是十分不利的,如不及时治理将会导致更进一步的破坏.

4 结 论

本文采用PFC-CFD数值模拟方法,基于南京长江大堤管涌案例,建立了土体管涌的数值模拟模型,研究了管涌渗透流体变化规律,分析了土体孔隙度及位移随管涌发展的变化趋势.主要得出以下结论:

1)采用PFC-CFD的模拟方法可以较好地用于土体管涌现象的模拟,可以较好地呈现渗透水流与土体的相互作用.

2)随着管涌发展,靠近渗流入口处的土体因为管涌渗透作用而率先流失,沿渗流路径上渗透流速随着管涌发展而不断加大,水压力和压力梯度出现分布变化,土体流失量迅速增加.

3)堤坝位移主要发生于堤坝内部渗流通道附近,随着土体不断流失,堤内坡脚处的位移不断增大,将导致堤坝破坏.