夏 霄，甘鹏路

（中国电建集团华东勘察设计研究院有限公司，浙江 杭州 311122）

1 概述

管涌作为一种土体内部侵蚀现象是基坑工程中的常见灾害，它引起的基坑挡墙倒塌事故时有发生[1]。例如，2015 年杭州某滨江基坑开挖中发生了管涌诱发的挡墙倒塌事故，承受较高水力梯度的砂质土体从挡墙裂缝处涌出，造成坑外土体大幅沉降和挡墙的严重变形[2]。类似基坑挡墙薄弱处（如桩或地连墙连接处）破裂引发的管涌灾害屡有发生[3-4]。

正是由于管涌细观机制的复杂性，相应的数值分析研究也逐步展开。一类是利用有限元将土体视为连续性介质并与渗透水流进行耦合的唯象分析[5-7]，它很难捕捉管涌中土颗粒流失，重新排列、沉积等重要特征。因此，另一类基于流体力学（CFD）和散体介质（DEM）的多相耦合仿真分析方法（CFD-DEM）在管涌研究中越发常见[8]。这种方法将土体视为不同相的组合，如土骨架以及土骨架孔隙中的水流与细颗粒，它可考虑内部侵蚀或管涌过程中细颗粒从土骨架中的流失。CFD-DEM耦合仿真分析方法可捕捉基坑管涌过程中土体的主要力学特征，也可用于解决实际基坑管涌问题[9]。

本文利用CFD-DEM 耦合仿真分析方法研究细粒含量（Fc）与水力梯度（i）对间断级配土基坑中管涌灾害的影响。首先介绍了CFD-DEM 流固耦合仿真模拟方法，然后通过建立基坑管涌模型，考虑不同的Fc和i进行一系列参数分析，得到不同工况下基坑土体颗粒速度场和流场特征，验证了该方法在基坑管涌模拟中的优越性。

2 CFD-DEM流固耦合仿真分析方法

CFD-DEM 方法由三个部分构成，分别是离散元法（DEM），计算流体力学（CFD）及CFD-DEM耦合算法。DEM主要用于对土颗粒的仿真，CFD则是仿真土中的渗透水流。CFD-DEM耦合算法主要计算流体—颗粒间相互作用力，以及实现CFD 与DEM 间信息交互。开源 DEM 软件 LIGGGHTS 和 CFD 软件 Open⁃FOAM 是本研究中具体使用的多相耦合仿真分析工具。下面简要介绍这三个仿真分析方法的基本原理和控制方程。

2.1 DEM方法控制方程

DEM 中颗粒运动遵循牛顿第二定律,颗粒i的平动和转动由下式计算：

式中：mi、It——颗粒i的质量和转动惯量；

Ui、ωi——颗粒i的平动速度与转动角速度；

本研究中颗粒间接触模型采用Hertz 非线性接触定律[10]。对于切向上的滑动阻力，利用Mohr-Coulomb准则计算。

2.2 CFD方法控制方程

CFD 域中，流体运动主要由连续性方程和动量方程控制，如下所示：

式中：Uf、p——流体网格的平均速度和压强；

εf——流体网格中颗粒体的孔隙率；

fp——单位体积颗粒与流体间平均相互作用力，即

ΔV——流体网络体积；

k——颗粒数；

ρf、μf——流体密度和粘度；

g——重力加速度，本文中为50倍的常重力速度。

2.3 流体—颗粒相互作用力计算

在管涌问题中，流体—颗粒相互作用力Ff主要包括拖曳力、粘滞力和动压差力。本研究中也主要考虑上述三种流固相互作用力，忽略其它相互作用力对管涌特性的影响。Ff的计算如公式（3）所示：

流体对土颗粒产生的拖曳力采用Buijtenen 等（2011）[11]方法计算，如公式（4）所示，它适用于颗粒体密度较大的情况，在流体雷诺数变化较大时也成立：

式中：dp、Up——颗粒的直径与速度；

Vp——单个颗粒的体积；

εp——流体网格中颗粒所占体积百分比，εp+εf=1；

Rep——单个颗粒的雷诺数。

作用于颗粒上的动压差力和粘滞力计算见公式（5）和（6）[12]：

式中：τ——流体剪切力。

3 基坑管涌仿真分析工况与参数设置

3.1 仿真分析工况

本文主要通过不同的土体细粒含量Fc与水力梯度i工况来研究其对基坑管涌的影响。一般认为Fc直接决定土体组构，采用Fc=15%和35%的试样研究两种不同的典型土体组构对基坑管涌的影响。

为研究基坑管涌启动前后土体力学性质演变，在仿真过程中，水力梯度由i=0.8 增至1.8（i= Δp/ρgL，Δp是基坑内处水压差，L是模型中流体最短流动路径，本研究中是0.5m，g是模型中采用的重力加速度值）。当i=0.8时，管涌并未启动，而i增至1.8时，管涌开始发生并持续至仿真结束。表1总结了本研究考虑的四个管涌仿真工况。图1是本研究中采用的土体的级配曲线，粗粒组和细粒组均为单一直径颗粒，分别为2mm和0.5mm。考虑到CFD-DEM 耦合仿真计算效率，每一种工况的仿真时长均为1.0s。结果表明，在此较短的仿真时长内，基坑变形也可达到稳定状态。

图1 试样级配曲线

表1 仿真分析工况

3.2 模型尺寸与参数设置

图2 是基坑DEM 模型的尺寸示意图。由于本研究主要针对管涌过程中颗粒与流体力学特性演化，因此挡墙假定为钢性体，在仿真过程中不发生变形。为研究挡墙裂缝引起的基坑管涌灾害，在挡墙10cm高度处开有2.5cm 高度的孔洞。挡墙右侧和左侧分别是基坑的开挖侧与支护侧，高度分别为50cm 和25cm，长度分别均为25cm和50cm。整个模型厚度为5cm，将实际的基坑管涌简化成平面应变问题。结果表明这一尺寸设置可以减小边界效应的影响。仿真模型的CFD域与DEM域范围大小相同，保证所有颗粒均受渗透水流的作用。表2列出了模型中颗粒与流体的力学性质参数。

图2 基坑管涌仿真DEM模型

表2 CFD与DEM参数汇总

3.3 边界条件设置

为模拟实际中基坑挡墙两侧水压差，在CFD 域的开挖侧顶面施加一定水压，使得基坑内外水力梯度在整个仿真过程中均维持在0.8或1.8。CFD域其余面均为不透水面。对于颗粒与DEM边界的接触，假定它们之间的摩擦系数等于0，法向接触刚度是颗粒弹性模量的10倍。

3.4 仿真分析步骤

首先在基坑区域上方生成一定数量的颗粒并使其在重力作用下自由下落。在填至模型1/5高度时，将颗粒体竖向缩至一定的相对密实度，本研究中将土体相对密实度设定为Re=65%～67%。重复上述两个步骤直至颗粒填满基坑支护侧与开挖侧。在DEM模型准备完成后，立即在CFD域的支护与开挖侧顶面上施加固定的水压力差以开始基坑管涌的仿真模拟。所有颗粒及它们间的接触信息在整个仿真过程中每隔0.01s输出。

4 基坑管涌仿真分析主要结果

图3是细粒含量分别为Fc=15%和35%情况下基坑管涌启动时刻（t=0.2s，i=1.8）地层变形与颗粒速度场。在管涌启动时刻，挡墙裂缝处和基坑开挖侧坑底土体首先开始流动，具有较大的初始速度，而支护侧土体相对流动速度较小。图4 是在Fc=15%和35%时，基坑管涌仿真结束时刻地层变形与颗粒速度场。此时，支护侧土体出现明显沉降，开挖侧坑底土体明显隆起。在挡墙裂缝前方出现一条管涌通道，支护侧土颗粒经此流入开挖侧并逐渐趋于稳定。当土体中含较多细粒时，支护侧土体沉降较大。这主要由于裂缝宽度与颗粒粒径之比较大时，颗粒不易堵塞，更易通过裂缝，因而在相同水力条件下流动性更强，导致最终支护侧颗粒流失量增大。

图4 最终时刻(t=1.0s)基坑变形与颗粒速度场

图5 是管涌过程基坑流场的CFD-DEM 仿真结果。基坑支护侧水体在水头差作用下通过挡墙裂缝流入开挖侧，裂缝周围水流速度达到最大。

图5 基坑管涌过程流场仿真结果

此外，通过比较不同工况下支护侧土体流失量，可以发现，对于土体细粒含量较低（Fc=15%）情况，土体流失量较多。如前分析，这主要由于裂缝宽度与颗粒粒径之比较大时，颗粒在裂缝处不易堵塞，更易流动导致的。

进一步对比Fc=15%和35%两种工况下的支护侧地表沉降演化，可以发现，地表沉降都在靠近挡墙处最大，且在管涌启动时刻（t=0.2s）不明显，而在t=0.4～0.8s 时刻内发展较快，随后沉降增量趋于减小。对于Fc=15%的情况，由于其支护侧土体流失量较小，地表沉降范围主要局限在0

然而，开挖侧坑底土体隆起演化结果与支护侧地表沉降演化不同。在管涌启动时刻（t=0.2s），坑底隆起即达到最大值，随后下降并由于支护侧土体流入而最终略有上升。这一现象与图3中开挖侧土体在管涌启动时刻即具有较大速度相吻合。这说明坑底隆起主要是管涌流对坑底土颗粒向上的拖曳力，使得土体悬浮于管涌流中引起。

5 结论

本文介绍了一种可反映流体与砂土颗粒间相互作用的CFD-DEM 方法，并利用该法对深基坑开挖过程挡墙裂缝导致的管涌问题进行了仿真分析。通过砂土基坑管涌的宏细观仿真结果，可以发现CFD-DEM 方法可以较好重现基坑管涌的主要特征，如坑底隆起、地表沉降和管涌通道中流体流动特征。管涌过程中地表沉降逐步增大，而坑底隆起在管涌启动时即达到最大值。管涌灾害主要由渗流力对基坑开挖侧土体力学性质的削弱引起，在实际工程中为了防止管涌灾害的发生，除提高挡墙或止水帷幕施工质量外，仍需对开挖侧土体实施良好的加固措施。