考虑动态流体网格的颗粒-流体耦合算法

何金辉，李明广，陈锦剑，夏小和

(上海交通大学 土木工程系， 上海 200240)

由于天然土具有散粒性及含水性，采用微观流固耦合方法进行岩土材料不排水特性分析已成为越来越多学者的选择.目前，微观耦合法主要有基于离散-连续耦合理论的管域法和粗网格法[1].管域法按照初始颗粒的分布划分流体域，不同流体域通过管道进行流体交换[2].该方法简单直观，易于实现，近年来得到了广泛应用.Zhang等[3]采用管域法模拟了黄土中裂隙注浆及其发展过程，并同模型试验结果进行了对比分析，发现最终都会形成“Y”形注浆裂缝.Zeng等[4]通过管域法研究了横向各向同性岩体中水压致裂的机理，并对注水速率、地应力比和倾角等参数进行了敏感性分析.Wang等[5]在管域法基础上，根据流动守恒原则重新建立了流体压力更新准则，并通过稳态饱和介质渗流验证了方法的可靠性.尽管管域法可简化耦合模拟过程并获得较为准确的结果，但流域结构的分布只与颗粒初始位置有关，而不会根据颗粒位移进行更新，因此仅适用于小变形的情形.

粗网格法对颗粒群划分初始流体网格，通过求解N-S方程来更新流体信息并计算流体对颗粒的作用力[6].因计算方便，许多学者采用该方法研究宏观现象机理，推动了相关理论的发展.如Zeghal等[7]采用粗网格法研究了饱和无黏性沉积土在动态激励作用下的液化机理，揭示了若干沉积质液化的微观力学机理和响应模式.Tao等[8]通过粗网格法对管涌的微观力学机制进行了研究，弥补了相关理论的空白.王学红[6]通过粗网格法研究了振动桩下沉的微观力学机理及环境响应规律，对完善相关理论、指导现场施工具有重要意义.粗网格法中的流体网格始终是固定的，不会随颗粒群边界的变化而移动，因此在大变形案例模拟中具有一定局限性.

随着微观耦合方法应用场景的扩大化、复杂化，考虑大变形下流固边界的匹配性具有重要的研究价值和现实意义.目前，已有学者基于传统耦合法进行了改进，考虑了流固边界相适应的因素，并进行了验证.如Zhang等[9]通过把动网格节点的速度引入N-S方程，从而考虑动边界的影响；Zhang等[10]在研究潜蚀问题时，根据大颗粒初始分布建立四面体流体网格，后续模拟过程中结合颗粒位置进行网格重新划分，从而实现流体网格的实时更新.当前动网格方法虽解决了部分应用场景下流固边界不匹配的问题，但在推广过程中仍存在一定局限性.比如，N-S方程迭代运算效率低，方程因简化而未考虑部分项(比如流体-颗粒相互作用力)，四面体流体网格算法空间搜索耗时长，应用场景受限制(需大小颗粒共存).因此，亟需一种原理简单、易于实现、适用范围广的考虑动边界的微观双向耦合算法.

本文在离散元商业软件PFC2D的基础上，引入流体网格动态更新方法，推导动网格下的达西渗流方程和颗粒-流体相互作用方程，并通过Python语言和C语言混编的方法将算法嵌入软件，开发了可考虑流体动边界的微观耦合模块.将该算法用于模拟饱和土不排水剪切双轴试验，并通过和常体积法的结果对比验证了算法的有效性.该算法扩大了微观耦合法的适用范围，在大变形算例的模拟中有助于提高仿真精度，同时通过技术优化，在一定程度上解决了运算效率低的问题，因此对微观尺度上的大变形流固耦合模拟具有一定的参考价值.

1 颗粒-流体耦合算法原理及程序 实现

1.1 颗粒场分析

在不考虑流场作用下，颗粒一般受到体力、阻尼力及粒间相互作用力的影响.本文考虑了流体相和颗粒相之间的双向耦合，因此，颗粒还受到流体对其施加的作用力.在流场中，颗粒一般受到流体静态作用力和流体动态作用力的影响[11].本文采用该耦合算法模拟双轴试验，因加载过程中流体和颗粒速度均较慢，因而可以忽略流体动态力的影响.根据牛顿第二定律，颗粒受力情况可以表示为

(1)

(2)

1.2 流体场分析

从流体角度进行分析时，主要考虑体积应变引起孔压增量和不同孔隙之间发生渗流作用两个因素.每经历一个耦合间隔时，对流体网格位置及流体信息进行更新，然后返回对颗粒的作用力矩阵.本文参考文献[9]中的动网格实现方法，对网格的速度和位置进行更新，如下式：

(3)

(4)

对于更新之后的网格及此时的颗粒分布状态，根据文献[12]中孔隙结构变化引起孔压累积的原理方法进行孔压增量计算.首先对各个网格有：

(5)

(6)

(7)

(8)

(9)

根据式(7)和(9)，预测相对速度vpre和实际相对速度vrel间的误差δ可以表示为

(10)

(11)

(12)

Δp=Bfεv

(13)

式中：Bf为流体体积模量.在空间流体网格之间，由于流体压力梯度的存在，会发生流体的渗流作用.为了避免繁重的计算任务，本文不考虑采用N-S方程进行模拟，而是采用达西渗流定律进行计算.如图1所示，每个流体单元在压差作用下，都会与周围其他单元发生流体交换(I，J表示流体网格在二维空间内的位置坐标，取1，2，…).任意流体单元因压力梯度而与相邻单元发生流体交换(流入/流出)的速度可以表示为[13]

图1 某网格与相邻网格间渗流作用示意图Fig.1 Schematic diagram of seepage action between a certain mesh and adjacent meshes

(14)

式中：g=1, 2表示二维空间某方向；h=1, 2表示g方向坐标轴的负向或正向；vgh为某网格g方向上h侧的边界渗流速度；K和Kh分别为该网格和g方向h侧相邻网格的渗透系数；p和ph分别为上述相邻两网格的孔隙水压力；Ug和Ugh分别为相邻两网格中心位置g方向上的坐标分量；γw为流体的容重.

因为网格是时刻变化的，所以需要根据网格的速度对流体交换速度进行修正，获得真实的渗流流速.结合式(3)，(14)，可得边界修正渗流速度为

(15)

(16)

(17)

式中：Vcell为网格容积.由上可得新的孔隙水压力pnew为

(18)

式中：pold为该流体网格上一时刻的总孔压.

1.3 两相耦合场分析

如1.1节所述，本文只考虑流体对颗粒的静态作用力.如图2所示，当颗粒完全浸入网格时，颗粒处于平衡状态；而当颗粒处于边界位置时，颗粒处于不平衡状态，将受到流体对其的作用力.由此可知，流体网格内某不平衡状态颗粒所受的静态作用力为[14]

图2 流体网格中颗粒位置及受力情况示意图Fig.2 Schematic diagram of positions of particles and their stress states in a fluid mesh

fstatic=pnewAcrn

(19)

式中：Acr为流体网格内流场对不平衡状态颗粒的等效作用面积；n为颗粒所受流体作用力方向上的单位向量.

1.4 颗粒-流体耦合算法程序实现

该耦合算法的实现思路为：① 进行力学计算，在未达到程序终止条件的前提下，当该力学进程时间Δtcur等于设定的耦合间隔时间Δt时，根据颗粒边界的速度和位移，对网格位置及节点速度进行更新.② 在更新后的流体网格下，求解因颗粒运动引起的体积应变，并更新孔隙水压力.③ 根据动网格下的修正达西渗流方程对流场进行更新.④ 判断流场内的颗粒浸润状态，计算流体对颗粒的作用力，并以外力形式施加于各个颗粒，然后进入下一个循环.具体流程如图3所示.

图3 耦合流程示意图Fig.3 Schematic diagram of coupling process

2 颗粒-流体耦合算法验证

2.1 试样制备

首先在宽度为10 cm、高度为20 cm的容器内，按照初始孔隙率0.16生成试样，生成的颗粒数目为 1 004.试样的颗粒级配曲线示意图如图4所示.图中：dp为颗粒直径，w为试样中直径小于某数值的颗粒所占的质量分数.采用线性接触本构模拟密实粒状土不排水剪切特性，本构微观参数取值如表1所示.

表1 颗粒参数Tab.1 Particle parameters

图4 颗粒级配曲线图Fig.4 Diagram of grain size distribution

试样生成后，对材料施加100 kPa的有效应力，从而实现试样的等向固结.之后，保持100 kPa的围压，在耦合条件下进行不排水双轴试验，其中流体模量为2.0 GPa，密度为 1 000 kg/m3.双轴试验模型示意图如图5所示，加载过程中对上下刚性墙体施加5 mm/s的轴向加载速度，同时通过伺服机制保持侧向围压恒定，记录加载过程中偏应力、超孔隙水压力及轴向应变等用于后续的对比和分析.

图5 双轴试验模型示意图Fig.5 Schematic diagram of biaxial test model

2.2 常体积法原理

虽然常体积法不能考虑土水之间相互作用机制，但因其计算效率高、模拟结果基本能反映不排水剪切规律，常被用于耦合算法的验证.

当采用常体积法时，通过控制试样体积不变来模拟材料的不排水行为[15].因此，当轴向加载速度为va时，计算得到的径向速度vr应保证试样在任意时间间隔Δτ内体积保持不变，即体积应变εv为0，具体为

(20)

式中：Wx，Wy分别为容器x，y方向的长度.由式(20)可以得出：

(21)

2.3 结果对比验证

图6所示为常体积法和耦合法的超孔隙水压力对比图，图中：εa为轴向应变，pexc为侧墙附近产生的超孔隙水压力.可以看出，两种方法下的孔压规律较为一致，数值比较接近.初始时刻均生成正孔压，然后试样发生体积剪胀，孔压逐渐消散，变为负孔压.该趋势与较密实砂土的不排水剪切行为一致.

图6 耦合法和常体积法下的超孔隙水压力对比Fig.6 Comparisons of excess pore water pressure in coupling method and constant volume method

图7所示为两种方法下的偏应力曲线，图中q为偏应力.由图可得，轴向应变为4%以内时，两种方法下的偏应力数值基本没有差别.继续加载，轴向应变超过4%时，二者数据虽有差异但偏差不大，均表现为随着加载过程而偏应力不断增大的趋势，同时轴向应变在10%附近，两者趋于相近.

图7 耦合法和常体积法下的偏应力结果对比Fig.7 Comparisons of deviatoric stress results in coupling method and constant volume method

事实上，常体积法理论中流体是不可压缩的，也未考虑土水的相互作用，通过间接方式考虑流体信息的计算和更新，这种简化使得其模拟结果与实验数据不能高度吻合，但其结果显示的规律性同实验较为一致.从图6、7可以看出，本文提出的微观耦合法同常体积法数值比较接近，模拟规律一致，都体现出密实砂土在不排水剪切时的“硬化”趋势.同时该耦合法结合了流固耦合及动网格技术，模拟过程更为贴合实际，因此可为类似研究提供参考.

3 不同围压影响下的双轴试验结果 分析

为了进一步验证方法的有效性，采用耦合法模拟了围压分别为100，200，300 kPa时的饱和土不排水剪切双轴试验，并对偏应力、孔隙水压力及有效应力比的变化情况进行了对比分析.

图8所示为不同围压下的偏应力变化情况.从图8可以看出，当围压增大时，相同应变时偏应力更大.原因为大的围压对土体的侧向膨胀限制作用更加明显，颗粒间的摩擦力更大，表现为土体模量较大，因此剪切强度也就越大.

图8 不同围压下偏应力结果对比Fig.8 Comparison of deviatoric stress results at different confining pressures

图9所示为不同围压下的超孔压发展情况.可以看出，不同围压下试样的超孔压变化趋势均为先升高后下降，表明3组试样在应变较小时孔隙受到压缩，产生正的孔隙水压力，随着加载的进行，试样发生了剪胀，孔压逐渐降低，出现负值.对比3组试样的孔压曲线可知，当围压增大时，正孔压幅值增大，峰值点延后，负孔压值较小.因为随着围压的增大，侧墙对试样剪胀的抑制作用更明显，所以试样受压缩程度较高，因而正孔压的累积作用也就较强.

图9 不同围压下的超孔隙水压力结果对比Fig.9 Comparisons of excess pore water pressure results at different confining pressures

图10所示为不同围压下有效应力比η′的变化趋势图.三者有效应力比峰值均出现在轴向应变为1%左右，之后应力比逐渐降低.在1%附近时，100，200及300 kPa围压(后文默认为此顺序)对应的有效摩擦角分别为31.9°，31.3°及30.3°.3组试样在轴向应变9%附近的摩擦角差别相对较大，分别对应为18.0°，23.1°，22.3°.由图8可以看到，100 kPa围压的试样在轴向应变8%～9%时偏应力有一个突变，此时可能为该试样局部力链发生断裂，颗粒随机重新分布.而在加载终止时，三者的有效应力比趋于一致，其对应的摩擦角分别为20.3°，22.0° 及21.4°.综上，围压变化对摩擦角影响不大，这与文献[13]的研究结论一致.

图10 不同围压下有效应力比结果对比Fig.10 Comparisons of effective stress ratio results at different confining pressures

4 结论

本文研究了一种考虑动态流体网格的颗粒-流体耦合算法，将算法嵌入离散元商业软件PFC2D中，在大变形微观耦合算例模拟中可提高仿真精度.将该算法用于模拟饱和土不排水剪切双轴试验，并通过和常体积法的结果对比验证了算法的有效性.最后，采用该算法研究了不同围压影响下粒状土双轴剪切特性，模拟规律与室内试验结果较为吻合.主要结论如下：

(1) 在引入流体网格动态更新方法的前提下，改进并推导动网格下的达西渗流方程和颗粒-流体相互作用方程，在离散元商业软件PFC2D基础上开发了考虑动态流体网格的颗粒-流体耦合模块.该算法原理简单，易于实现，弥补了传统微观耦合方法的不足，为大变形微观耦合模拟提供了参考，但在流场可视化、算法效率等方面仍需进一步的改进和研究；

(2) 将算法用于模拟饱和土不排水剪切双轴试验，并通过和常体积法的结果对比，验证了算法的有效性.结果显示，两种方法模拟得到的偏应力曲线和孔压曲线规律一致，数值较为接近.