上海交通大学学报

上海交通大学学报, 2026, 60(3): 463-472 doi: 10.16183/j.cnki.jsjtu.2024.180

船舶海洋与建筑工程

土体中水力劈裂大规模数值模拟及应用

王风1,2, 陈铁林,1, 樊容1, 李跃鹏2

1 北京交通大学 城市地下工程教育部重点实验室, 北京 100044

2 国能朔黄铁路发展有限责任公司, 河北 肃宁 062350

Large-Scale Numerical Simulation and Application of Hydraulic Fracturing in Soil

WANG Feng1,2, CHEN Tielin,1, FAN Rong1, LI Yuepeng2

1 Key Laboratory for Urban Underground Engineering of the Ministry of Education, Beijing Jiaotong University, Beijing 100044, China

2 Shuohuang Railway Development Co., Ltd., Suning 062350, Hebei, China

通讯作者: 陈铁林,教授,博士生导师,电话(Tel.):010-51688113;E-mail:tlchen1@bjtu.edu.cn.

责任编辑: 王一凡

收稿日期: 2024-05-17   修回日期: 2024-07-21   接受日期: 2024-09-4  

基金资助: 国家重点研发计划资助项目(2021YFB2300805)

Received: 2024-05-17   Revised: 2024-07-21   Accepted: 2024-09-4  

作者简介 About authors

王风(1981—),博士生,高级工程师,从事隧道及岩土工程研究.

摘要

在岩土工程领域,分析土体水力劈裂规律和机理具有重要意义,可为劈裂注浆脉体的扩散和分布规律研究提供依据,指导实际工程.通过使用块对角预处理、预处理的对称拟最小残量(PSQMR)迭代法及改进的稀疏矩阵矢量乘并行算法,开发土体中水力劈裂三维有限元程序,可实现大规模Biot固结有限元计算.实现了基于中央处理器(CPU)串行计算和图形处理器(GPU)并行计算平台的有限元求解,大幅提高计算规模,在个人计算机上可进行工程尺度的大规模有限元计算.通过数值模拟、试验模拟和理论分析与计算结果的对比,验证了所述计算方法的正确性,为研究土体中水力劈裂提供了有力的数值模拟工具.

关键词： 水力劈裂; 块对角预处理; 预处理的对称拟最小残量迭代法; 稀疏矩阵; 并行计算

Abstract

In the field of geotechnical engineering, it is of great significance to analyze the law and mechanism of soil hydraulic fracturing, which can provide a basis for the study of the diffusion and distribution of splitting grouting veins and guide practical engineering. In this paper, a three-dimensional finite element program is developed by using block diagonal preprocessing, preconditioned symmetric quasi-minimal residual (PSQMR) iterative method, and improved sparse matrix vector multiplication parallel algorithm, which enables large-scale Biot consolidation finite element calculation, realizing the finite element solution based on CPU serial computing platform and GPU parallel computing platform, greatly improving the calculation scale, and accomplishing large scale finite element calculation on personal computer. By comparing the numerical simulation results, experimental simulation results, and theoretical analysis with the calculation results, the correctness of the calculation method is verified, which provides a numerical simulation tool for the study of hydraulic fracturing in soil.

Keywords： hydraulic fracturing; block diagonal preprocessing; preconditioned symmetric quasi-minimal residual (PSQMR) iterative method; sparse matrix; parallel computing

PDF (8967KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

王风, 陈铁林, 樊容, 李跃鹏. 土体中水力劈裂大规模数值模拟及应用[J]. 上海交通大学学报, 2026, 60(3): 463-472 doi:10.16183/j.cnki.jsjtu.2024.180

WANG Feng, CHEN Tielin, FAN Rong, LI Yuepeng. Large-Scale Numerical Simulation and Application of Hydraulic Fracturing in Soil[J]. Journal of Shanghai Jiaotong University, 2026, 60(3): 463-472 doi:10.16183/j.cnki.jsjtu.2024.180

水力劈裂是指在高水头压力作用下,流体注入岩土体中,导致岩土体起裂、扩展、贯通,最终在岩土体中形成裂缝或裂缝网络的过程[1].依据此原理开发的水力劈裂技术已被广泛应用于石油[2]、天然气[3]、地热[4]等能源的开发和利用中.此外,在岩土工程领域中,结合此原理发展而来的劈裂注浆技术也广泛应用于岩土工程加固、堵水与灾害防治[5].目前,水力劈裂研究主要涉及理论研究、试验研究以及数值模拟等方式.理论研究方面,由于水力劈裂的物理机制较为复杂,同时岩土介质的变异性较大,往往需要约定一系列假设条件,从而难以反映真实岩土介质水力劈裂的过程[6].而试验研究,一方面受到室内试验与工程实际的比例尺换算制约,开展水力劈裂的试验研究极为复杂[7];另一方面水力劈裂具有隐蔽性[8],需要通过人工开挖的方式获取其发展规律,造成试验过程的复杂性,试验数据难以获取,通过较少次数的试验得到水力劈裂的试验规律并不真实.因此,数值模拟方法已成为研究水力劈裂规律的重要方法.

目前与水力劈裂相关的数值方法可分为以下几类:有限元法[9-10]、扩展有限元法(XFEM)[11-12]、相场法[13-15]、光滑粒子流体动力学(SPH)法[16]、无网格法[17]、离散元方法[18-19]、连续-非连续耦合方法[20-22]等.由于水力劈裂所涉及的计算量巨大,计算规模的突破已成为水力劈裂数值模拟的一大难点.以离散单元法为例,模拟水力劈裂工程的算例通常需要百万级单元量,而岩土领域的商业颗粒离散元软件在个人计算机上运行的计算单元数通常在十万单元以内[23] ;使用相场法完成水力劈裂模拟往往需要在裂纹扩展处加密网格,由于水力劈裂所涉及的变量较多,在三维模拟时所耗费的计算资源急剧增大,通常只能完成试件尺度的水力劈裂模拟[24].现有的数值模拟技术在计算水力劈裂相关问题时,难以进行工程尺度的大规模数值模拟,极大地限制了数值模拟与工程应用之间的联系.另外在对实际工程的模拟中,不仅要考虑最大的计算规模,还要考虑多步骤分析的耗时.为了提高研究方法的实用性,亟需开发一种可以在个人计算机上进行工程尺度大规模水力劈裂求解的数值模拟方法.

近年来,有研究表明图形处理器(GPU)极其适用于海量矩阵运算,其计算效率相较于中央处理器(CPU)可提高数十倍[23,25].因此,在本研究中,采用CPU串行计算和GPU并行计算的高性能计算技术进行水力劈裂有限元求解,从而有效地提高科学计算的速度.同时,在数值计算方法上,本研究利用广义雅克比(GJ)预处理的(preconditioned)对称拟最小残量(symmetric quasi-minimal residual, SQMR)迭代法,即PSQMR迭代法求解对称非正定线性系统,采用稀疏矩阵存储数据,从而大幅提高了在个人计算机上的计算自由度,实现了工程尺度大规模水力劈裂问题的数值模拟.最后,本研究利用上述计算方法,求解了5个水力劈裂案例,模拟结果证明了该方法的正确性.

1 计算方法简介

水力劈裂是一个复杂的过程,包括4个过程的耦合[9]:① 裂隙的起裂和扩展;② 水在裂隙内的流动;③ 裂隙水与周围土体的力学耦合作用,称为“水-土相互作用”;④ 土体中孔隙水与应力的耦合作用,称为“渗流-应力耦合”.实际工程中水力劈裂的应用十分广泛,受到研究人员的广泛关注,并已经形成一系列的算法体系,总体上包括两类.第1类是采用平面裂缝假设,针对均质岩土体进行解析或半解析的计算;第2类是基于连续或者离散的数值模拟方法进行计算.

本文主要是在文献[9]的基础上,将二维算法推广至三维情况,并开展土体中三维水力劈裂过程的数值模拟,耦合应力场-流体场-劈裂破坏,模拟土体中的水力劈裂的发展过程.采用GPU-CPU混合计算平台在个人计算机上进行大规模水力劈裂问题的有限元求解[26],通过GPU并行计算求解运算量大的矩阵矢量相乘部分,从而加速线性方程组的求解.求解方法上采用GJ预处理的SQMR迭代法和压缩稀疏列(CSC)稀疏矩阵存储法,此部分在CPU上运行.通过GPU-CPU混合计算大幅提高在个人计算机上的计算自由度.具体而言:通过Biot固结理论计算得出劈裂过程中土体中任一点的孔隙流体压力和位移;采用莫尔-库伦准则判断土体中裂缝的启裂和扩展;裂缝单元的渗透性采用立方定律.更多的算法细节可参考文献[9],此处限于篇幅不作展开.

Biot固结方程离散后,采用增量形式可表达为

$\left[\begin{array}{ll}K& L\\ {L}^{T}& G\end{array}\right]\left[\begin{array}{l}\Delta u\\ \Delta p\end{array}\right]$=$\left[\begin{array}{l}\Delta f\\ \Delta {p}_{t}\end{array}\right]$

式中:K为对称正定的土体的刚度矩阵;G为对称半正定的孔隙水刚度矩阵;L为连接固-液耦合矩阵;Δu为位移增量;Δp为孔隙水压力增量;Δf为节点力增量;Δpt为当前节点孔隙水压.

流固耦合及水力劈裂过程的有限元计算中,通过空间离散和时间差分,Biot固结方程产生一个对称非正定线性大型稀疏方程组,简化为

Ax=b

式中:A∈Rn×n,是对称非正定的刚度矩阵;x是位移u和孔隙水压力p组成的向量;b∈Rn是荷载向量.此方程的求解是有限元计算的关键,计算时间占比较大,因此该方程组的高效求解一直是有限元计算优化方面的研究热点.

1.1 块对角预处理

由于Biot固结线性方程组的系数矩阵中孔隙水压力刚度矩阵部分的对角元素和土体刚度矩阵差异很大,求解难以收敛.Phoon等[27]提出的GJ法可以在一定程度上解决计算难以收敛的问题,从而优化求解过程.GJ预处理因子可以表示为

MGJ=$\left[\begin{array}{ll}diag K& 0\\ 0& \hat{aS}\end{array}\right]$

式中:a是一个实数,Phoon等[27]发现当其小于1时具有较好的收敛性;$\hat{S}$为Schur补矩阵,

$\hat{S}$=diag(G-LTdiag K-1L)

1.2 PSQMR迭代法及并行计算原理

目前,共轭梯度法是求解对称正定线性方程组最为有效的方法,但是由于Biot固结所离散出的线性方程组为对称非正定的线性方程组,所以采用共轭梯度法受到限制,一般比较常用的是最小残量(MINRES)方法和对称LQ(SYMMLQ)方法.但最小残量法和对称LQ法所使用的预处理子只能是对称正定的.Freund[28]发展了SQMR方法来求解对称非正定线性系统,且该迭代法可以结合任意对称预处理子.

PSQMR迭代法可采用如下方式描述,令M=MT=MLMR为任一可逆的n阶对称矩阵,ML和MR分别为左右侧预处理因子.

选择一个初始值,x0∈Rn,令s0=b-Ax0,t=${M}_{L}^{-1}$s0,q0=${M}_{R}^{-1}$t,r0=‖t‖2,v0=0,p0=sT0q0,d0=0.

从1到最大迭代次数进行迭代计算:

① t=Aqk-1;

② ck-1=qTk-1t;ak-1=$\frac{{p}_{k-1}}{{c}_{k-1}}$;sk=sk-1-ak-1t;

③ t=${M}_{L}^{-1}$sk;vk=$\frac{{‖t‖}_{2}}{{r}_{k-1}}$;ck=$\frac{1}{\sqrt{1+{v}_{k}^{2}}}$;

④ rk=rk-1vkck;dk=${c}_{k}^{2}{v}_{k-1}^{2}$dk-1+${c}_{k}^{2}{a}_{k-1}^{2}$qk-1;

⑤ xk=xk-1+dk;

⑥ 检查计算结果是否收敛,若收敛则停止计算;

⑦ 计算uk=${M}_{R}^{-1}$t;ρk=sTkuk;βk=$\frac{{\rho }_{k}}{{\rho }_{k-1}}$;qk=uk+βkqk-1

其中:k为迭代次数;xk为第k步解向量,x0为任意初始解;sk为第k步残差向量,s0为初始残差向量;t为左侧预处理的中间向量;qk和uk为预处理后的方向向量,q0为经过左右预处理后的初始搜索方向向量;rk为残差范数;ak,ck,vk,βk为标量步长或比例因子;pk和ρk为内积,p0为初始内积;dk为解向量迭代增量.

在PSQMR迭代法的求解中,矩阵矢量乘积、矢量点积和矢量更新是计算的主要部分,其中稀疏矩阵矢量乘是迭代计算中计算量最大的核心部分,且耗时最多,因此对矩阵矢量的优化是整个有限元求解关键所在.

在上述的计算过程中第①步为矩阵矢量相乘,运算量最大,将其交给GPU并行计算处理.在GPU计算过程中,矩阵-矢量乘的计算过程中常产生不必要的数据传输.每一步都是矢量qk-1与矩阵A相乘得到的一个矢量t,A保持不变,因此将矩阵A一次存放到GPU显存中,避免了矢量的往返更替.

1.3 稀疏矩阵矢量乘并行算法的改进(SpMV-GPU)

串行的稀疏矩阵向量乘法SpMV计算过程中,A下三角矩阵中的非零元素与矢量的点积过程以行为单位逐次进行,每一次计算都是独立完成的.当采用GPU并行计算时,每一次A的多行非零元素同时进行计算,得到对应的多个点积结果,如图1所示.图中:tP为线程块(thread)划分数目;U、LA分别为A矩阵的上三角部分和下三角部分;向量b为A和向量x的乘积;wi为向量b的第i个元素.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   GPU并行计算过程

Fig.1   GPU parallel computing process


上三角采用CSC存储,A矩阵与x矢量相乘再相加才能得到最终的b向量,此种方式彼此相关联而不便于实现并行计算.为此,将三角的CSC存储格式转化成压缩稀疏行(CSR)的存储格式,这样就可以保证上三角和下三角都可进行并行计算.采用Fortran语言在PGI Visual Fortran 编译器上编写三维有限元程序,矩阵矢量乘在编程算法上可以采用内积法或外积法两种计算方式.关于编程方法可参照文献[9,26],内容较多不展开叙述.

2 算例

算例1和土体固结理论解比较,验证流固耦合计算的准确性;算例2和水力劈裂理论解比较,验证求解劈裂的准确性;算例3和水力劈裂试验比较,验证三维情况下模拟劈裂形状的准确性;算例4验证三维多步骤模拟中劈裂顺序的影响的合理性;算例5验证工程尺度条件下三维多步骤模拟的能力.

2.1 算例1

饱和土体3 m×3 m×3 m,顶部施加300 kPa的水压力、位移自由,其他边界不透水、法向位移固定.

该问题的控制方程为

Cv$\frac{{\partial }^{2}p}{\partial {z}^{2}}$=$\frac{\partial p}{\partial t}$

边界条件如下:

z=0, p=pc; z=H, $\frac{\partial p}{\partial z}$=0 (t>0)

式中:Cv为固结系数,Cv=$\frac{{k}_{c}}{{m}_{v}{\gamma }_{w}}$,kc为渗透系数,γw为水的容重,mv=$\frac{(1+\nu )(1-2\nu )}{E(1-\nu )}$为土的体积压缩系数, E为弹性模量,ν为泊松比;z为深度;t表示施工时间;pc为边界孔压;H为土层的总厚度或最大排水路径.

初始条件如下:

t=0, p=0; t=∞, p=pc
(0≤z≤H)

控制方程的解为

p=-$\frac{4{p}_{c}}{\pi }\sum _{m=1}^{\infty }\frac{1}{m}$sin$\frac{m\pi z}{2H}$exp(-m2π2Tv/4)+pc
(m=1, 3, 5, …)

式中:Tv=$\frac{{C}_{v}t}{{H}^{2}}$,表示时间因数.

位移的计算公式如下:

ds=mv$\frac{\partial p}{\partial t}$dzdt

式中:ds为微小位移增量.

不同固结时间下,孔压-深度曲线和位移-深度曲线如图2及图3所示.数值模拟计算耗时随自由度的变化与常规求解方法对比如图4所示,可以发现本文所计算程序是常规求解方法1~1.5倍,且两种方法提随着自由度的增加,并行效率都会降低,特别是自由度大于300万以后并行效率下降较快.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   不同固结时间下孔压-深度曲线

Fig.2   Pore pressure versus depth under different consolidation durations


图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   不同固结时间下位移-深度曲线

Fig.3   Displacement versus depth under different consolidation durations


图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   自由度和计算耗时

Fig.4   Degree of freedom and computational time


2.2 算例2

平面应变水力劈裂理论(KGD)模型[29]如图5所示.图中:r为从裂缝尖端到研究点的径向距离;w为裂缝宽度;x为沿裂缝长度方向的坐标;h为模型的高度.KGD模型适用于裂缝长高比较小的情况,其假设条件是: ① 缝高恒定;② 裂缝在水平面上处于平面应变条件;③ 裂缝顶端是尖的;④ 压裂液为牛顿流体;⑤ 裂缝断面为椭圆形,最大缝宽在裂缝中部;⑥ 岩石遵循弹性力学理论;⑦ 不考虑压裂液滤失.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   KGD模型

Fig.5   KGD model


基本方程为

L=0.841 3${\left[\frac{8G{Q}^{3}}{(1-\nu )\mu }\right]}^{1/6}$t2/3
W0=0.44${\left[\frac{8(1-\nu ){Q}^{3}\mu }{G}\right]}^{1/6}$t1/3
Pw=σmin+0.242${\left[\frac{2{G}^{3}Q\mu }{{(1-\nu )}^{3}{L}^{2}}\right]}^{1/4}$

式中:L为缝长;G为切变模量;Q为泵入排量;μ为流体黏度;W0为缝宽; Pw为井底水压力;σmin为最小水平主地应力.

KGD模型计算参数为:G=10 GPa, ν=0.2, μ=1×10-3 Pa·s, Q=1×10-3 m3/(s·m).计算结果如图6所示.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   计算结果

Fig.6   Computed result


模型尺寸长×宽×高为20 m×2 m×10 m.水平向初始应力为1 MPa,垂向初始应力为2 MPa,上下边界的垂直向位移固定.模型的左侧施加水压力.

长度、宽度和孔口压力模拟结果和理论计算结果对比,具有很好的一致性,如图7~9所示.KGD模型在水平面上采用的平面应变假设,不适用于长裂缝,所以当裂缝长度较长时,模拟结果和解析解逐渐偏离.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   裂缝长度-时间曲线

Fig.7   Comparison of fracture versus time curves


图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   孔口缝宽-时间曲线

Fig.8   Comparison of orifice slit width versus time curves


图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   孔口水压-时间曲线

Fig.9   Comparison of orifice water pressure versus time curves


2.3 算例3

Chang[30]通过改变颗粒材料和浆液的性质、边界条件、初始应力状态、注入量和注入速率等控制参数,进行了一系列试验,研究无黏性沉积物中水力劈裂的物理机制.采用本文提出的方法对其试验进行数值模拟,结果如下.

图10~12是模拟与试验结果的对比.针对硅粉、红黏土、细砂和硅粉混合物3种材料的数值模拟的浆脉形状均与试验结果相符,分别形成了叶片状、单一片状和团状浆脉,进一步验证了该数值模拟方法的适用性.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   硅粉中的劈裂

Fig.10   Splitting in silicon powder


图11

新窗口打开| 下载原图ZIP| 生成PPT

图11   红黏土中的劈裂

Fig.11   Splitting in red clay


图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   细砂和硅粉混合物中的劈裂

Fig.12   Splitting in the mixture of fine sand and silica fume


2.4 算例4

劈裂注浆是通过施加压力把浆液注入土层中以改善土层性质,在注浆过程中,浆液压力使周围地层发生裂缝,劈入土体中的浆体便形成了加固土体的网络或骨架.

在岩土工程领域,劈裂注浆技术虽然被大规模采用,但对浆脉扩散和分布的研究一直没有多大的进展.本文对劈裂注浆的过程进行模拟.

在数值模拟中,土体应力应变本构关系采用邓肯-张EB模型[31],材料参数如表1所示.表中:C为土体凝聚力;ϕ为内摩擦角;K为土体的固结指数;Ku为极限剪切模量;ns为应变硬化系数;Rf为屈服函数系数;Kb为初始刚度;mc为非线性系数.

表1   土体的邓肯-张E-B模型参数

Tab.1  Parameters of E-B model for soil

参数取值
C/kPa60
ϕ/(°)30
K50
Ku100
ns0.6
Rf0.95
Kb/kPa-120
mc0.5

新窗口打开| 下载CSV


此外,土体的密度取为1 800 kg/m3,渗透系数取为1×10-7 m/s,拉伸强度取为0.6 kPa.

模型尺寸长×宽×高为6 m×6 m×10 m,计算规模为76万个自由度.设4个注浆孔,顺次注浆.每个注浆孔4个出浆口,由下往上后退式注浆,每个出浆口注浆10 min.

因为浆脉的三维形状较复杂,为了尽可能清楚呈现立体形态,每张图以透视图和剖面图组合的方式表示,如图13所示.不同颜色代表在不同出浆口注浆所形成的浆脉.由图13可以看出各出浆口产生的浆脉形成了复杂的分支状结构.为更清晰展示每个注浆孔所形成浆脉的几何特征,将4个注浆孔所形成的浆脉单独展示如图14~17所示.

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13   注浆孔位置和整体浆脉透视图

Fig.13   Perspective view of grouting hole position and overall grouting vein


图14

新窗口打开| 下载原图ZIP| 生成PPT

图14   1号注浆孔形成的浆脉

Fig.14   Grout vein formed by grouting hole No.1


图15

新窗口打开| 下载原图ZIP| 生成PPT

图15   2号注浆孔形成的浆脉

Fig.15   Grout vein formed by grouting hole No.2


图16

新窗口打开| 下载原图ZIP| 生成PPT

图16   3号注浆孔形成的浆脉

Fig.16   Grout vein formed by grouting hole No.3


图17

新窗口打开| 下载原图ZIP| 生成PPT

图17   4号注浆孔形成的浆脉

Fig.17   Grout vein formed by grouting hole No.4


在水力劈裂过程中,每个张开的裂缝对围岩、临近的裂缝产生的附加应力场被称为“应力影”.多个临近的水力劈裂裂缝之间相互作用,这种“应力影”效应称为“应力干扰”,由于应力状态复杂,每条裂缝不仅影响自身的浆脉形态,也能影响临近裂缝的浆脉形态.

图14中,1号孔水平向和竖直向都呈现多分支,以水平向更明显.初始条件下水平向的应力均布,因而水平向分支均匀分布,且每条分支的扩展方向基本不变化.而在垂直方向上呈螺旋形扩展,类似于风扇的叶片,表明劈裂过程造成小主应力围绕着垂向的旋转,表现出复杂的应力变化.

图15中,2号孔以竖直向扩展为主.水平向基本是单支分布,且其扩展方向指向1号孔.这是因为1号孔注浆对周围土体产生挤压,使得1号孔径向方向的水平应力增大,所以2号孔的浆脉在水平方向上沿着1号孔的径向扩展.在接近地表的部分,出现实际工程中经常出现的“串浆”现象,即2号孔的裂缝与1号孔的裂缝导通,并在1号孔浆脉的扩展方向上继续劈裂.

图16中,3号孔浆脉的扩展趋势和2号孔类似.但其浆脉的扩展范围要小于2号孔.表明经过1号和2号孔的挤压后,浆脉更难扩展.

图17中,4号孔以水平向扩展为主.受到1号、2号和3号孔的挤压,大主应力方向已变为水平向.1号孔浆脉的扩展没有明显的优势方向,2号和3号孔以竖向为主要扩展方向,4号孔以水平向为主要扩展方向,说明前序注浆可以对后序注浆浆脉的扩展产生影响.

4个注浆孔的注浆压力随时间变化曲线如图18所示,具体数据如表2所示.

图18

新窗口打开| 下载原图ZIP| 生成PPT

图18   注浆压力-时间曲线

Fig.18   Grouting pressure versus time curve


表2   出浆口起裂压力和终压统计表

Tab.2  Statistics of slurry outlet initiation pressure and final pressure

出浆口起裂压力/kPa,终压/kPa
1号孔2号孔3号孔4号孔
a175,120181,129181,123188,143
b169,109146,115159,118209,137
c152,113157,119166,143186,126
d141,103145,110172,137209,122

新窗口打开| 下载CSV


对同一注浆孔来说,随着埋深的减少,起裂压力和终压都减小.前期的注浆提高了地层中的压应力,并挤密了地层,使得后期的起裂压力和终压逐渐升高.相同深度时,后序孔注浆压力大于前序注浆压力.

2.5 算例5

土体参数同算例4.模型尺寸长×宽×高为20 m×20 m×10 m.注浆孔的间排距都设为2 m.共81个注浆孔,324个出浆口,计算规模640万自由度.注浆顺序为先外后内.

浆脉形成复杂的交错形状,如图19所示.为了观察浆脉的形状,挖去中部的土体,保留3个相互垂直的剖面,如图20所示.

图19

新窗口打开| 下载原图ZIP| 生成PPT

图19   浆脉透视图

Fig.19   Perspective view of pulp pulse


图20

新窗口打开| 下载原图ZIP| 生成PPT

图20   开挖揭示浆脉分布图

Fig.20   Distribution of grout veins revealed by excavation


地表2 m范围内整体向上变形,其他范围向上和向下交错,大致规律是在注浆孔附近出浆口以上土体向上变形,出浆口以下土体向下变形,如图21所示.

图21

新窗口打开| 下载原图ZIP| 生成PPT

图21   开挖面上垂直位移云图

Fig.21   Vertical displacement contour on the excavation surface


图22为两注浆孔之间沿深度竖向变形分布及对应的浆脉分布.图中:蓝色折线对应绿色虚线处的数据.可以看出地表2 m以下为压缩层,地表2 m以内抬升层.水平浆脉为分界线,其上土体整体上抬,其下土体整体压缩.从抬升效果来说,水平浆脉作用明显,且其下要有持力层.

图22

新窗口打开| 下载原图ZIP| 生成PPT

图22   沿深度竖向变形分布及对应的浆脉分布

Fig.22   Vertical deformation distribution along the depth and the corresponding slurry vein distribution


图23为注浆前后小主应力云图,可以看出小主应力的大小有明显的提升,但由于抬升作用,土体表面出现小的拉应力.

图23

新窗口打开| 下载原图ZIP| 生成PPT

图23   注浆前后小主应力

Fig.23   Small principal stress before and after grouting


3 结语

本文开发了三维有限元程序,能够在个人计算机上实现较大规模的土体中水力劈裂的数值模拟,对工程尺度的算例进行计算.通过采用GJ预处理的SQMR迭代法,联合CPU串行计算平台和GPU并行计算平台进行水力劈裂有限元数值模拟求解,提高了计算速度.采用了CSC稀疏矩阵存储法减小了对内存的需求.因而,大幅提高了在个人计算机上的计算规模,可达到百万自由度,为研究土体中水力劈裂提供了数值模拟工具.通过算例的分析对比,验证了程序上的计算能力和正确性.

参考文献

LI K R, QI C Z, WANG M Y, et al.

Research on the influence of rock fracture toughness of layered formations on the hydraulic fracture propagation at the initial stage

[J]. Geohazard Mechanics, 2024, 2(2): 121-130.

DOI:10.1016/j.ghm.2024.03.004      [本文引用: 1]

Deep underground rocks exhibit significant layered heterogeneity due to geological evolution and sedimentation. Rock fracture toughness, as one of the important indicators of hydraulic crack propagation, also exhibits heterogeneous distribution. In order to investigate the influence of non-uniform fracture toughness of layered rocks on hydraulic crack propagation, this paper establishes a planar three-dimensional hydraulic crack propagation model. The model is numerically solved using the 3D displacement discontinuity method (3D-DDM) and the finite difference method. The calculation results indicate that when the distribution of the fracture toughness of layered rocks changes from uniform to non-uniform, the fracture morphology develops from a standard circular crack to an elliptical crack. When the difference of the rock fracture toughness between adjacent rock layers and the middle rock layer (pay zone) is large enough, the fracture morphology will develop towards a rectangular shape. In addition, when the fracture toughness of rock layers is non-uniformly distributed, the hydraulic crack not only rapidly expand in the softening layer (rock layer with lower fracture toughness), but also slowly propagate in the strong layer (rock layer with higher fracture toughness). However, the propagation speed in the softening layer is much faster than that in the strong layer. The results indicate that the heterogeneity of rock fracture toughness has an important impact on the morphology, propagation speed, and direction of hydraulic fractures.

AZAD M, GHAEDI M, FARASAT A, et al.

Case study of hydraulic fracturing in an offshore carbonate oil reservoir

[J]. Petroleum Research, 2022, 7(4): 419-429.

DOI:10.1016/j.ptlrs.2021.12.009      URL     [本文引用: 1]

LI X Y, WANG Y, LI X S, et al.

Study on the production of gas hydrates and underlying free gas by horizontal well under different directions of hydraulic fracturing

[J]. Energy, 2024, 290: 130199.

DOI:10.1016/j.energy.2023.130199      URL     [本文引用: 1]

TAN P, PANG H W, JIN Y, et al.

Experiments and analysis of hydraulic fracturing in hot dry rock geothermal reservoirs using an improved large-size high-temperature true triaxial apparatus

[J]. Natural Gas Industry B, 2024, 11(1): 83-94.

DOI:10.1016/j.ngib.2024.01.002      URL     [本文引用: 1]

LIANG J S, DU X M, FANG H Y, et al.

Intelligent prediction model of a polymer fracture grouting effect based on a genetic algorithm-optimized back propagation neural network

[J]. Tunnelling and Underground Space Technology, 2024, 148: 105781.

DOI:10.1016/j.tust.2024.105781      URL     [本文引用: 1]

LI G, ZHAO Y, HU L H, et al.

Simulation of the rock meso-fracturing process adopting local multiscale high-resolution modeling

[J]. International Journal of Rock Mechanics and Mining Sciences, 2021, 142: 104753.

DOI:10.1016/j.ijrmms.2021.104753      URL     [本文引用: 1]

WANG F, LIU W, DENG J G, et al.

Hydraulic fracture propagation research in layered rocks based on 3D FEM modeling and laboratory experiments

[J]. Geoenergy Science and Engineering, 2024, 234: 212670.

DOI:10.1016/j.geoen.2024.212670      URL     [本文引用: 1]

程少振, 陈铁林, 郭玮卿, 等.

土体劈裂注浆过程的数值模拟及浆脉形态影响因素分析

[J]. 岩土工程学报, 2019, 41: 484-491.

[本文引用: 1]

CHENG Shaozhen, CHEN Tielin, GUO Weiqing, et al.

Numerical simulation of fracture grouting and influencing factors for morphology of grout veins

[J]. Chinese Journal of Geotechnical Engineering, 2019, 41: 484-491.

[本文引用: 1]

CHEN T L, ZHANG L Y, ZHANG D L.

An FEM/VOF hybrid formulation for fracture grouting modelling

[J]. Computers and Geotechnics, 2014, 58: 14-27.

DOI:10.1016/j.compgeo.2014.02.002      URL     [本文引用: 5]

CHEN T L, PANG T Z, ZHAO Y, et al.

Numerical simulation of slurry fracturing during shield tunnelling

[J]. Tunnelling and Underground Space Technology, 2018, 74: 153-166.

DOI:10.1016/j.tust.2018.01.021      URL     [本文引用: 1]

DONG Y, TIAN W, LI P C, et al.

Numerical investigation of complex hydraulic fracture network in na-turally fractured reservoirs based on the XFEM

[J]. Journal of Natural Gas Science and Engineering, 2021, 96: 104272.

DOI:10.1016/j.jngse.2021.104272      URL     [本文引用: 1]

ZHOU Y, YANG D S.

A fast simulation method for hydraulic-fracture-network generation in fractured rock based on fully coupled XFEM

[J]. Computers and Geotechnics, 2022, 150: 104892.

DOI:10.1016/j.compgeo.2022.104892      URL     [本文引用: 1]

WANG Z Y, LIANG W G, LIAN H J, et al.

Numerical study of multiple hydraulic fractures propagation in poroelastic media based on energy decomposition phase field methods

[J]. Computers and Geotechnics, 2024, 170: 106259.

DOI:10.1016/j.compgeo.2024.106259      URL     [本文引用: 1]

SONG X K, LIU Y T, YANG X W, et al.

Quanti-fying the impact of geological and construction factors on hydraulic fracture dynamics in heterogeneous rock layers using the phase field method

[J]. Computers and Geotechnics, 2024, 169: 106223.

DOI:10.1016/j.compgeo.2024.106223      URL    

WANG F Y, ZHOU M L, SHEN W Q, et al.

Fluid-solid-phase multi-field coupling modeling method for hydraulic fracture of saturated brittle porous materials

[J]. Engineering Fracture Mechanics, 2023, 286: 109231.

DOI:10.1016/j.engfracmech.2023.109231      URL     [本文引用: 1]

YU S Y, ZHOU Y, YANG J, et al.

Hydraulic fracturing modelling of glutenite formations using an improved form of SPH method

[J]. Geoenergy Science and Engineering, 2023, 227: 211842.

DOI:10.1016/j.geoen.2023.211842      URL     [本文引用: 1]

FENG R F, FOURTAKAS G, ROGERS B D, et al.

A general smoothed particle hydrodynamics (SPH) formulation for coupled liquid flow and solid deformation in porous media

[J]. Computer Methods in Applied Mechanics and Engineering, 2024, 419: 116581.

DOI:10.1016/j.cma.2023.116581      URL     [本文引用: 1]

LI M L, WU J F, LI J F, et al.

Modeling of hydraulic fracturing in polymineralic rock with a grain-based DEM coupled with a pore network model

[J]. Engineering Fracture Mechanics, 2022, 275: 108801.

DOI:10.1016/j.engfracmech.2022.108801      URL     [本文引用: 1]

TOMAC I, GUTIERREZ M.

Coupled hydro-thermo-mechanical modeling of hydraulic fracturing in quasi-brittle rocks using BPM-DEM

[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2017, 9(1): 92-104.

DOI:10.1016/j.jrmge.2016.10.001      URL     [本文引用: 1]

ABUAISHA M, EATON D, PRIEST J, et al.

Hydro-mechanically coupled FDEM framework to investigate near-wellbore hydraulic fracturing in homogeneous and fractured rock formations

[J]. Journal of Petroleum Science and Engineering, 2017, 154: 100-113.

DOI:10.1016/j.petrol.2017.04.018      URL     [本文引用: 1]

CUI W J, LIU Q S, WU Z J, et al.

Fracture deve-lopment around wellbore excavation: Insights from a 2D thermo-mechanical FDEM analysis

[J]. Engineering Fracture Mechanics, 2024, 295: 109774.

DOI:10.1016/j.engfracmech.2023.109774      URL    

NI T, PESAVENTO F, ZACCARIOTTO M, et al.

Hybrid FEM and peridynamic simulation of hydraulic fracture propagation in saturated porous media

[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 366: 113101.

DOI:10.1016/j.cma.2020.113101      URL     [本文引用: 1]

刘春, 乐天呈, 施斌, 等.

颗粒离散元法工程应用的三大问题探讨

[J]. 岩石力学与工程学报, 2020, 39: 1142-1152.

[本文引用: 2]

LIU Chun, LE Tiancheng, SHI Bin, et al.

Discussion on three major problems of engineering application of the particle discrete element method

[J]. Chinese Journal of Rock Mechanics and Engineering, 2020, 39: 1142-1152.

[本文引用: 2]

KAR S, CHAUDHURI A, SINGH A, et al.

Phase field method to model hydraulic fracturing in satura-ted porous reservoir with natural fractures

[J]. Engineering Fracture Mechanics, 2023, 286: 109289.

DOI:10.1016/j.engfracmech.2023.109289      URL     [本文引用: 1]

覃金帛, 曾志强, 梁藉, 等.

GPU并行优化技术在水利计算中的应用综述

[J]. 计算机工程与应用, 2018, 54(3): 23-29.

DOI:10.3778/j.issn.1002-8331.1711-0336      [本文引用: 1]

水利计算是水利规划、设计以及运行的基础，提高水利计算的效率对水利信息化和水资源管理具有重要实践意义。通过技术对比发现，GPU（Graphics Processing Unit）并行优化技术是性价比较高的提速策略。系统概述了GPU并行优化技术在水利计算中的应用进展；简要介绍了当前应用较多的几种并行技术；建设性提出了该项技术在水库调度、中长期水文预报和水文模型计算中的应用前景和优势；详细总结了应用该项技术的一般方法，为技术推广提供指导。最后从学科发展和应用需求的角度，有针对性的提出了技术应用难点和今后发展趋势，以期为GPU并行优化技术在水利计算中的应用提供借鉴。

QIN Jinbo, ZENG Zhiqiang, LIANG Ji, et al.

Review of application GPU technology in hydraulic parallel optimization calculation

[J]. Computer Engineering and Applications, 2018, 54(3): 23-29.

DOI:10.3778/j.issn.1002-8331.1711-0336      [本文引用: 1]

Water conservancy calculation is the basis of water conservancy planning, design and operation. Improving the efficiency of water conservancy calculation is of great practical significance for water informatization and water resources management. The technical comparison shows that GPU（Graphics Processing Unit） parallel optimization technology is a cost-effective speed-up strategy. This paper systematically summarizes the application progress of GPU parallel optimization in water conservancy calculation, briefly introduces several parallel techniques currently used in the proposed calculation; constructively puts forward the application prospects and advantages of this technology in reservoir scheduling, mid-and?long term hydrologic forecasting and hydrological model calculation; summarizes the general methods with the application of this technology and provides guidance for technology promotion. Finally, from the perspective of discipline development and application demand, the difficulties and future trends of technology application are pointed out, in order to provide reference for the application of GPU parallel optimization technology in water conservancy calculation.

李路涛. GPU加速的大规模岩土工程有限元计算中的迭代求解[D]. 北京: 北京交通大学, 2012.

[本文引用: 2]

LI Lutao. GPU accelerated iterative solution of large-scale geotechnical finite element computations[D]. Beijing: Beijing Jiaotong University, 2012.

[本文引用: 2]

PHOON K K, TOH K C, CHAN S H, et al.

An efficient diagonal preconditioner for finite element solution of Biot’s consolidation equations

[J]. International Journal for Numerical Methods in Engineering, 2002, 55(4): 377-400.

DOI:10.1002/nme.v55:4      URL     [本文引用: 2]

FREUND R W.

Quasi-kernel polynomials and their use in non-Hermitian matrix iterations

[J]. Journal of Computational and Applied Mathematics, 1992, 43(1/2): 135-158.

DOI:10.1016/0377-0427(92)90263-W      URL     [本文引用: 1]

GEERTSMA J, DE K F.

A rapid method of predicting width and extent of hydraulically induced fractures

[J]. Journal of Petroleum Technology, 1969, 21(12): 1571-1581.

DOI:10.2118/2458-PA      URL     [本文引用: 1]

With the design charts presented here, and nothing more elaborate than aslide rule, it is possible to predict the dimensions of either a linearly or aradially propagating, hydraulically induced fracture around a wellbore.

CHANG H. Hydraulic fractures in particulate materials[D]. Atlanta, USA: Georgia Institute of Technology, 2004.

[本文引用: 1]

DUNCAN J M, CHANG C Y.

Nonlinear analysis of stress and strain in soils

[J]. Journal of the Soil Mechanics and Foundations Division, 1970, 96(5): 1629-1653.

DOI:10.1061/JSFEAQ.0001458      URL     [本文引用: 1]

/

〈 〉