李双 夏晓舟 裴磊

摘要：推导了含裂隙饱和多孔介质的水压致裂全耦合问题的扩展有限元增量迭代求解格式。通过对含中心裂纹无限大平板的裂尖位移场解析解的分析，按I型和Ⅱ型变形特征进行提取，构建了一种新型的内聚裂尖加强函数，使传统扩展有限元中受裂尖加强影响的节点由每个节点8个裂尖加强自由度缩减到2个（平面问题），大大减小了计算规模。以含中心斜裂纹的多孔介质平板为例，验证了内聚裂尖加强模式的正确性，并对两种裂尖加强模式的雅克比矩阵条件数进行了计算。结果表明：内聚裂尖加强模式下扩展有限元的计算精度与传统扩展有限元的计算精度同阶，且网格较密时，本文处理方法的计算精度高于传统模式的。最后将改进的内聚裂尖加强模式应用到混凝土重力坝的水压致裂模拟中，取得了良好的计算效果。

关键词：扩展有限元;饱和多孔介质;渗流一应力耦合;水压致裂;裂尖加强函数

中图分类号：TV139.14

文献标志码：A

doi：10.3969/j .issn. 1000- 1379.2019.05.023

多孔介质[1-2]是由不同大小和形状的固体颗粒骨架构成的，并至少有一种流体介质填充在介质的孔隙中。全部或部分饱和多孔介质的力学行为一般由固体骨架与孔隙流体的相互作用决定，因此如果考虑在恒定荷载作用下的土壤固结，通常会遇到流体流动和固体骨架变形的耦合问题。水压致裂是油藏最常用的增产措施[2-3]。在传统有限元法中，裂纹路径局限于单元间边界，通常遇到网格依赖性问题。扩展有限元法（ Extended Finite Element Method，簡称XFEM） [4-5]是一种强大而精确模拟不连续介质的方法，该方法使整个裂纹的几何形状独立于网格，并且完全避免了裂纹扩展时网格的重新划分。对于裂纹扩展、孔洞、复合材料等含不连续问题的介质，可以用扩展有限元法有效模拟。XFEM通过引入具有特殊性质的加强函数及与之相应的加强自由度，模拟含间断问题。

本文将多孔介质耦合控制方程与XFEM结合，完成对含裂隙饱和多孔介质的流固耦合数值模拟。同时对裂尖位移场加强模式进行改进，依据线弹性断裂力学裂尖位移场的解析解形式，按I型和Ⅱ型裂纹张开变形特征进行提取，构建了一种新型的内聚裂尖加强函数，模拟裂尖的奇异性[6]，同时将原来裂尖加强模式中每个裂尖影响节点的8个裂尖加强自由度减到2个。通过算例分析证明，内聚裂尖加强函数不仅减小了矩阵规模，而且能更好地模拟含裂隙多孔介质的流固耦合问题，并运用内聚裂尖加强模型分析了混凝土大坝的流固耦合问题。

1 含裂隙多孔介质扩展有限元建模

1.1 单元位移及压力模式

XFEM基于单元分解的基本思想，在传统有限元法的近似解中引入合适的加强函数项来反映研究对象的某些特征，如裂纹或材料界面等的不连续性，同时使研究对象独立于网格划分。

1.1.1 XFEM单元位移模式

为了描述裂隙多孔介质中裂隙处的位移跳跃，假定位移场在裂纹边缘上是强不连续的，并用Heaviside函数和裂尖渐进加强函数，分别描述裂纹面的不连续性及裂尖的奇异性，得到单元内任意一点的位移：

1.1.2 XFEM单元压力模式

为了表示流体流量在垂直于裂缝方向上的跳跃，假设流体压力场是弱不连续的，即压力场连续，而其垂直于裂缝的梯度在裂缝边缘是不连续的，并用修正的水平集函数描述压力梯度的不连续性，单元中压力表达式为[1]

1.2 裂尖加强函数的改进

由式（2）中的基函数组成的裂尖渐进函数F1（x），每个受裂尖影响的节点都有8个裂尖附加自由度[7]，增大了整体刚度矩阵的组装及方程求解的计算量。为了减少裂尖附加自由度个数，本文以断裂力学裂尖位移场的解析解为基函数，按I型和Ⅱ型裂纹张开的变形特征进行提取，在局部坐标系下重新构造能反映裂尖奇异性的内聚裂尖位移加强模式。整体坐标系、裂尖局部坐标及极坐标系示意见图1。

1.3 含裂隙多孔介质的扩展有限元离散

1.3.1 多孔介质的控制方程

变形多孔介质的控制方程包括流固混合相及各相线性动量平衡方程和连续性方程，运用Biot理论得到[9]。

忽略流体相的相对加速度，流固混合相的线性动量平衡方程可表示为

1.3.3控制方程的XFEM空间离散

用XFEM的位移场及压力场模式对弱形式的控制方程（23）进行空间离散，根据Bubnov Galerkin方法，可以得到积分方程即式（23）的离散形式：

2 含裂隙多孔介质扩展有限元程序计算

笔者依据含裂隙多孔介质的流固全耦合方程，编写了裂尖加强模式改进前后的扩展有限元程序，计算了两种裂尖加强模式下雅克比矩阵的条件数（2范数）。

2.1 程序可靠性验证

为验证编写的多孔介质流固全耦合扩展有限元程序的可靠性，将程序的计算结果与已有文献计算结果进行对比。计算模型（见图3）为10 mxlo m的正方形断裂域，中心含一长度为2m的裂纹，为平面应变模型。模型的位移边界条件是左右两边约束U=0，底部约束U_v=0，顶部不约束;压力场边界条件是左右两边不排水，顶部为排水边界，压力设置为0，底部为流量人口边界，q= lxl0^-4m/s。弹性模量E=9 GPa，泊松比为0.4，渗透率k=lxl0^-5m/s，孔隙度为0.3，固相压缩模量为lx10²⁷Pa，流体压缩模量为2.lx10⁹Pa，固相密度为2 000 kg/m3，流体相密度为1 000 kg/m³，流体黏度μU=lxl0^-3Pa.s。时间步长为0.01 s，计算总时间T= 10 s。

T=8 s时含30°斜裂纹的多孔介质竖向位移云图如图4所示。可以看出，裂尖加强模式改进前后的扩展有限元程序所得云图分布规律与文献[1]结果吻合，证明了程序的正确性。两种裂尖加强模式下的压力云图如图5所示。由位移和压力云图可以看出，两种裂尖加强模式都能很好地模拟含30°斜裂纹的多孔介质流固全耦合问题。

2.2 雅克比矩阵条件数对比

为了比较两种裂尖加强模式的优越性，分别以水平中心裂纹、30°及45°中心斜裂纹为例，对比分析两种加强模式下竖向位移和水平位移云图，见图6-图8。可以看出：对于水平裂纹和30°斜裂纹，两种裂尖加强模式都能很好地模拟多孔介质的耦合问题。但对于45°斜裂纹，原来的加强模式所得结果出现上下裂纹面靠近的趋势，而改进的加强模式所得结果并没有出现，所以改进后模拟效果更好。

还可通过程序在4组不同的网格密度下分别查看裂尖加强模式改进前后雅克比矩阵的条件数。水平中心裂纹、30°及45°中心斜裂纹的雅可比矩阵条件数分别见表1-表3。

从雅克比矩阵的条件数来看，网格密度为11 x11时，改进前的雅可比矩阵条件数大于改进后的条件数，网格密度为2lx21时结果相反。网格密度小于3lx31时，两种裂尖加强模式的条件数并无规律，但当网格密度增大时，裂尖加强模式改进后的条件数比加强模式改进前的条件数小，并且网格越大，两者的差距越明显。

3 混凝土重力坝水压致裂模拟

将内聚裂尖加强模式运用到混凝土重力坝水压致裂的渗流一应力全耦合分析中，并针对不同流体黏度和渗透系数对水压致裂的裂缝影响进行对比，重力坝模型及参数见图9、图10及表4。将混凝土视为多孔介质材料，库水按牛顿流体处理。

3.1 流体黏度对孔隙水压力的影响

不同流体黏度下混凝土大坝裂隙的孔隙水压力分布如图11所示。

可以看出，流体黏度从μ= lx 10^-1Pa.s到μ=1×10^-6Pa.s，裂隙注水口附近的孔隙水壓力在逐渐变小，即流体黏度对注水口的孔隙水压力有影响，水压力随黏度的增大而变大。

3.2 渗透系数对孔隙水压力及张开度的影响

混凝土大坝的渗透能力由裂隙的渗透系数表示，渗透系数是影响混凝土重力坝水压致裂的主要因素之一，可以通过改变介质的渗透系数来控制裂隙的萌生及裂缝扩展过程。不同渗透系数下大坝的孔隙水压力及裂隙张开位移随渗透系数的变化见图12和图13。

可以看出，在注水开始一段时间后，渗透系数为30x 10^-9m/s时的孔隙水压力最小，为lx 10^-9m/s时孔隙水压力最大，可见在这一段时间里孔隙水压力与渗透系数成反比关系[11]。同样，在图13中，裂隙张开位移随渗透系数的增大而减少，与文献[11]的结果相符。渗透系数越大，孔隙水压力越小，对裂隙的张开作用越小，则相对裂隙张开位移越小。

4 结论

多孔介质的渗流应力耦合问题一直是研究的热点，本文对含裂隙饱和多孔介质流固全耦合问题进行了数值建模，并结合扩展有限元法建立了含裂隙多孔介质的有限元模型。为了减少裂尖附加自由度个数，对裂尖位移场加强模式进行改进，以线弹性断裂力学裂尖位移场的解析解为基函数，构建内聚裂尖加强模式，将裂尖加强自由度由8个减到2个。通过含中心裂隙的无限平板算例验证了模型的可靠性，查看了改进前后雅克比矩阵的条件数，比较了两种裂尖加强模式的优越性，同时将改进的裂尖加强模式运用到了含裂隙混凝土大坝中，对大坝的参数进行了简单分析，证明了改进的裂尖加强在运用中的可行性。

（1）由位移云图和压力云图可知，在裂纹倾角较小时，改进后的加强模式和改进前的加强模式都能很好地模拟含裂隙多孔介质的流固耦合问题。但当裂纹倾角较大时，改进后的内聚裂尖加强模式模拟效果更好。

（2）由条件数可知，当网格密度较稀疏时，两种裂尖加强模式的条件数并无规律，但当网格密度增大时，2个附加自由度的加强模式条件数比8个附加自由度的加强模式条件数小，并且网格越大，两者的差距越明显。

（3）运用内聚裂尖加强对混凝土大坝的裂隙进行模拟，得出孔隙水压力受流体黏度及渗透系数的影响，同时随着水压力一起变化的还有裂隙张开位移，孔隙水压力越小，裂隙张开位移越小。

综上，内聚裂尖加强模式下的扩展有限元计算精度与传统扩展有限元计算精度同阶，且网格较密时内聚裂尖加强模式的计算精度高于传统模式的，这种方法不仅缩小了矩阵规模，而且能更好地模拟含裂隙多孔介质的流固全耦合问题，将其应用到混凝土重力坝的水压致裂模拟能得到良好的计算效果。

参考文献：

[1] KHOEI A R.Extended Finite Element Method：‘rheory andApplications[M].[S.l.]：[s.n.]， 2014： 409-460.

[2] 柳占立，庄茁，孟庆国，等，页岩气高效开采的力学问题与挑战[J].力学学报，2017，49（3）：507-516.

[3] 杨克己，韩理安，桩基工程[M].北京：人民交通出版社，1992：2- 10.

[4]

MOHAMMADI S.Extended Finite Element Method： for FractureAnalysis of Structures[J].Xfem Fracture Analysis of Composites，2007，45（22）：5675-5687.

[5]李录贤，王铁军，扩展有限元法（XFEM）及其应用[J].力学进展，2005（1）：5-20.

[6]夏晓舟，赵睿祺，章青，基于扩展有限元方法的含裂纹结构在动力荷载作用下的数值仿真[ C]//南方计算力学联络委员会，江苏省力学学会，第十一届南方计算力学学术会议（ SCCM-II）摘要集，南方计算力学联络委员会，江苏省力学学会，2017：1.

[7] 江守燕，杜成斌，一种XFEM断裂分析的裂尖单元新型改进函数[J].力学学报，2013，45（1）：134-138.

[8] 王荣华，章青，夏晓舟，含裂隙饱和多孔介质流一固耦合的扩展有限元分析[J].岩土力学，2017，38（5）：1489-1496.

[9] 王忠福，刘汉东，周敏，等，黄河大堤丁坝流固耦合分析研究[J].人民黄河，2008，30（8）：16-17.

[10] RETHORE J，DE BORST R，ABELLAN M A.ATwo-Scale Approach for Fluid Flow in Fractured Porous Media[J] .Intemational Joumal for Numerical Methods in Engi-neering， 2007， 71（7）：780- 800.

[11]

MOHAMMADNEJAD T，KHOEI A R.An Extended FiniteElement Method for Hydraulic Fracture Propagation in De-formable Porous Media with the Cohesive Crack Model[J].Finite Elements in Analysis and Design， 2013，73：77-95.