完整井与非完整井抽水引起的地裂缝活动量计算模型

2022-11-29郭乐凡刘雪弟成盛林韩宁博

河北地质大学学报 2022年6期

郭乐凡, 刘雪弟, 刘超,2,3, 成盛林, 韩冰, 韩宁博

1. 河北地质大学城市地质与工程学院, 河北石家庄 050031; 2. 河北省高校生态环境地质应用技术研发中心, 河北石家庄 050031; 3. 河北省地下人工环境智慧开发与管控技术创新中心, 河北石家庄 050031

0 引言

地裂缝是指浅表岩土体在地球内外动力作用和人类活动作用下发生破裂或者错断的地质现象。中国作为世界上地裂缝灾害最为严重的国家之一, 我国地裂缝主要分布在汾渭盆地、河北平原、苏锡常平原等地[1,2]。截止到2015 年, 已在全国共28 个省、直辖市的4 000 余处累计发现地裂缝 5 000 余条[2]。地裂缝在人类生活中会对房屋建筑物[3]、农田等造成毁坏, 在工程建设中地裂缝对隧道、地铁[4]、大坝等工程的施工造成严重的影响。

关于地裂缝的成因机制大体有3 种说法: 地下水开采成因机制、构造成因机制、复合成因机制[5]。人为过量开采地下水是诱发地裂缝活动的一个重要因素。在抽水作用下, 地裂缝上下两盘产生沉降差异[6,7], 诱发地裂缝活动, 抽水引起地裂缝活动是地层中渗流场和应力场之间复杂的耦合过程[8-10]。

模拟抽水引起变形的方法主要有“两步走”、部分耦合、完全流固耦合3 大类思路[11]。其中“两步走” 模型是将土体变形和渗流两个过程分别考虑, 而部分耦合模型则在“两步走” 模型基础上考虑到计算参数的变化, 然而力学模型和渗流模型仍是分开模拟, 而完全耦合模型是基于Biot 固结理论考虑土体变形和地下水流动的相互作用, 将土体变形和地下水流动模型统一于相同的物理空间[12], 因此以比奥固结理论为基础的三维完全耦合模型是较为完善的模型[13]。

本文以抽取承压含水层引起地裂缝的机制为基础, 基于理想的地层模型, 通过FLAC3D6.0 软件编码程序, 采用三维完全耦合的方法, 模拟不同流量下抽水引起的地裂缝活动, 同时对比和讨论了完整井和非完整井抽水对地裂缝活动的影响。

1 基本理论

FLAC3D6.0 软件在模拟计算时所采用的数学力学计算原理为有限差分法, 即将空间离散点处的控制方程组中每一个导数直接由位移、应力等含场变量的代数表达式替换, 因此可以在每一步重新生成有限差分方程, 通过“显式的” 时程方法求解代数方程。 “显式的” 时程方法首先根据运动方程由应力和外力导出新的速度和位移, 进而得到应变速率和新的应力, 构成下一个循环的初始状态。

那么采用FLAC3D6.0 软件模拟三维完全耦合计算时, 需要考虑平衡方程、边界条件方程、流量平衡方程、应力应变方程等。

取土中微小单元体, 则平衡微分方程为:

其中,ρ为密度,σij为应力分量,u·i为速度分量,t为时间,x为方向。

在设置边界应力参数时, 边界条件方程为:

其中ni为边界段外法线方向单位矢量, Δs为应力,作用在边界段的长度,Fi为施加在边界的力。在进行完全耦合模拟时, 水的渗流满足达西定律, 即:

其中,q为流量,k为渗透系数,A为过水面积,i为水力坡降。

再取微单元体为研究对象, 则整体的渗流平衡方程为:

其中,dQ为水的体积变化量,q为单位时间内通过单元体的流量,z表示竖直方向。

本文所建模型为具有各向同性的弹性连续性介质摩尔—库伦地层模型, 则在所建本构模型中, 应力应变依据胡克定律, 用应变增量形式表示

其中:

αi为比奥固结系数,K为体积模量,G为剪切模量,Δσij为应力,Δeij为应变。

而其中主应力和主应变的表达式为

其中Δσi为i方向的主应力,为i方向的弹性主应变。

2 计算模型

如图1 和图2, 选取计算模型, 建立抽取承压水引起地裂缝活动计算模型。计算模型中承压含水层为10 m 厚, 隔水层为40 m 厚, 整体土体模型为50 m厚, 计算模型的长度和宽度均为100 m。地裂缝倾角约80°, 上盘分布两个土层, 上层为粉质黏土隔水层,下层为砂土承压含水层, 下盘为粉质黏土隔水层, 抽水井设置在模型中心, 在地面距地裂缝出露处约30 m, 抽水井选取为完整井和非完整井两种类型, 选取不同的流量工况, 见表1。

图1 完整井模型Fig.1 The model of complete well

图2 非完整井模型Fig.2 The model of incomplete well

表1 模拟工况表Table 1 Simulation conditions

各土层本构模型选为弹性模型, 利用FLAC3D6.0编写代码建立模型, 如图3 所示, 模型参数见表2。

图3 计算模型图Fig.3 The model of calculation diagram

表2 模型参数表Table 2 Model parameters

模型四周为定水头边界, 底面为隔水边界, 力学边界条件中, 上边界为自由边界, 四周和下边界固定, 计算模型地裂缝错断了含水层, 设为不透水的结构面单元。

计算模型首先生成初始应力场和孔压场, 含水层顶面初始承压水头为50 m, 承压含水层的渗透系数远大于隔水层。

3 计算结果及分析

3.1 完整井抽水计算结果及分析

3.1.1 孔压分析

计算模拟的结果分别考虑隔水层顶面孔压变化、地表面变形, 以及竖直方向上变形。

(1) 含水层顶面孔压分布

在承压含水层顶面, 孔压分布的云图见图4 (以工况2 为例)。

图4 工况2 承压含水层顶面孔压Fig.4 The top pressure of confined aquifer of condition 2

图5 为承压含水层顶面孔压分布, 图4 和图5 反映出地裂缝的上盘承压含水层顶面, 抽水引起了孔压降落漏斗, 漏斗中心在抽水井附近, 降落漏斗中心的孔压随抽水流量的增大而降低。地裂缝的隔水作用使得下盘孔压分布与上盘不连续, 地裂缝的下盘孔压从地裂缝向远处略微增大, 从云图也可以反映出, 下盘孔压分布也呈漏斗形, 漏斗中心处大致在地裂缝附近, 随着上盘抽水流量的增加, 孔压略微增加。

图5 承压含水层顶面孔压Fig.5 Top pressure of confined aquifer

(2) 竖直方向孔压分布

在竖直方向孔压分布的云图见图6 (以工况2 为例)。

图6 工况2 抽水井处竖直方向孔压Fig.6 The pore pressure in vertical direction at pumping well of condition 2

图7 为地裂缝上盘沿抽水井竖直方向上的孔压变化, 由图6 和图7 可以看出, 由于地裂缝的隔水作用, 因此地裂缝两侧孔压不连续, 上盘的分布特征与没有地裂缝时类似。从云图中可以看出下盘的孔压变化量小于上盘, 孔压沉降中心位于地裂缝附近。

图7 抽水井处竖直方向孔压Fig.7 The pore pressure in vertical direction at pumping well

3.1.2 变形分析

(1) 地表竖直方向沉降分析

图8 显示了地面沉降分布的云图(以工况2 为例)。

图8 工况2 地面沉降变形Fig.8 Ground settlement and deformation of condition 2

图9 为沉降变形曲线, 从图8 和图9 中可以看出, 地裂缝上盘地裂缝附近沉降量最大, 随着距地裂缝的距离增加, 沉降量逐渐减小, 地裂缝的下盘略向上发生隆起变形。以上下盘沉降差作为地裂缝的活动量, 在抽水量400 cm3/s 时, 地裂缝活动量约为29 mm。

图9 地面沉降变形Fig.9 The settlement and deformation of ground

(2) 地表水平方向变形分析

图10 显示了地表水平方向变形的云图(以工况2为例)。

图11 为地表水平x方向的变形曲线, 从图10 和图11 中可以看出, 由于地裂缝的存在, 上下两盘水平方向的变形不连续, 下盘土体在水平方向上主要受压力作用, 而上盘则表现为受拉和受压区域相间分布, 随着抽水流量的增加, 变形量在逐渐增加。

图10 工况2 地面水平x 方向位移Fig.10 The ground horizontal X-direction displacement of condition 2

图11 地面水平x 方向位移Fig.11 The ground horizontal X-direction displacement

图11 显示, 在地裂缝的下盘, 土体向着背离地裂缝方向挤压变形, 地裂缝附近的变形量最大, 约为3.6 mm, 在地裂缝上盘, 在距地裂缝约5 m 的范围内, 土体向着背离地裂缝的方向挤压变形, 在地裂缝附近变形量最大, 约为1 mm, 向远处逐渐减小, 距地裂缝5 m 以外, 地表土体向着地裂缝方向运动, 随着距离的增加位移量先增大后减小, 直到恢复天然状态。地裂缝出露地表处, 上下盘均背离地裂缝方向变形, 地裂缝处于拉张状态, 加剧了地裂缝出露的地表宽度, 随着抽水量的增加裂缝宽度增大。上盘地表土体在距地裂缝60 m 范围内处于压缩状态, 在0 ～5 m范围内, 由于受到地裂缝产状的倾向上盘的影响, 地面沉降使得上盘土体背离地裂缝方向压缩, 5 ～60 m范围内, 受到抽水和地裂缝产状的影响, 水平变形向着地裂缝方向产生压缩, 压缩速率在地裂缝与抽水井之间的范围内相对较大, 在抽水井外侧略减小。在60 m 以外, 处于沉降漏斗的边缘区域, 随着距离的增加, 变形量逐渐减小, 地表土体处于拉伸状态。

(3) 竖直方向变形分析

图12 为竖直方向上变形的云图(以工况2 为例)。为了能够明显显示出上下盘差异, 竖直方向的位移放大了100 倍, 从云图中明显看出, 地裂缝上下两盘发生了错动, 上盘向下运动, 下盘相对上盘上升, 下盘地表处略隆起, 随着抽水流量的增加, 上下两盘的差异变形量增加, 地裂缝的活动量加剧。

图12 工况2 抽水井处竖直方向位移Fig.12 Vertical displacement at pumping well of condition 2

图13 为沿抽水井竖直方向的位移曲线, 在地裂缝上盘, 隔水层的沉降量远大于承压含水层的沉降量, 土体竖直方向上的变形趋势与没有地裂缝时计算结果的变化趋势是一致的。随着抽水流量的增加, 隔水层和承压含水层的沉降量都增加。

图13 抽水井处竖直方向位移Fig.13 Vertical displacement at pumping well

3.2 非完整井抽水计算结果及分析

采用非完整井, 井在含水层中的部分为5 m, 是完整井的一半, 抽水流量为400 cm3/s, 模型其他条件不变。

承压含水层顶面孔压分布的云图见图14, 竖直方向孔压分布的云图见图15, 可以看出地裂缝上盘的孔压分布仍呈漏斗型, 与完整井的分布形式一致。

图14 非完整井工况承压含水层顶面孔压Fig.14 Top pressure of confined aquifer under incomplete well condition

图15 非完整井工况抽水井处竖直方向孔压Fig.15 Vertical pore pressure at pumping well under incomplete well condition

图16 为地面沉降分布的云图, 图17 为地表水平方向变形的云图, 图18 为竖直方向上变形的云图(放大倍数100), 从变形的云图上可以看出变化的趋势与完整井抽水是一致的, 变化的数值较完整井抽水时小。

图16 非完整井工况抽水地面沉降变形Fig.16 Ground subsidence deformation of pumping under incomplete well conditions

图17 非完整井工况抽水地面水平x 方向位移Fig.17 Horizontal X-direction displacement of pumping ground under incomplete well condition

图18 非完整井工况抽水井处竖直方向位移Fig.18 Vertical displacement at pumping well under incomplete well condition

4 完整井与非完整井抽水计算对比

图19 为沿抽水井竖直方向孔压分布曲线。从图19 可以得到以完整井和非完整井两种形式同样抽水流量进行抽水, 孔压变化的趋势的是一致的, 从数值上, 隔水层的孔压变化差别小, 而承压含水层孔压变化相差较大。完整井400 cm3/s 抽水埋深50 m 处孔压约为310 kPa, 非完整井400 cm3/s 抽水埋深50 m 处孔压约为473 kPa, 完整井引起的孔压降落大。

图19 沿抽水井竖直方向孔压分布Fig.19 Pore water pressure distribution along the vertical direction of pumping well

图20 为承压含水层顶面孔压分布曲线。从图20可以得到以完整井和非完整井两种形式同样抽水流量进行抽水, 受抽水井和地裂缝位置关系的影响, 在抽水井两侧各20 m 范围内孔压有较大的差别, 在抽水井处差别最大。

图20 承压含水层顶面孔压分布Fig.20 Top pressure distribution of confined aquifer

图21 为沿抽水井竖直方向的位移分布曲线, 从图21 可以看出, 完整井与非完整井产生的沉降变化趋势是一致的, 在隔水层内从地表随着深度的增加,变形量先增大, 后减小, 在含水层内随深度增加, 变形量减小, 完整井产生的变形量较非完整井大。在相同深度处竖直方向完整井与非完整井产生的沉降量的差值, 在承压含水层比隔水层大。在承压含水层顶面, 沉降的差值达到最大, 约为0.72 mm, 随着深度增加逐渐减小, 在隔水层内, 完整井和非完整井产生的沉降量的差值从地表向下逐渐增加。