一种CFD-DEM 流固耦合方法在渗流导致城市地面沉降问题中的应用

2020-12-10李晓蛟武亚军

上海大学学报（自然科学版） 2020年5期

李晓蛟，陆烨，武亚军

(上海大学土木工程系，上海200444)

城市局部区域地面沉降会对周边构建筑物及使用人群造成影响，带来巨大经济损失，而地下水渗流是造成地面沉降的主要原因. 对于城市局部区域地面沉降这一问题，现场监测有较大难度，且现场外部环境随时间变化，影响测量结果的因素较多. 模型试验一方面难以揭示地面沉降的细观机理，另一方面也不能得到土体全场位移规律. 因此，数值模拟方法可作为现场监测和模型试验的补充，对沉降过程及机理进行深入研究.

在以往的研究中，对于城市地面沉降的数值模拟主要采用有限元方法(finite element method， FEM)[1-2]或离散元方法(discrete element method， DEM). 有限元方法能够较好地模拟宏观尺度下连续介质的变形，但是在模拟细观尺度下离散介质的运动规律时有所不足，即有限元方法要满足变形协调条件，其内部细观结构的尺寸要远小于模拟对象的尺寸[3]，这对于研究渗流作用下颗粒的迁移规律显然不再适用. 而Cundall等[4]提出的离散元方法则非常适合用于模拟细观尺度下土颗粒的迁移运动，并应用于多种岩土问题[5-8]. 但是，当模拟对象中涉及渗流或孔压时，如何处理流体是需要解决的问题. Li等[8]曾提出一个简化的流体-颗粒耦合方法，该方法中颗粒之间的孔隙被认为是通过颗粒接触之间的管子连接的流域，通过建立网格来表征流域的体积以及管子的长度和直径. 这种方法对少量颗粒系统比较简便，但在三维情况下，孔隙体积计算将变得繁琐.

计算流体动力学(computational fluid dynamics， CFD)技术是目前广泛使用的模拟流体的计算方法，不仅可以进行三维流场的计算，还可以模拟复杂形状的流场[9]. Zeghal等[10]首次在岩土工程问题的研究中引入了CFD-DEM 流固耦合方法，分析了土坡渗流和饱和土体震动液化问题，取得了较好的结果. 周健等[11]和罗勇等[12]也运用CFD-DEM 方法开展了土体渗流等问题的研究. 蒋明镜等[13]利用PFC2D，借助Fish 函数进行二次开发研究流固耦合问题，通过单颗粒水下自由沉降和一维固结试验算例来验证CFD-DEM 耦合的可行性，取得了较好的结果. 王胤等[14]提出了一种考虑粒间滚动阻力的CFD-DEM 模拟方法，能够较好地反映颗粒的形状效应，建立了流体与固体颗粒的动量传递模型，能够很好地描述岩土问题中水-土结合的力学性质，但该方法的应用范围具有一定局限性，例如DEM 模型不能识别具有复杂形状的CFD 模型. 从以上研究可以看出，已有的CFD-DEM 流固耦合模拟方法虽然实现了二维到三维的跨越，在岩土工程中的应用也越来越多，但大部分研究工作基于结构化网格，所模拟的流场形状十分规则单调(如方形或圆柱形)，对于大型的、具有复杂形状的实际工程，模拟起来仍有一定难度.

鉴于此，本工作采用一种 CFD-DEM 方法，即颗粒流程序 (particle flow code， PFC)与计算流体力学程序(Fluent)流固耦合方法，通过单颗粒在水中自由沉降算例验证本方法的正确性；依据本方法建立具有复杂流场形状的基坑临空面渗漏与河道边坡渗流两个算例，在设置不同地下水位的前提下，通过Fluent 求解渗流场，然后将渗流场导入PFC 土颗粒模型中，实现耦合计算，最终模拟出不同水力梯度作用下渗流导致地面沉降的整个动态过程. 为便于区分，本工作提出的CFD-DEM 方法在下文中统一称为PFC-Fluent 方法.

1 PFC-Fluent 计算理论基础

本工作中的流固耦合程序采用PFC3D 模拟固相颗粒的细观行为，采用均化流体CFD 技术模拟液相流体的运动. 对于固相颗粒，采用牛顿第二定律求解颗粒的平移和转动方程; 对于液相流体，通过求解平均的Navier-Stokes 方程模拟孔隙流体的运动; 流体对颗粒的作用力由Ergun 方程表达[15]. 每个流体单元包含一定数量的土颗粒，流体流动产生的渗透力以体力的形式作用于土颗粒.

PFC-Fluent 耦合程序是单向耦合，对于地下水渗流导致的地面沉降问题，土体在地下水渗流力作用下被冲刷破坏，而地下水位基本保持不变，流体对土颗粒的影响远远大于土颗粒对流体的影响，因此在这种情况下，采用PFC-Fluent 耦合程序是可行的.

1.1 固相运动方程

固相颗粒的运动由颗粒离散元来模拟. 作用在单个颗粒上的主要有体力(颗粒的自重)、其他与之接触的颗粒通过接触点施加的接触力，以及流体作用在该颗粒上的力. 在这些力的作用下，颗粒的运动满足牛顿第二定律，平动和转动方程为

式(1)和(2)中：u 为土颗粒速度矢量; t 为时间; fmech为颗粒重力与接触力的合力; ffluid为流体对单个土颗粒的合力; m 为土颗粒质量; g 为重力加速度; ω 为角速度矢量; M 为力矩; I 为惯性矩.

1.2 液相运动方程

对于密度不变不可压缩流体的运动，可以由平均的Navier-Stokes 连续方程和动量方程来表达，即

式(3)和(4)中：n 为孔隙率; ∇为哈密顿算子; v 为流体速度矢量; ρf为流体密度; ∇p 为压力梯度; τ 为黏性应力张量; fd为单位体积内流体与颗粒之间的作用力.

流体网格内孔隙率n 的计算方法有两种：质心法和多面体法. 本工作采用质心法，即当颗粒球心位于网格内时即识别为整个颗粒都处于网格内. 质心法的优点是计算效率高，缺点是随着颗粒的运动，孔隙率会产生跳跃，因此在模拟中设定的网格尺寸相对颗粒来说足够粗，以削弱这种影响，相对于整个流场，网格又是足够细的，这就在保征准确性的前提下提高了计算效率.

1.3 流体对土颗粒的作用力

在固-液两相饱和介质中，流体对颗粒的作用力从大类上可分为静水力和动水力. 静水力即为颗粒静置于流体中的浮力作用，流体对颗粒的浮力大小与该颗粒周围的流体压力梯度有关，表达式为

式中：fb表示单个颗粒受到的浮力; r 表示颗粒半径.

动水力包括拖拽力、上举力和虚拟质量力，后两种力的大小和拖拽力相比可以忽略. 拖拽力的表达式为

式中：fdrag为土颗粒受到的拖拽力; Cd为拖拽力系数，适用于整个区域的近似公式

其中Re 为雷诺数，

μ 为黏滞系数.

2 PFC-Fluent 耦合流程与程序实现

采用PFC 软件求解固体颗粒运动方程，用Fluent 软件求解流体动力方程. 将Fluent 求得的三维流场导入PFC 程序中，将颗粒与流体相互作用力作为耦合纽带，实现相应的流固耦合计算. PFC-Fluent 耦合计算流程如图1 所示.

图1 PFC-Fluent 耦合方法流程图Fig.1 Flow chart of PFC-Fluent coupled method

3 数值模拟验证

首先，通过3 个不同尺寸颗粒在水下自由沉降的算例，验证PFC-Fluent 耦合程序的准确性. 然后，通过不同水位差的两组算例来验证PFC-Fluent 耦合计算的可行性. 两组数值算例如下：基坑开挖时，围护结构出现渗流导致坑外地层沉降; 河道边坡在渗流作用下发生破坏.

3.1 算例一 3 个不同尺寸颗粒水下自由沉降

3.1.1 自由沉降速度理论解

斯托克定律表明，小球在水中重力作用下的沉降速度是一定的. 采用PFC-Fluent 耦合程序模拟3 个不同尺寸颗粒在水中重力作用下的自由落体运动，以验证本方法的正确性. 在黏性流体中，球形颗粒在重力作用下下落的运动方程如下：

当下落的球体达到一个稳定的速度，加速度为0，式(9)变为

3.1.2 PFC-Fluent 数值模拟

将半径为1.0， 1.5， 2.0 mm 的颗粒置于黏滞系数为0.001 的流体中，在重力作用下自由下落. 颗粒与液体的密度分别为2 650， 1 000 kg/m3，重力加速度为9.8 m/s2. 将上述参数代入式 (7)， (8)， (10)，并用 Matlab 数值求解得到自由沉降速度理论解分别为 0.245， 0.324，0.391 m/s.

耦合程序中颗粒速度变化过程如图2 所示，最终速度达到稳定值，与理论解一致，表明本工作采用的耦合程序可以实现对颗粒施加流体作用.

图2 不同时刻颗粒下落速度Fig.2 Drop velocity of particle at different time

3.2 算例二基坑临空面渗漏

本算例考虑基坑开挖完成后还未浇筑底板时，由于降水坑内外水位差较大，同时基坑围护结构在临近坑底处出现裂缝，水土在水力梯度作用下向基坑临空面渗漏，导致坑外土体流失从而产生地面沉降.3.2.1 数值建模

基于颗粒流软件PFC3D 平台建立基坑边坡土体模型的基本流程如下.

第 1 步，采用 wall 命令生成一个尺寸为 6.0 m×1.0 m×9.0 m(长×宽×高)的模型槽，在模型槽内按线性接触模型生成高度为6.5 m，孔隙率为0.4 的砂土层. 土体通过一系列球形颗粒来模拟，对所有土颗粒施加重力，在重力作用下颗粒沉积、固结，并达到平衡状态. 在砂土层上部按平行黏结点接触模型生成高度为1.5 m，孔隙率为0.4 的黏土层，施加重力，同样达到平衡状态. 基坑模型中颗粒总数为67 721 个，颗粒最大半径为0.06 m，最小半径为0.02 m，平均半径为0.04 m，颗粒与墙体间的摩擦系数取值为0，以减小边界效应的影响. 模型中参数取值参考类似地质条件[16]，如表1 所示.

表1 PFC3D 模型细观参数Table 1 Mesoscopic parameters of PFC3D model

第2 步，由于本工作主要考虑临空面渗漏对地面沉降的影响，不考虑基坑围护结构变形的影响，故采用wall 命令在距模型槽右侧边界2.5 m 处生成一面无厚度墙作为基坑的围护结构，围护结构总长7 m，删除围护结构右侧高度为3.5 m 区域内的土颗粒，围护结构插入比为1. 针对实际工程中基坑围护结构产生裂缝的情况，在围护结构离坑底0.5 m 处开一圆孔. 若圆孔半径太小，土颗粒会堵塞裂缝，故经过不断试算，得到圆孔直径为5 倍土粒平均粒径可引起渗流，因本工作重点不对圆孔尺寸进行对比，故只取这一种尺寸进行计算. 基坑土体模型如图3 所示.

第 3 步，打开 PFC 流体模块，将 Fluent 软件计算出的流场 (3.2.2 节详述)导入 PFC 软件中，即可进行流固耦合计算.

3.2.2 流场计算

对于多孔介质的孔隙率及相关阻力系数，孔隙率取值与土体孔隙率一致; 对于土中水渗流这种流速较低的层流状态，阻力系数只考虑黏性阻力系数，根据孔隙率与土颗粒粒径计算出黏性阻力系数值[13].

本算例模拟了3 种不同地下水位，其中坑底地下水位保持不变，均为坑底地面以下0.5 m，坑外地下水位依次取地面以下1.5， 0.5 和0 m. 流场进口采用压力进口边界，进口压力为0 kPa，出口采用压力出口边界，出口压力为0 kPa，其余边界设置为墙边界，水在重力作用下流动. 计算出的流场如图4 所示，将其导入PFC 建立的土体模型中进行耦合计算.

图4 基坑开挖流场流速矢量图Fig.4 Velocity vector diagram of flow field for excavation

3.2.3 模拟结果

在算例二中，造成地面沉降的原因有两方面，一方面是土颗粒在地下水渗流作用下从基坑围护结构破损处流出(流出的土颗粒被程序自动删除)，继而上部土体往下运动，最终造成坑外地面沉降; 另一方面是土颗粒在水力梯度作用下从高水位向低水位运动，导致坑外地面沉降.本算例分别模拟了3 种不同地下水位，通过不同时步土颗粒位置坐标变化得出地表沉降位移如图5 所示. 从图中可以看出，随计算时步的增加，地表沉降值越来越大; 基坑边缘处地表沉降值最大，沉降可达0.35 m; 坑外水位越高，地表沉降值越大. 水位高度在地面以下0.5 m 与在地面以下0 m 的地表沉降比较接近，原因可能是二者水力梯度比较接近.

图5 基坑临空面渗漏引起的地表沉降Fig.5 Ground settlements caused by leaking on retaining structure of excavation

在3 种不同水位高度下，土颗粒位移云图比较一致，故仅选取水位高度在地面以下0 m 这一种情况下的不同时步的土颗粒位移云图(见图6)，图中各颜色表示土颗粒位移大小，蓝色表示位移为0 m，红色表示位移为1.0 m 及以上，中间每隔0.1 m 用不同颜色表示. 模拟结果表明，随着时步不断增加，基坑围护结构开孔处土颗粒不断流失，上部土体不断往下移动，最终形成地表沉降. 同时，基坑底部在渗流作用下有轻微隆起.

图6 水位高度在地面以下0 m 时的颗粒位移Fig.6 Soil particle displacements for underground water level of 0 m

3.3 算例三河道边坡渗流

本算例考虑河道边坡在水力梯度作用下发生渗流破坏诱发边坡外侧的地面沉降.

3.3.1 数值建模

基于颗粒流软件PFC3D 平台建立河道边坡土体模型的基本流程如下.

第1 步，采用wall 命令生成一个尺寸为12.5 m×1.0 m×9.0 m(长×宽×高)的模型槽，在模型槽内按线性接触模型生成高度为6.5 m，孔隙率为0.4 的砂土层. 土体通过一系列球形颗粒来模拟，对所有土颗粒施加重力，在重力作用下颗粒沉积、固结，并达到稳定状态. 在砂土层上部按平行黏结点接触模型生成高度为1.5 m，孔隙率为0.4 的黏土层，施加重力，达到稳定状态.模型槽中颗粒总数为131 308 个，土颗粒尺寸与算例二一致. 颗粒与墙体间的摩擦系数取值为0，以减小边界效应的影响. 模型中参数取值如表1 所示.

第 2 步，河道边坡是一个倾斜角为 45°，长 6.4 m 的护坡墙，由 2 000 个半径为 0.1 m 的圆球按平行黏结接触模型生成. 借助Fish 语言删除边坡右侧高度为3.5 m 内的土颗粒(见图7).

第 3 步，打开 PFC 流体模块，将 Fluent 软件计算出的流场导入 PFC 软件中，即可进行流固耦合计算.

3.3.2 流场计算

河道边坡模型的流场同样用Fluent 进行计算，由于算例三与算例二的土体颗粒级配、孔隙率一样，故本算例中流体网格尺寸和黏性阻力系数与算例二一致. 本算例模拟了河道水位的3 种高度，其中地下水位保持不变，为地面以下0.5 m，河道水位高度依次取为2， 1 和0 m. 流场进口采用压力进口边界，进口压力为0 kPa，出口采用压力出口边界，根据水位深度，出口压力依次取为0， 9.81， 19.62 kPa，其余边界设置为墙边界. 计算出的流场如图8 所示，将其导入PFC 建立的土体模型中进行耦合计算.

图 7 河道PFC 模型Fig.7 PFC model for river embankment

图8 河道渗流流场流速矢量图Fig.8 Velocity vector diagram of flow field for river embankment

3.3.3 模拟结果

算例三的地表沉降位移如图9 所示. 当河道水位高度为2 m 时，运行80 万步后地表沉降位移最大值为0.34 m; 当河道水位高度为1 m 时，运行80 万步后地表沉降位移最大值为0.58 m; 当河道水位高度为0 m 时，运行80 万步后地表沉降位移最大值为1.14 m. 可见，边坡两侧水位高度相差越大，地表沉降值也越大; 距河道坡顶越近，地表沉降值越大.

图9 河道边坡渗流引起的地表沉降Fig.9 Ground settlements caused by seepage under river embankment

在3 种不同河道水位下，运行80 万步时不同水位土颗粒位移云图如图10 所示(图中颜色表示意义与图6 一致). 模拟结果表明，河道边坡两侧水位高度相差越大，土颗粒运动越迅速，地表沉降值也越大.

图10 河道边坡渗流引起的土颗粒位移Fig.10 Soil particle displacements caused by seepage under river embankment

河道边坡渗流算例主要用来观察在初始流场作用下渗流对土颗粒位移的影响，最终模拟结果符合河道渗流规律. 河道边坡两侧水位高度差和渗透系数是影响流场分布的重要因素，在本算例中二者无明显变化. 土颗粒虽然发生了较大位移，但仍然处于现有流场范围内，且流体网格尺寸相对颗粒尺寸是足够粗的，当有土颗粒流出时，周围其他土颗粒也会流入，因此孔隙率变化并不是特别明显. 本算例主要研究的是固体颗粒的位移规律，流体对颗粒的影响远大于颗粒对流体的影响，因此采用PFC-Fluent 是有一定合理性的.