基于近场动力学与有限体积法耦合的裂隙岩体渗流模拟

李卓徽，周宗青，高成路，张道生，白松松

（1. 浙江大学建筑工程学院，浙江 杭州 310058；2. 山东大学齐鲁交通学院，山东 济南 250002；3. 陆军勤务学院军事设施系，重庆 401311）

近几十年来介质中的流固耦合数值研究已成为力学、土木和环境工程等领域的主要研究课题，该研究涉及到固体变形和裂隙演化扩展以及裂隙区流体流动。为了求解这类问题，出现了各种复杂的数值模型，可简单分为3 类：连续、离散和混合方法。典型的连续方法主要是基于有限元法，但此类方法难以解决裂纹尖端应力场奇异问题；离散方法通常采用离散元法，在宏细观参数标定方面存在复杂性和局限性；为了结合连续方法和离散方法的优点，一些混合方法被提出［1］，但这些方法均在处理多裂纹、裂纹分叉问题时仍面临挑战，特别是在三维条件下。

近场动力学（PD）理论最早是由美国Sandia 国家实验室的Silling［2］在2000年提出。其基本思想是以非局部作用的积分模型代替传统理论的微分模型，可以解决有限元方法中尖端奇异性的问题，同时，兼具分子动力学方法和无网格方法的优点，又突破了经典分子动力学方法在计算尺度上的局限［2］。这使得近场动力学理论可以描述从连续到不连续、从微观到宏观的一系列力学行为。Zhou等［3］提出了一种考虑压缩荷载作用下岩体峰后阶段的微观弹塑性本构模型。Gao 等［4-6］将近场动力学应用于隧道开挖，分析了在不同强度折减系数下不同节理倾角的节理岩体的破坏模式，以及模拟了在隧道开挖过程中岩体渐进破坏导致突水通道形成的演化过程。与此同时，近场动力学理论还能模拟流体渗流问题，在统一框架下实现流-固耦合模拟是近场动力学理论的又一大突出优势。因此，越来越多的学者将近场动力学理论用于岩石材料流固耦合问题的模拟研究中。寿云东［7］针对岩土工程中的温度场-渗流场-应力场耦合问题，引入最大拉应力强度准则、莫尔-库仑强度准则和双剪强度准则，建立了基于近场动力学理论的裂隙岩体温度场-渗流场-应力场耦合数值模型，实现了巷道围岩失稳破坏过程的模拟。王允腾［8］专注键型近场动力学理论研究，针对岩体热-水-力-化（THMC）耦合作用机理，提出了THMC 耦合作用下的键型近场动力学数值模型，从而实现了THMC 耦合作用下岩体破裂过程模拟。然而，由于近场动力学涉及到非局部的相互作用，与经典模型相比，近场动力学理论在离散化时需要更多的计算资源，这一点在模拟流固耦合问题时尤为明显。为了在合理的计算时间内模拟各种多尺度物理过程，国内外学者提出了近场动力学与经典模型耦合的不同方案。Macek和Silling［9］利用嵌入节点和单元将PD 网格与有限元（FEM）网格耦合。Kilic 和Madenci［10］在PD 域和有限体积（FV）域的界面处引入了一个重叠区域，在 该 区 域 内 求 解2 个 方 程。 Agwai 等［11］和Oterkus［12］采用了FEM 进行整体建模，又采用PD进行预测材料失效的子建模。Lubineau 等［13］采用基于能量等效的变形函数进行耦合。这些函数只影响本构参数，因此允许模型作为纯粹非局部的、纯 粹 局 部 的 或 混 合 的。Liu 和Hong［14］将FEM 和PD 域与界面元耦合，设计了2 种不同的耦合方案。Ni 等［15］提出了一种基于PD-FEM 耦合的模拟饱和多孔介质中水力裂缝扩展的混合建模方法。可以看出利用近场动力学和连续介质理论耦合的方法在不降低求解精度的情况下有效地减少计算量已经成为了目前研究的一大热点问题。

受Ni 等［15］启发，提出一种有限体积法与近场动力学耦合的方法来模拟饱和孔隙裂隙介质中水力裂缝扩展问题，以此提高计算效率，并开展室内试验尺度下多工况岩体水力压裂模拟，验证该方法在流体驱动下模拟饱和裂隙孔隙介质中裂纹扩展的能力。

1 基本原理

1.1 普通态型近场动力学

近场动力学根据物质点间相互作用的不同分为：键型（BB-PD）、普通态型（OSB-PD）和非普通态型（NOSB-PD）3 种不同类型［16-19］。键型近场动力学模型中对点力只考虑物质点间相对位置，力的大小和方向均相同，而在描述物体力学性质的过程中，键型近场动力学最棘手的问题在于材料的泊松比固定，即二维情况下所描述的材料泊松比恒等于1/3［19］；三维情况下所描述的材料泊松比恒等于1/4［19］，这极大地限制了近场动力学在工程领域的应用范围。

Silling等［17］在2007年提出了态型近场动力学理论。在该理论中，Silling 指出计算区域内离散物质点的弹性应变能可分解为各向同性膨胀应变能和形状改变应变能［17-18，20］，提出了“状态”的概念，类似于传统连续介质力学张量概念，并通过引入变形矢量状态和力密度矢量状态使近场动力学能够与传统连续介质力学理论有效结合，从而有效解决泊松比限制问题。

普通态型近场动力学的运动方程为［17］

根据普通态型近场动力学理论模型［12］，可将力矢量状态写为

式中：θ为近场动力学理论中的体积膨胀；-ed为伸长状态的偏张量；k、α为近场动力学材料参数，二者对应于传统弹性理论中的材料体积模量和剪切模量；m为加权体积，m5。

二维平面应变问题下的普通态型近场动力学力密度标量状态表达式为［20］

其中

式中：K和G分别为体积模量和剪切模量；ν为泊松比。

1.2 有限体积法概述

有限体积法（Finite Volume Method，FVM）借鉴了有限元的部分优点，以有限差分法（FDM）为基础发展而来。同FEM 和FDM 一样，将计算区域划分成有限体积大小的离散网格，即进行离散化处理，每个网格节点按照一定的方式形成表征该节点所对应的控制体积，如图2 所示。其关键思想在于针对每一个控制体积将待解的微分方程进行积分形式处理，从而得出一组离散方程。其中的未知数是网格点上的因变量的数值。

图2 有限体积法示意Fig.2 Schematic diagram of finite volume method

有限体积法的基本思路易于理解，具有具体的物理含义。离散方程实际上表示的是控制体积的通量平衡，也就是通量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。有限体积法得出的离散方程使得因变量的积分守恒对任意控制体积都满足，也就是说有限体积法的局部表述可以实现局部守恒，那么整个计算区域自然也得到满足，进而能够以自然、直接的方法来解决对流占主导的流动问题的离散化，则特征变量φ在控制体积内的守恒关系为

式中：φt为φ随时间的变化量；φc为边界对流进入控制体积的量；φd为边界扩散进入控制体积的量；φs为源项产生的量。

对比有限元法，可以发现有限元法的一个缺点是，对于连续试函数和基函数，没有局部守恒，只能保证全局守恒。换句话说，只能保证域边界上的净通量是平衡的。另一个缺点是无法控制局部通量，这意味着稳定对流占主导地位的流动的离散化并不简单。此外，有限体积法推导其离散方程是通过控制体积的积分方程作为出发点，这与有限差分法直接从微分方程推导完全不同。所以，有限差分法仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也能显示出准确的积分守恒。

综上，有限体积法的特点可以总结为以下四点：①控制方程以积分形式表示，这有别于有限差分法。②积分方程体现了局部守恒性，即变量在任意控制体积内守恒，这与有限元法不同。③离散化后，各个离散方程的每一项均有相应的物理解释，这是有限体积法最具优势的地方。④计算区域内离散的各个网格节点之间控制体积不重叠，从而离散的局部变量守恒保证了整个计算区域内场变量的守恒，实现整体守恒。

2 近场动力学与有限体积法耦合理论

2.1 固体模块计算

为了模拟固体变形与孔隙压力变化之间的相互作用，需要考虑流体压力的影响和变形的影响，对近场动力学理论的运动方程进行修正。参考在单相流体饱和多孔介质理论［21-22］中的有效应力原理，其表达式为

式中：σeff为有效应力张量；σtotal为总应力张量；α为Biot系数；p为孔压；I为单位张量。

在近场动力学理论中的标量力密度状态可看作是一个膨胀项和一个偏差项的加和，而且该膨胀项对应于经典连续介质力学理论中的体积应力，利用式（11）所示的有效应力原理，将孔隙压力项直接添加进标量力密度状态的膨胀项中，引入标量有效力密度状态，则标量力密度状态与标量有效力密度状态之间的关系可写成类似于式（11）所给出的总应力张量与有效应力张量之间的关系。其中，-ttotal在三维情况下对应于式（3）中的-t，在二维平面应变情况下对应于式（4）中的-t，在二维平面应力情况下对应于（7）中的-t。

参考文献［15，23-25］思路，将整个计算区域划分为固体域（Rs）、过渡域（Rt）和裂隙域（Rf）三部分。同时，运用近场动力学损伤场作为划分依据，设置2个指标用来标识3个区域，如图3所示。所得线性指标函数为

图3 计算区域划分Fig.3 Division of calculation area

则利用线性指标函数对固体域和裂隙域的渗透率进行插值，可以得到过渡域的固体渗透率

另外，式（12）中的标量有效力密度状态取决于物质点对所处的计算区域，可分为以下3种情况［26］：

（1）物质点对处于固体域（Rs）。物质点对之间的键未断裂，且各物质点邻域内的其余键未出现断裂，即各物质点损伤值为零，则三维和二维情况下的标量有效力密度状态如式（12）所示。

（2）物质点对处于过渡域（Rt）。物质点对之间的键出现断裂，而各物质点邻域内的其余键未出现断裂。此时，式（12）中的物质点与物质点之间的相互作用力项就会消失，但是由于物质点对断裂之后孔隙空间是连续的，因此孔隙压力项是存在的，则三维和二维情况下的标量有效力密度状态表示为

（3）物质点对处于裂隙域（Rf）。裂隙空间的产生需要满足2 个条件：一是物质点之间的键的能量密度超过了临界能量密度，换句话说就是物质点对之间的键断裂；二是各物质点之间的损伤值超过了临界损伤值。通过研究统计可知，当物质点损伤值超过0.4时，裂隙开始生成。故而，三维和二维情况下的标量有效力密度状态表示为

2.2 流体模块计算

用达西定律来描述饱和孔隙裂隙介质中的流场，则可表示为

式中：v为体积流速；μf为流体黏度；P为流体压力。对于均匀各项同性体，K=k I。

在经典理论中，流体计算节点只与相邻节点相互作用，且有限体积单元稳态达西流遵循流量守恒方程。因此，控制方程如式（19）所示：

式中：ρf为流体密度；v为体积流速；S为源项；∇为散度算子。

对于裂隙渗流，采用立方定律来计算裂缝域的渗透率［27-28］

式中：b为裂隙开度。如图1所示，根据有限体积法，网格节点的局部通量守恒方程可以离散为

图1 影响域示意Fig.1 Schematic diagram of area influenced

式中：aP、aW、aS、aE、aN分别为网格节点与相邻网格节点的离散方程系数；φP、φw、φS、φE、φN分别为网格节点与相邻网格节点的特征通量；sP为网格节点源项数值。aP=aW+aS+aE+aN。

对所有网格节点均可列出对应如式（21）的离散方程，对计算域边界处的离散方程按边界条件修正各个系数，再将其写成矩阵形式，可得到

2.3 耦合计算方案

如图4所示，在平面离散中，固体部分由近场动力学物质点离散，流体部分由结构性网格离散。考虑流体压力和固体变形的相互影响，需要在2 个离散网格之间建立过渡层，从而通过过渡层实现流体压力与固体变形之间的相互传递。

图4 耦合过程中的双向影响示意Fig.4 Schematic diagram of bidirectional influence in coupling process

如图4 所示，当流体压力向固体部分传递时，过渡层由近场动力学物质点控制，通过物质点与流体部分网格节点位置来确定物质点所受到的流体压力由哪一个网格节点控制。当固体部分变形向流体网格传递时，此时过渡层由流体部分网格节点所控制，如流体压力向固体部分传递一般，通过确定流体部分网格节点控制体积内所包含的近场动力学物质点，利用这些物质点间的渗透系数向流体网格传递固体变形。算法流程图如图5所示。

图5 算法流程Fig.5 Flowchart of algorithm

3 流固耦合模拟应用与验证

通过3个不同的数值模拟算例说明本文所提出的近场动力学与有限体积法耦合方法适用于模拟岩石材料流固耦合及裂纹扩展破裂过程，揭示相应的裂纹扩展演化力学机制。同时，将该耦合方法的数值模拟结果与其他数值或解析结果相对比，进一步说明所提出的耦合方法能够有效模拟岩石流固耦合的力学扩展破裂过程。

3.1 多孔介质渗流模拟

为了验证近场动力学与有限体积法耦合方法能对基质岩石中的渗流问题进行有效的数值模拟，将基质岩石视作多孔介质，对其渗流过程进行了模拟对比。二维多孔岩石介质模型的几何尺寸如图6所示，将试件简化为一定厚度的矩形，其几何尺寸L、W分别为1.0m、0.2m。在二维多孔岩石介质模型的左右两侧分别施加恒定水压，即PL0=9.5MPa，PR0=4.5MPa。另外多孔岩石介质的渗透率为kr=1.0×10-12m2。

图6 模型的几何尺寸Fig.6 Geometric dimensions of model

在固体近场动力学数值模型中，采用8 000个近场动力学物质点表征相应的多孔岩石介质的计算域。其中，相邻物质点的间距为Δx=0.000 5 m，邻域半径δ=3.015Δx，材料的弹性模量为E=22GPa，密度为ρ=2 462kg·m-3，在此次模拟中可不用设置虚拟物质点作为固体层的边界区域。在流体层设置8 000个网格表征流体的计算域，将固体层中的物质点与流体层中的网格节点一一对应，另流体密度为ρ=1 000kg·m-3，流体黏度系数为μ=1.0×10-3Pa·s。

多孔岩石介质渗流条件下的孔隙水压分布如图7 所示。可以明显从图中看出，多孔介质的水压力从左右2 个边界向中心汇聚，且岩石左侧水压高于岩石右侧水压。为了验证模拟结果的正确性，将多孔岩石介质渗流模拟结果中的水压分布数值模拟结果与解析解［29］对比，如图8 所示。可以观察到数值模拟结果与已有的解析解结果基本吻合，且通过与用近场动力学方法进行渗流模拟的计算时间［8］对比，可以看出近场动力学与有限体积法耦合方法能够高效进行岩石材料相关渗流问题的模拟。

图7 模拟结果云图Fig.7 Cloud chart of simulation results

图8 沿中心轴水压分布PD-FVM数值计算结果与解析解对比Fig.8 Comparison of PD-FVM numerical calculation results and analytical solutions of water pressure distribution along central axis

3.2 流体驱动的裂缝扩展模拟

3.2.1 预制水平裂隙试样

为了验证近场动力学与有限体积法耦合方法能对基质岩石中的渗流问题进行有效的数值模拟，将对含有预制水平裂隙的岩石基质中的渗流过程进行了模拟验证。二维含水平预制裂隙的岩石介质模型的几何尺寸如图9 所示，将试件简化为一定厚度的矩形，其几何尺寸L、W分别为1.0m、1.0m。在该岩石介质模型的中间位置预制一个水平裂隙，该预制裂隙长度为0.2m。将岩石介质模型的4 个边界均设为排水边界，即P0=0MPa，并且在预制水平裂隙内部施加水压p=1×104Pa·s-1。另外多孔岩石介质的渗透率为kr=1.0×10-12m2，而裂隙中的等效渗透率为kf=1.333×10-6m2。

图9 含有预制水平裂隙模型的几何尺寸Fig.9 Geometric dimensions of model with a preexisting horizontal flaw

在固体近场动力学数值模型中，采用10 000个近场动力学物质点表征相应的多孔岩石介质的计算域。其中，相邻物质点间距为Δx=0.000 5 m，邻域半径δ=3.015Δx，材料的弹性模量为E=22GPa，密度为ρ=2 462kg·m-3，泊松比ν=0.25。在此次模拟中可不用设置虚拟物质点作为固体层的边界区域，且各物质点的渗透率各个方向均相同。在流体层设置10 000 个网格表征流体的计算域，将固体层中的物质点与流体层中的网格节点一一对应，另流体密度ρ=1 000kg·m-3，流体黏度系数μ=1.0×10-3Pa·s。因为假设各物质点的渗透率各向同性，故各网格节点的渗透率也各向相同。

用解析解对数值结果进行验证。Sneddon 和Lowengrub［30］给出了在y=0 平面上长度为2lc的单个预制裂隙，在平面应变假设下，y方向压力驱动的位移为

图10 数值模拟结果与解析解沿裂纹的垂直位移对比Fig.10 Comparison of vertical displacement along crack of numerical simulation results and analytical solution

3.2.2 预制倾斜裂隙试样

图11 含有预制水平裂隙试样在不同时间步下裂纹扩展与水压分布Fig.11 Crack propagation and water pressure distribution of specimen containing a pre-existing horizontal flaw at different time steps

为了验证近场动力学与有限体积法耦合方法能对基质岩石中的渗流问题进行有效的数值模拟，对含有预制倾斜裂隙的岩石基质中的渗流过程进行了模拟验证。二维含预制倾斜裂隙的岩石介质模型的几何尺寸如图12 所示，将试件简化为一定厚度的矩形，其几何尺寸L、W分别为0.15m、0.15m。在该岩石介质模型的中间位置预制一个角度为θ=30°的倾斜裂隙，该预制裂隙长度为0.015m。将岩石介质模型的4 个边界均设为排水边界，即P0=0MPa，并且在预制倾斜裂隙内部施加水压p=1×104Pa·s-1。另外多孔岩石介质的渗透率为kr=1.0×10-12m2，而裂隙中的等效渗透率为kf=1.333×10-6m2。

图12 含有预制倾斜裂隙模型的几何尺寸Fig.12 Geometric dimensions of model with a preexisting inclined flaw

在固体近场动力学数值模型中，采用10 000 个近场动力学物质点表征相应的多孔岩石介质的计算域。其中，相邻物质点的间距为Δx=0.001 5 m，邻域半径δ=3.015Δx，材料的弹性模量为E=15GPa，密度为ρ=2 685kg·m-3，泊松比ν=0.15。在此次模拟中可不用设置虚拟物质点作为固体层的边界区域。在流体层设置10 000个网格表征流体的计算域，将固体层中的物质点与流体层中的网格节点一一对应，另流体密度为ρ=1 000kg·m-3，流体黏度系数为μ=1.0×10-3Pa·s。为了能更好地体现倾斜裂隙的渗透特性，此时需要考虑各物质点的渗透率为各向异性，故各网格节点的渗透率也各向异性。

含有预制倾斜裂隙的岩石介质在注水条件下当施加压力达到9.006MPa 时预制裂纹起裂，则裂纹扩展过程中的裂纹分布和压力分布的变化如图13所示。可以明显从图中看出，岩石介质中的水压力从注水孔中心向四周扩散，且岩石注水孔处的水压高于四周水压。综合来看，与预制水平裂隙模拟结果相似，新生的水力裂纹从预制裂隙尖端起裂，随着水压增大，裂纹沿着预制裂隙方向扩展，且裂纹宽度无明显变化。同时，最大水压力存在于裂纹内部，水压在沿裂纹扩展方向的压力梯度远小于垂直于裂纹扩展方向的压力梯度。可以将裂纹内部的水压力视作均匀分布。通过与已有的水力压裂裂纹扩展试验［30］对比，模拟计算中裂纹的起裂压力9.006MPa与试验中试件起裂压力8.908MPa 的结果基本接近，且模拟的裂纹扩展形态与文献中的试验结果基本吻合，如图14 所示，说明了近场动力学与有限体积法耦合方法在考虑渗透率各向异性的情况下，能够更加贴合实际地对岩土工程中的水力压裂问题进行高效模拟。

图13 含有预制倾斜裂隙试样在不同时间步下裂纹扩展与水压分布Fig.13 Crack propagation and water pressure distribution of specimen containing a pre-existing inclined flaw at different time steps

图14 实验室预制30°倾斜裂隙水力压裂试验结果［31］Fig.14 Hydraulic fracturing test results of laboratory prefabricated 30°inclined fracture[31]

依据现有参数开展了不同围压对水力压裂裂纹扩展的影响研究。一共设置了5 种工况，具体工况情况如表1所示。

表1 不同工况的围压设置情况Tab.1 Confining pressure setting under different conditions

模拟结果如图15～19 所示，对比5 个工况可以知道，围压的存在会改变裂隙起裂压力，而围压差的存在对水力压裂裂缝扩展具有诱导作用，裂纹总是沿着最大主应力方向延伸，这一现象与试验结果几乎吻合。对比起裂压力，可以看出，当围压差相同的情况下，围压越大，岩石的起裂压力就越大；围压越小，则岩石越容易被裂隙内水压力压裂。

图15 工况1裂纹扩展和水压分布Fig.15 Crack propagation and water pressure distribution under Condition 1

图16 工况2裂纹扩展和水压分布Fig.16 Crack propagation and water pressure distribution under Condition 2

图17 工况3裂纹扩展和水压分布Fig.17 Crack propagation and water pressure distribution under Condition 3

图18 工况4裂纹扩展和水压分布Fig.18 Crack propagation and water pressure distribution under Condition 4

图19 工况5裂纹扩展和水压分布Fig.19 Crack propagation and water pressure distribution under Condition 5

4 结语

针对在近场动力学统一框架下进行流固耦合模拟存在的计算效率问题，开展了普通态型近场动力学与有限体积法耦合模拟分析方法研究，建立了一种用于模拟饱和裂隙孔隙介质中水力裂缝扩展的混合有限体积法和普通态型近场动力学建模分析方法。普通态型近场动力学用来描述固体骨架的变形和捕捉裂纹的发展，而有限体积法则用来以一种计算高效的方式描述基于达西定律的流体渗流问题。

首先，通过多孔介质渗流模拟算例验证了所提方法的有效性，其数值结果与解析解吻合较好，且计算效率大幅提高，能够很好改善近场动力学进行流固耦合模拟存在的计算效率问题。

其次，通过几个数值算例进一步验证了该方法在流体驱动下模拟饱和裂隙孔隙介质中裂纹扩展的能力和主要特点，在流体驱动的水力压裂模拟中观察到的现象与实验结果一致。

最后，通过对5个工况的模拟，展现了围压差对水力压裂裂缝扩展具有诱导作用的现象，而且模拟结果表明裂纹总是沿着最大主应力方向延伸，这一现象与试验结果吻合良好。

作者贡献声明：

李卓徽：程序研发，模型构建、数据分析及论文撰写。

周宗青：研究思路指导，论文审阅与修改。

高成路：研究内容及方案指导。

张道生：提供编程帮助。

白松松：提供理论帮助。