岩体水力压裂应力-渗流耦合近场动力学模拟

摘要

裂隙岩体应力-渗流耦合机制是油气开采、地应力测量与地质灾害防控等岩土工程活动的理论基础。基于近场动力学非局部作用思想提出了物质点双重覆盖理论模型，通过将近场动力学在模拟固体材料变形损伤与地下水渗流两方面的优势相结合，采用“混合”时间积分方案，构建了流体压力驱动条件下裂隙岩体应力-渗流耦合的常规态型近场动力学模拟方法，并将其应用于空心圆柱体注水试验模拟，揭示了水力裂隙起裂、扩展和贯通的作用机制，通过与室内试验及传统数值方法计算结果对比验证了模拟方法的有效性。模拟结果显示，空心圆柱体注水试验过程中岩体的变形和破坏完全是由水力驱动的，水力裂隙的产生是随机的，不需要指定裂隙扩展路径，并且水力压裂过程中致使试件破裂的能量存在积蓄-释放过程，应用近场动力学方法可以较好地捕捉该现象。

关键词

裂隙岩体; 近场动力学; 水力压裂; 应力-渗流耦合; 数值模拟

方法与软件研发。E-mail： [e-mail]

地下水是导致工程岩体失稳破坏的主要原因之一^［1］，地下水渗流与地应力耦合作用下工程岩体的变形、损伤、破坏和失稳机制是众多研究人员关心的课题，尤其是岩体水力压裂的干预和控制成为当前研究热点。岩体水力压裂是指由于水压升高引起岩体中裂隙发生与扩展的物理现象，在提高油气产量方面具有重要的应用价值^［2］，同时，水力压裂法也是国际岩石力学学会建议的地应力测量方法之一。^［3］。在岩土工程领域，岩体水力压裂又作为一种地质灾害广泛存在，如大坝漏水、隧道突涌水等，往往造成严重的伤亡事故。因此，科学描述岩体水力压裂过程、揭示岩体应力-渗流耦合破坏机制，对于油气开采和灾害防控等均具有重要意义。

针对岩体水力压裂机理，国内外学者开展了大量力学分析，基于弹性理论、损伤理论、断裂理论和分形理论等先后提出了多种描述方法^［4-6］，但是大量的现场实践表明，现有理论模型很难对大型水力压裂中观测的多缝起裂和扩展路径进行合理的解释。在室内试验方面，国内外学者针对水力压裂以及裂缝扩展规律开展了大量研究，其中真三轴模拟压裂试验是最为常见的一种形式。陈勉等^［7］采用大尺寸真三轴模拟试验系统结合声发射监测装备对天然岩样和人造岩样进行水力压裂裂缝扩展机理模拟实验，揭示了地应力、断裂韧性、节理和天然裂缝等因素对水力裂缝扩展的影响规律。赵益忠等^［8］利用真三轴模拟压裂实验系统对多种岩样进行了水力压裂裂缝起裂及裂缝扩展模拟实验，得到了压后裂缝几何形态和压裂过程中压力随时间的变化规律。郭印同等^［9］、衡帅等^［10］、李芷等^［11］、侯振坤等^［12］采用真三轴岩土工程模型试验机结合工业CT扫描技术建立了一套页岩水力压裂物理模拟与压裂缝表征方法，并开展了页岩水力压裂物理模拟试验研究，实现了水力压裂裂缝的空间形态描述，并揭示了水力裂缝与天然裂缝的相互作用规律。虽然真三轴压裂测试系统结合工业CT扫描、声发射装置、X衍射等先进技术成为研究复杂裂缝网络形成机理的主要试验手段，但是由于岩石材料内部结构和组成成分的复杂性，想要直观展示水力裂缝的扩展过程还是比较困难的。

理论、试验和模拟是研究岩石力学问题的三大途径。随着数值模拟理论与方法的丰富和完善，研究人员采用有限元、边界元和颗粒离散元等方法针对水力压裂过程开展了大量研究。谢兴华^［13］开发了三维水力劈裂有限元计算程序，研究了开裂对渗透性的影响规律。薛炳等^［14］采用渗流-应力耦合模型对水力压裂过程进行了三维有限元研究。袁志刚等^［15］采用ANSYS软件模拟了穿层钻孔水力压裂裂缝的扩展。彪仿俊等^［16］使用ABAQUS软件针对水力压裂水平裂缝影响参数开展了数值模拟研究。Francisco等^［17］利用扩展有限元模拟了多孔岩石中水力裂缝和天然裂缝的交互作用，揭示了地应力差和接近角等参数在裂缝网络发展中的作用规律。王涛等^［18］利用扩展有限元法进行了页岩水力压裂裂缝扩展模拟。Olson^［19］利用边界元理论建立了水平井多段压裂的裂缝扩展模型。王理想等^［20］基于连续-非连续单元法和中心型有限体积法提出了解决水力压裂流固耦合问题的二维混合数值计算模型。Yang等^［21］利用颗粒离散元软件PFC开展了大量水力裂缝扩展规律研究，揭示了射孔角度和射孔布置对同步压裂裂缝扩展机理和岩石破坏模式的影响规律。

虽然传统数值方法在模拟岩体水力压裂过程中取得了较好的研究结果，但是岩体水力压裂过程是典型的连续-非连续应力-渗流耦合过程，传统连续变形分析方法和非连续变形分析方法将岩体材料视为完全连续体或者完全非连续体，均与实际情况有所出入。基于非局部作用思想的近场动力学不再采用连续变形假设，而是通过积分方程求解物理力学问题，突破了传统方法裂纹扩展路径的限制，成为求解连续-非连续问题的一种有效方法^［22-23］。该方法不仅可以模拟固体材料的损伤演化机制，还能够描述地下水在孔隙介质、裂隙介质中的渗流规律。本文将近场动力学在模拟固体材料变形与地下水渗流两方面的优势相结合构建流体压力驱动条件下裂隙岩体应力-渗流耦合的常规态型近场动力学模拟方法，通过与室内试验及传统数值方法计算结果对比验证模拟方法的有效性，揭示天然岩体水力压裂过程中水力裂隙的扩展规律。

1 近场动力学基本原理

近场动力学通过将计算区域离散为有限数量带有物性信息的、具有一定体积的物质点，并且利用所有物质点空间位置关系的变化和对外界扰动信息的响应来表征材料的变形与破坏。根据近场动力学的非局部作用思想（示意图见图1），计算区域R_x内的任意物质点i的位置矢量${\mathit{x}}_i$与其一定范围内所有物质点的位置矢量${\mathit{x}}_j$（${\mathit{x}}_j\in R_x,‖{\mathit{x}}_i-{\mathit{x}}_j‖\le \delta$）存在相互作用，其中$\delta$为邻域（horizon），为一个正常数。${\mathit{x}}_i$只与其邻域范围内的物质点存在相互作用，而与超出这个局部区域的物质点之间的相互作用为零。也就是说，近场动力学研究的是物质点与其邻域范围内的所有物质点相互作用的物理现象，非局部作用思想使近场动力学在模拟裂隙岩体破坏问题时能够实现跨裂纹计算^［24-25］。

图1  近场动力学非局部作用示意

Fig. 1  Schematic diagram of nonlocal action of peridynamics

在键型近场动力学模型中，物质点对之间存在大小相等、方向相反的力密度矢量。假设${\mathit{x}}_i$与其邻域内任意其他物质点的${\mathit{x}}_j$之间的相互作用力分别为

$${\mathit{T}}_{ij}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)=\frac{\mathrm{1}}{\mathrm{2}}\mathit{f}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)$$

（1）

$${\mathit{T}}_{ji}({\mathit{u}}_i-{\mathit{u}}_j,{\mathit{x}}_i-{\mathit{x}}_j,t)=-\frac{\mathrm{1}}{\mathrm{2}}\mathit{f}({\mathit{u}}_i-{\mathit{u}}_j,{\mathit{x}}_i-{\mathit{x}}_j,t)$$

（2）

式中：T为力密度矢量，N·m^-6；u为物质点位移矢量，m；t为时间，s；f为本构力函数（constitutive function），包含材料的所有本构信息，不再以传统的应力-应变关系形式出现，不再需要位移场的连续性假设及对位移场进行空间求导。因此，无论位移场的连续与否，都不再影响基本方程的求解，近场动力学模型在求解不连续问题时不再出现数值奇异性等病态特征。

近场动力学基本方程^［26］为

$$\rho ({\mathit{x}}_i)\ddot{\mathit{u}}({\mathit{x}}_i,t)=\underset{H_{x_j}}{\int }({\mathit{{\rm T}}}_{ij}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)-{\mathit{T}}_{ji}({\mathit{u}}_i-{\mathit{u}}_j,{\mathit{x}}_i-{\mathit{x}}_j,t))d{V}_{x_j}+\mathit{b}({\mathit{x}}_i,t])$$

（3）

式中：$\rho$为物质点密度，kg·m^-3；$\ddot{\mathit{u}}$为物质点加速度，m·s^-2；$V$为物质点体积；m³；H为物质点的非局部作用范围；b为外载荷密度，N·m^-3。据此可得由本构力函数表示的${\mathit{x}}_i$的运动方程，如式（4）：

$$\rho ({\mathit{x}}_i)\ddot{\mathit{u}}({\mathit{x}}_i,t)=\underset{H_{x_j}}{\int }\mathit{f}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)\mathrm{d}{V}_{x_j}+\mathit{b}({\mathit{x}}_i,t)$$

（4）

在近场动力学中，物质点键的伸长率s可以由键的相对位置$\mathit{\xi }={\mathit{x}}_j-{\mathit{x}}_i$和相对位移$\mathit{\eta }={\mathit{u}}_j-{\mathit{u}}_i$表示，见式（5）：

$$s=\frac{|\mathit{\eta }+\mathit{\xi }|-|\mathit{\xi }|}{|\mathit{\xi }|}$$

（5）

在近场动力学本构关系中通过引入标量函数$\mu (t,\mathit{\xi })$表征材料的破坏特性，如图2所示，即当键的伸长率s超过极限值$s_{\mathrm{0}}$时，则认为对应的键断裂，此时物质点间相互作用力消失。

图2  物质点间相互作用力与伸长率的关系^［27］

Fig. 2  Relationship between interaction force and stretch of material points^[27]

因此，对于微观弹脆性（PMB）材料^［27］，其本构力函数为

$$f(t,\mathit{\xi })=\mu (t,\mathit{\xi })cs$$

（6）

式中：c为微弹性模量，Pa·m^-4；$\mu (t,\mathit{\xi })$为具有历史依赖性的标量值函数，表示物质点键的完整性，即当$\mu (t,\mathit{\xi })=\mathrm{1}$时键是完整的，而当$\mu (t,\mathit{\xi })=\mathrm{0}$时键是断裂的，其表述为

$$\mu (t,\mathit{\xi })=\left\{\begin{array}{cc}\mathrm{1},& s<s_{\mathrm{0}}\\ \mathrm{0},& \mathrm{其}\mathrm{他}\end{array}\right.$$

（7）

式中：$s_{\mathrm{0}}$为临界伸长率，表征物质点键所具有的最大伸长量。需要注意的是，尽量PMB材料在初始状态下是各向同性的，但是变形过程中某些特定方向键的断裂将导致后续响应的各向异性。

据此，可以通过统计物质点断裂键的数量与初始状态的比值定义物质点的局部损伤，如式（8）：

$$\phi (x,t)=\mathrm{1}-\frac{{\int }_{H_x}\mu (x,t,\mathit{\xi })\mathrm{d}{V}_{\mathit{\xi }}}{{\int }_{H_x}\mathrm{d}{V}_{\mathit{\xi }}}$$

（8）

由此可见，物质点的局部损伤是一个从0到1之间变化的值。当$\phi =\mathrm{1}$时，所有与物质点相互作用的键全部断裂；当$\phi =\mathrm{0}$时，所有与物质点相互作用的键全部完整。局部损伤代表了材料内部裂隙的形成过程，近场动力学本构力函数中包含了材料的损伤与断裂描述，在模拟材料的损伤破坏时，不再需要附加的断裂准则。

为了描述材料的损伤程度，将计算区域围岩损伤状态与完整状态之比定义为损伤指数^［28］，如式（9）：

$$\chi =\frac{\int \phi (x,t)\mathrm{d}{V}_x}{\int \mathrm{d}{V}_x}$$

（9）

因此，当损伤程度越大时，$\chi$值越高。

键型近场动力学模型为解决材料损伤破坏等不连续问题提供了有效的模拟手段，但是依然存在3个主要问题：一是过度简化，泊松比受到了限制；二是键型近场动力学以对点力描述材料的本构特征，无法与传统应力-应变本构关系直接建立联系；三是无法区分材料的体积变形和形状变形，因此无法描述材料的塑性不可压缩性。因此，Silling等^［29］通过引入“状态”的数学概念建立态型近场动力学模型。依据对点力是否沿键长方向，态型近场动力学模型又可以分为常规态型和非常规态型。

在常规态型近场动力学模型中，物质点间的相互作用力与键长方向一致，其力密度矢量也满足角动量守恒定律，因此，力密度矢量可定义为

$${\mathit{T}}_{ij}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)=\frac{\mathrm{1}}{\mathrm{2}}A\frac{{\mathit{y}}_j-{\mathit{y}}_i}{|{\mathit{y}}_j-{\mathit{y}}_i|}$$

（10）

$${\mathit{T}}_{ji}({\mathit{u}}_i-{\mathit{u}}_j,{\mathit{x}}_i-{\mathit{x}}_j,t)=-\frac{\mathrm{1}}{\mathrm{2}}B\frac{{\mathit{y}}_j-{\mathit{y}}_i}{|{\mathit{y}}_j-{\mathit{y}}_i|}$$

（11）

式中：y为变形后物质点的位置矢量；A、B为近场动力学参数，可以由应变能密度变分得到，如式（12）、式（13）：

$$A=\mathrm{4}{\omega }_{ij}[ad{\theta }_i\frac{{\mathit{y}}_j-{\mathit{y}}_i}{|{\mathit{y}}_j-{\mathit{y}}_i|}\cdot \frac{{\mathit{x}}_j-{\mathit{x}}_i}{|{\mathit{x}}_j-{\mathit{x}}_i|}+b(|{\mathit{y}}_j-{\mathit{y}}_i|-|{\mathit{x}}_j-{\mathit{x}}_i|)]$$

（12）

$$\begin{array}{l}B=\mathrm{4}{\omega }_{ij}[ad{\theta }_j\frac{{\mathit{y}}_i-{\mathit{y}}_j}{|{\mathit{y}}_i-{\mathit{y}}_j|}\cdot \frac{{\mathit{x}}_i-{\mathit{x}}_j}{|{\mathit{x}}_i-{\mathit{x}}_j|}+b(|{\mathit{y}}_i-{\mathit{y}}_j|-\\ |{\mathit{x}}_i-{\mathit{x}}_j|)]\end{array}$$

（13）

式中：ω为量纲一化影响函数，控制${\mathit{x}}_i$周围质点对其作用力的大小；θ为体积应变；a、b、d为近场动力学参数，可以由近场动力学应力能密度与经典连续介质力学应变能密度等效求出。依据Madenci和Oterkus^［30］的研究，不同维度条件下近场动力学参数值为

$$\left\{\begin{array}{c}\begin{array}{cc}a=\frac{\mathrm{1}}{\mathrm{2}}(K-\frac{\mathrm{5}G}{\mathrm{3}}),b=\frac{\mathrm{15}G}{\mathrm{2}\mathrm{\pi }{\delta }^{\mathrm{5}}},d=\frac{\mathrm{9}}{\mathrm{4}\mathrm{\pi }{\delta }^{\mathrm{4}}}& \mathrm{三}\mathrm{维}\end{array}\\ \begin{array}{cc}a=\frac{\mathrm{1}}{\mathrm{2}}(K-\mathrm{2}G),b=\frac{\mathrm{6}G}{\mathrm{\pi }h{\delta }^{\mathrm{4}}},d=\frac{\mathrm{2}}{\mathrm{\pi }h{\delta }^{\mathrm{3}}}& \mathrm{二}\mathrm{维}\end{array}\\ \begin{array}{cc}a=\mathrm{0},b=\frac{E}{\mathrm{2}A{\delta }^{\mathrm{3}}},d=\frac{\mathrm{1}}{\mathrm{2}A{\delta }^{\mathrm{2}}}& \mathrm{一}\mathrm{维}\end{array}\end{array}\right.$$

（14）

式中：h为二维结构厚度，m；A为物质点横截面积，m²；K和G分别为体积模量和剪切模量，Pa。

$$\left\{\begin{array}{ccc}\begin{array}{c}K=\frac{E}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mu )},G=\frac{E}{\mathrm{2}(\mathrm{1}+\mu )}\\ K=\frac{E}{\mathrm{2}(\mathrm{1}-\mu )},G=\frac{E}{\mathrm{2}(\mathrm{1}+\mu )}\end{array}& \begin{array}{l}\begin{array}{c}\mathrm{三}\mathrm{维}\\ \end{array}\\ \mathrm{二}\mathrm{维}\end{array}& \begin{array}{c}\\ \end{array}\end{array}\right.$$

（15）

式中：E为弹性模量，Pa；$\mu$为泊松比；${\omega }_{ij}$利用量纲分析可设定其大小为

$${\omega }_{ij}=\frac{\delta }{|{\mathit{x}}_j-{\mathit{x}}_i|}$$

（16）

同样，利用近场动力学力密度表示柯西应力，并得到物质点${\mathit{x}}_i$体积应变的一般形式为

$${\theta }_i=d{\displaystyle \sum _{j=\mathrm{1}}^N}{\omega }_{ij}{s}_{ij}\frac{({\mathit{y}}_j-{\mathit{y}}_i)\cdot ({\mathit{x}}_j-{\mathit{x}}_i)}{|{\mathit{y}}_j-{\mathit{y}}_i|}{V}_j$$

（17）

近场动力学不仅可以描述固体材料的损伤破坏，也可以用来描述地下水渗流。在近场动力学模型中，计算区域中任意一个物质点${\mathit{x}}_i$与其邻域内其他所有物质点${\mathit{x}}_j$通过键的形式相互作用，通过将物质点看作是储存水的容器，而将键看作是传递水的通道，可以模拟地下水在多孔介质中的渗流过程^［31］，其中键两端物质点存在的水力势能差被认为是驱动水运输的动力。在键型近场动力学理论中，认为每个键上的流体运输是相互独立的，与其他键不产生相互影响。对于二维均匀各向同性多孔介质的单向流动过程，其近场动力学基本方程为^［32］

$$\begin{array}{l}\frac{\partial }{\partial t}({\rho }_{\mathrm{w}}({\mathit{x}}_i)\phi ({\mathit{x}}_i))=\\ \underset{H_x}{\int }\underline{\omega }({\mathit{x}}_i,{\mathit{x}}_j)\frac{{\rho }_{\mathrm{w}}({\mathit{x}}_i,{\mathit{x}}_j)}{{\mu }_{\mathrm{w}}}·\frac{\mathrm{2}k}{\mathrm{\pi }\delta }\frac{\Phi ({\mathit{x}}_j)-\Phi ({\mathit{x}}_i)}{{‖{\xi }_{x_i{x}_j}‖}^{\mathrm{2}}}\mathrm{d}{A}_{x_j}+R({\mathit{x}}_i)\end{array}$$

（18）

式中：$\Phi$为物质点的有效水力压强，Pa；${\mu }_{\mathrm{w}}$为流体黏度，Pa$·$s；$k$为材料渗透率，m²；${\rho }_{\mathrm{w}}$为流体密度，kg·m^-3；$\phi$为孔隙率；R为源汇项，kg·s^-1；$\underline{\omega }$为量纲一化影响函数，代表某个物质点的${\mathit{x}}_i$的渗流状态受其邻域内其他物质点的${\mathit{x}}_j$影响程度的大小，取$\underline{\omega }({\mathit{x}}_i,{\mathit{x}}_j)=\mathrm{1}/‖{\xi }_{x_i{x}_j}‖$。在进行二维均匀各向同性孔隙介质稳态渗流模拟时，地下水渗流场水力压强可以表示为

$$\Phi ({\mathit{x}}_i)=P({\mathit{x}}_i)+g{\rho }_{\mathrm{w}}({\mathit{x}}_i)z({\mathit{x}}_i)$$

（19）

式中：P为流体压强，Pa；g为重力加速度，m·s^-2；$z$为垂直方向水头高度，m。对于恒温条件下微小恒定压缩流，流体密度可以表示为

$${\rho }_{\mathrm{w}}({\mathit{x}}_i)={\rho }_{\mathrm{0}}(\mathrm{1}+\gamma (P({\mathit{x}}_i)-P_{\mathrm{0}}))$$

（20）

式中：${\rho }_{\mathrm{0}}$为流体密度在压力为$P_{\mathrm{0}}$时的密度，kg·m^-3；$\gamma$为流体的压缩系数，是描述流体压缩性大小的物理量，被定义为单位压力变化时引起的液体单位体积的变化量，水体压缩系数与压力和温度有关，Pa^-1。

当岩体渗流主要受裂隙控制时，一般采用离散裂隙网格模型来模拟岩体的渗流特征。该模型假设完整岩石的渗透性远远低于节理裂隙，因此岩石的渗透性可以忽略，认为流体只在网络状的裂隙中运动，因此该模型比等效连续介质模型更加符合实际情况。

由于裂隙是岩体渗流的主要通道，而岩体在地质构造的长期影响下，裂隙的形状差异十分严重，并且裂隙具有明显的方向性，平行于裂隙方向的渗透系数显然比垂直于裂隙方向的渗透系数大得多，因此裂隙岩体渗流具有明显的各向异性特征。影响裂隙岩体渗流特性的主要因素有：裂隙的分组、空间分布、密度和规模等；裂隙组及单条裂隙的产状、大小、粗糙度、连通性和充填特性等；裂隙之间的空间交互关系；岩块孔隙率和渗透性；水的流动状态等^［33］。由于这些因素的复杂性、随机性和不确定性，在数值计算和工程分析中通常难以建立准确的表征裂隙渗流特性的关系式，因此，通常将裂隙简单概化为由具有一定长度、倾向、倾角和宽度的空间离散裂隙网格。立方定律是岩体水力学的理论基础，将其写成经典达西定律的形式为

$$q=k_{\mathrm{f}}\frac{g{w}_{\mathrm{d}}J}{\nu }$$

（21）

式中：w_d为裂隙的宽度，m；$\nu$为流体运动黏滞系数，m²·s^-1；J为水力梯度；$k_{\mathrm{f}}$为缝隙的渗透率，m²，有

$$k_{\mathrm{f}}=\frac{{w_{\mathrm{d}}}^{\mathrm{2}}}{\mathrm{12}}$$

（22）

由质量守恒定律，根据式（18）同理可推导出当$\underline{\omega }({\mathit{x}}_i,{\mathit{x}}_j)=\mathrm{1}/‖{\mathit{\xi }}_{{\mathit{x}}_i{\mathit{x}}_j}‖$时，二维裂隙渗流的近场动力学基本方程为

$$\frac{\partial }{\partial t}({\rho }_{\mathrm{f}}({\mathit{x}}_i){\phi }_{\mathrm{f}}({\mathit{x}}_i))=\underset{H_x}{\int }\underline{\omega }\frac{{\rho }_{\mathrm{f}}({\mathit{x}}_i,{\mathit{x}}_j)}{{\mu }_{\mathrm{f}}}\frac{\mathrm{2}{k}_{\mathrm{f}}}{\mathrm{\pi }\delta }\frac{{\Phi }_{\mathrm{f}}({\mathit{x}}_j)-{\Phi }_f({\mathit{x}}_i)}{{‖{\xi }_{{\mathit{x}}_i{\mathit{x}}_j}‖}^{\mathrm{2}}}\mathrm{d}{A}_{x_j}+R_{\mathrm{f}}({\mathit{x}}_i)$$

（23）

式中：${\rho }_{\mathrm{f}}$为裂隙流体密度，即穿越裂隙的流体键两端物质点的流体密度，kg·m^-3，一般可取为平均值；${\phi }_{\mathrm{f}}$为裂隙空间体积与物质点体积的比值；${\Phi }_{\mathrm{f}}$为裂隙表面物质点的有效水力压强，Pa；$R_{\mathrm{f}}$为裂隙源汇项，kg·s^-1；${\mu }_{\mathrm{f}}$为裂隙流体黏度，Pa$·$s。

2 物质点双重覆盖理论模型

天然岩体中赋存有大量的孔隙和裂隙，不仅导致岩体的力学性质发生改变，也影响着其渗透特征。受地质构造的影响，地下工程围岩不仅处于一定的地应力环境中，同时还处于一定的渗流场中，应力场与渗流场之间存在相互影响，这种相互影响被称为应力-渗流耦合作用（或流-固耦合作用）。裂隙岩体流固耦合作用主要表现为2个方面：一是地下水在裂隙岩体中渗流时导致岩体中产生渗流作用力，主要包括静水压力和动力水压力，从而改变岩体原始应力状态；二是裂隙岩体应力状态的改变导致岩体结构发生变形甚至破坏，进而改变岩体的渗透性能及地下水渗流路径^［34］。

在经典土力学中，Terzaghi（太沙基）提出了著名的饱和土有效应力原理，该原理将饱和土体看作由土颗粒和孔隙水2种介质组成，土体的变形由土颗粒骨架所承担的有效荷载控制^［35］。Terzaghi有效应力原理在分析饱和土体的力学特性时得到了符合客观规律的认识，对促进土力学成为一门独立的学科起到了重要的推动作用。然而由于岩体材料的固相是连通的而非散体材料，因此Terzaghi有效应力原理不再适用。为此，通过引入考虑饱和材料体积压缩变形的修正系数$\alpha$可以建立适用于岩体材料的有效应力表达式^［36］，如式（24）：

$$\sigma \mathrm{\text{'}}=\sigma -\alpha P$$

（24）

式中：σ为应力，Pa；$\alpha$为Biot系数，代表了孔隙水压力对有效应力的贡献程度。程远方等^［37］采用Cross-plotting法、排水实验法与声波动态法测试了不同渗透性天然岩芯的Biot系数值，结果表明：高渗岩芯为0.90~1.00，中渗岩芯为0.80~0.90，而低渗岩芯为0.60~0.75。

利用有效应力原理，本文提出的模拟裂隙岩体应力-渗流耦合过程的近场动力学物质点双重覆盖理论模型如图3所示。即物质点不仅具有表征材料变形的固体力学性质，同时也是储水“容器”，地下水通过物质点间的键传递水压力。由于孔隙水压力的作用，物质点固体键发生了弱化作用，同时还需要承受裂隙表面的渗透压力。与固体键的状态完全相反，岩体损伤导致固体键的断裂正好激活了流体键，此时，断裂的键具有更高的渗流能力。

图3  应力-渗流耦合模拟的物质点双重覆盖理论模型^［28］

Fig. 3  Theoretical model of material point double repeated cap for stress-seepage coupling^[28]

考虑地下水渗流时，近场动力学物质点间的相互作用力（式（10）、式（11））可以改进为

$${\mathit{T}}_{ij}({\mathit{u}}_j-{\mathit{u}}_i,{\mathit{x}}_j-{\mathit{x}}_i,t)=(\frac{\mathrm{1}}{\mathrm{2}}A-\frac{\mathrm{2}\alpha P_i}{mh}{\omega }_{ij}|{\mathit{x}}_j-{\mathit{x}}_i|)\frac{{\mathit{y}}_j-{\mathit{y}}_i}{|{\mathit{y}}_j-{\mathit{y}}_i|}$$

（25）

$${\mathit{T}}_{ji}({\mathit{u}}_i-{\mathit{u}}_j,{\mathit{x}}_i-{\mathit{x}}_j,t)=(\frac{\mathrm{1}}{\mathrm{2}}B-\frac{\mathrm{2}\alpha P_j}{mh}{\omega }_{ij}|{\mathit{x}}_i-{\mathit{x}}_j|)\frac{{\mathit{y}}_i-{\mathit{y}}_j}{|{\mathit{y}}_i-{\mathit{y}}_j|}$$

（26）

式中：m为体积加权系数，其表达式为

$$m=\underset{H_x}{\int }{\omega }_{ij}({\mathit{x}}_j-{\mathit{x}}_i)({\mathit{x}}_j-{\mathit{x}}_i)\mathrm{d}{V}_{x_j}$$

（27）

对于裂隙表面，物质点间不再发生相互作用，但是由于流体渗透压力的存在，其运动状态依然受到渗透作用力的影响。因此，当物质点间的键断裂后，其力状态可表示为^［38］

$${\mathit{T}}_{ij}({\mathit{x}}_i,t)=-\frac{\mathrm{2}{P}_i}{mh}{\omega }_{ij}|{\mathit{x}}_j-{\mathit{x}}_i|\frac{{\mathit{y}}_j-{\mathit{y}}_i}{|{\mathit{y}}_j-{\mathit{y}}_i|}$$

（28）

$${\mathit{T}}_{ji}({\mathit{x}}_j,t)=-\frac{\mathrm{2}{P}_j}{mh}{\omega }_{ij}|{\mathit{x}}_i-{\mathit{x}}_j|\frac{{\mathit{y}}_i-{\mathit{y}}_j}{|{\mathit{y}}_i-{\mathit{y}}_j|}$$

（29）

可以看到，在有效应力原理中，裂隙表面在水压力的作用下相当于牵引边界，从而降低了材料的总应力状态。

3 数值迭代求解策略

近场动力学运动方程和渗流方程的求解均需要借助数值方法来实现，通过在有限数量离散的物质点上进行时间积分，利用当前时间步的位移增量与渗流增量预测下一时间步物质点的位置、变形和孔隙水压力，据此判断物质点键的损伤与材料的破坏。时间积分是近场动力学求解的核心，在计算过程中驱动着模拟的进行和信息的传递，与经典局部方法相比，近场动力学时间积分方案并无显著差异。

时间积分分为显式时间积分和隐式时间积分。其中，显式时间积分是通过大量的微小时间步长来确定系统对模拟过程的响应，通常用于瞬态模拟，适合于求解材料变形和破坏过程中的动态问题。隐式时间积分通常用于模拟具有相对较大的载荷步并且动态影响可以忽略不计的情况，通过将加速度项设为零来得到准静态模拟的控制方程。

当忽略裂隙岩体破坏过程中的动力现象时，可以将渗透压力作为驱动载荷，逐级施加到模型中，对平衡方程进行迭代求解，从而得到系统在当前荷载条件下的稳定状态。在隐式时间积分方案中，系统的平衡状态可以通过残差来确定，当残差值降低到指定阈值时，准静态加载步被认为是收敛的。因此，本文将近场动力学流-固耦合求解过程划分为固体模块和流体模块分别进行，采用“混合”时间积分方案，即固体模块采用隐式时间积分、流体模块采用显式时间积分，其计算流程如图4所示。

图4  近场动力学流-固耦合模拟计算流程

Fig. 4  Flowchart of peridynamics fluid-solid coupling simulation

4 水力裂隙萌生-扩展过程模拟

岩石注水试验是水力压裂研究最常见的手段之一^［39］。利用水力压裂研究中的经典试验案例对近场动力学流-固耦合方法进行可靠性分析和有效性验证。空心圆柱体注水试验是经典水力压裂试验之一，采用平面应力模型开展二维条件下空心圆环水力压裂过程模拟，如图5所示。试件半径r=25mm，其中间预留一个半径r₀=2mm的注水孔。模型采用环形等距方式离散为7 636个物质点，物质点间距Δ=0.5mm，邻域大小$\delta$=1.5mm。

图5  空心圆柱体水力压裂数值模型

Fig. 5  Numerical model of hydraulic fracturing of hollow cylinder

模型采用均匀、各向同性假设，岩体密度$\rho$=2 090kg·m^-3，弹性模量E=5GPa，泊松比$\mu$=0.3，临界伸长率s₀=0.001。由于岩体水力压裂过程耗时较短，岩体破坏过程中的渗流多以裂隙为主，因此采用裂隙介质模型对材料渗流场进行模拟，从而忽略岩石基质的渗流作用。岩体材料孔隙率$\phi =\mathrm{0.5}$，Biot系数$\alpha$=1.0，裂隙流体密度${\rho }_{\mathrm{f}}$=1 000kg·m^-3，黏度$\mu$=0.899Pa·s，压缩系数$\gamma$=0.048 5。岩体渗透率与裂隙的宽度的平方成正比，即符合立方定律。试验注水速率采用压力增量的形式进行设置，其大小为$\mathrm{\Delta }P$=10kPa·s^-1。通过固体模块与流体模块的交替迭代可以模拟空心圆柱体注水过程中水力裂隙的演化过程。

空心圆柱体水力压裂过程中岩体的损伤、破坏与渗流演化规律如图6所示。在注水过程中，圆孔内壁孔隙水压力逐级递增，圆柱体受到水压力的作用，在周向承受拉应力作用，并引起材料发生拉伸变形。当内孔水压力较低时，试件所受应力未超过其强度，只发生变形，未出现损伤和破裂。随着内孔水压力的升高，大约在P=6.6MPa时，内孔周向拉应力超过试件抗拉强度，裂纹由内壁产生并向外壁扩展，最终形成贯穿裂隙，试件劈裂破坏。

图6  空心圆柱体水力压裂模拟结果

Fig. 6  Simulation results of hydraulic fracturing of hollow cylinder

在该模拟中岩体的变形和破坏完全是由水力驱动的，水力裂隙的产生是随机的，不需要指定裂隙扩展路径。对比岩体损伤和孔隙水压力分布云图可以看出，流体沿着裂隙扩展路径渗流，并且仅在裂隙表面局部范围内进行。由于采用裂隙渗流模型，在裂隙尖端孔隙水压力为零，因此裂隙内孔隙水压力呈梯度分布。裂隙内孔隙水压力的升高，进一步驱动着裂隙由试件内壁向外壁延伸，并最终形成劈裂破坏。

为了更好地分析空心圆柱体水力压裂过程中岩体的损伤破坏程度，建立岩体损伤指数随注水压力的变化关系，如图7所示。可以看出，在低水压力条件下（P<6.5MPa），岩体保持良好的完整性，此时试件所受应力未超过岩体破坏强度，并没有局部损伤的产生。但是当注水压力超过临界值时，岩体迅速损伤破坏。这表明，空心圆柱体在水力压裂过程中，致使试件破裂的能量存在积蓄-释放过程。材料的损伤和破坏消耗掉了储存在岩体中的变形能，并随着注水压力的增大，进一步产生能量的累积。但是由于裂隙的存在，在裂隙尖端出现应力集中，导致材料更容易发生破坏。

图7  空心圆柱体水力压裂模拟过程中损伤指数变化规律

Fig. 7  Variation of damage index in hydraulic fracturing simulation of hollow cylinder

图8给出了采用超高压大流量渗流-应力耦合试验仪开展的空心圆柱体水力压裂试验结果^［39］以及应用岩石破裂失稳的渗流应力耦合分析系统（F-RFPA）开展的数值模拟结果^［40］。试验和数值模拟结果表明，注水过程中水力裂隙首先出现在内壁表面，并随着水压力的增大，裂隙由内壁向外壁扩展，并最终形成贯通裂隙，试件整体发生劈裂破坏。本文利用近场动力学流-固耦合模拟方法得到了相似的研究结论。但是，由于天然岩体材料的非均匀性，试验过程中裂隙首先出现在内壁面材料强度较弱的某点上，并且由于注水压力引起的应力集中导致材料局部产生高拉应力区，引起试件开裂。F-RFPA采用Weibull分布将模型的力学参数进行随机赋值，以此来模拟岩体材料的非均匀性，得到了单条裂隙贯穿试件表面，而本文建立的近场动力学模型对称性较好，且采用了均匀性、各向同性假设，因此裂隙分布也呈现较好的对称性。

图8  空心圆柱体水力压裂试验和RFPA模拟结果

Fig. 8  Results of hydraulic fracturing test and RFPA simulation for hollow cylinder

5 结语

基于有效应力原理建立近场动力学流-固耦合模拟的物质点双重覆盖理论模型，并开发了基于固体和流体时间步循环迭代的求解器，通过将水力压裂动态问题转化为准静态问题实现水压力驱动作用下裂隙岩体水力压裂过程的模拟。通过开展裂隙岩体水力压裂试验模拟研究，并与前人研究结果进行对比，验证该理论模型的准确性，并取得了如下结论：

（1）基于有效应力原理，通过将近场动力学在模拟固体变形与流体渗流两方面的优势结合起来可以构建表达裂隙岩体流-固耦合效应的物质点双重覆盖理论模型，即物质点不仅储存了关于固体介质变形和损伤的力学信息，同时也是储存流体介质的容器，键不仅反映物质间相互作用力，同时也是流体渗流的通道。

（2）通过开展空心圆柱体水力压裂试验模拟，并与室内试验及前人研究结果进行对比，验证了本文所建立的裂隙岩体应力-渗流耦合模拟方法的准确性和可靠性。

（3）研究结果显示，空心圆柱体水力压裂过程中裂纹首先由内壁产生并向外壁扩展，最终形成贯穿裂隙，试件劈裂破坏，该过程完全是由水力驱动的，水力裂隙的产生是随机的，不需要指定裂隙扩展路径，并且水力压裂过程中致使试件破裂的能量存在积蓄-释放过程，近场动力学流-固耦合模型可以有效模拟水力裂隙萌生扩展的客观规律，为油气开采、地下水灾害防控等提供了分析手段。

作者贡献声明

高成路：程序研发及论文撰写。

李术才：研究思路指导。

周宗青：研究内容制定。

李利平：研究方案指导。

张道生：数值建模。

厉明浩：算例验证。

刘冠男：论文修订。