车 平 刁志翔 王奇石 商翔宇

(1.江苏省有色金属华东地质勘查局,江苏 南京 210007;2.中国矿业大学 力学与土木工程学院,江苏 徐州 221116;3.中国矿业大学深部岩土与地下工程国家重点实验室,江苏 徐州 221116)

华东矿区是我国非常重要的煤炭能源基地。自1987年以来,该矿区相继发生井筒破裂100余处,严重影响了矿区的安全生产。为此,学者们围绕立井井筒的破坏机理开展了大量研究。目前,附加应力破坏理论已经得到广泛认可[1]:即底部含水层疏水引发周围土层固结和沉降,沉降地层在立井井筒上施加竖直向下的附加应力,附加应力与井筒原有的应力叠加,导致井筒破裂。在此基础上,学者们开发出多种有效技术措施来防治井壁破裂[2,3]。然而,经修复加固后的破裂井壁仅能在不确定的一段时间内安全运行,并且几乎不可避免地会再次发生破裂[4]。归根结底,这是由于井筒周围地层长期持续沉降造成的。

已有大量研究表明,底部含水层疏水是造成立井井筒周围深厚地层沉降的主要原因[1,5,6]。也有关于立井井筒周围不同深度地层变形对地表长期沉降贡献的研究。李文平和于双忠分析认为,深埋黏土天然含水量小于吸附结合水含量,因而可能不具有压缩性[7];但李定龙分析认为,底部含水层上、下的微弱隔水黏土层中部分弱结合水会随底部含水层和下覆裂隙基岩水位的下降而逐渐释放[6]。王秀艳等试验研究指出,深层黏土中以结合水为主,但仍能释水压缩[8];许延春等[5]的实测结果指出,底部含水层深部砂土及相邻深埋黏土的压缩变形是地表沉降的主要来源。柴敬和杨金宏等[9,10]基于光纤技术的实测数据也表明,厚度较大的深部砂土和黏土层均压缩明显。由以上文献分析可知,深部黏土地层的压缩变形很可能对井筒周围地表长期沉降有贡献,但具体贡献占比如何? 伴随着井壁的局部破裂,其在井筒整个服役周期内将如何变化? 这些问题的正确解答对于井筒治理工作具有重要意义[4],但目前尚未有相关研究报道。

笔者拟利用模拟方法,对华东矿区孔庄煤矿主井井筒周围地表长期沉降进行模拟对比研究,以期初步回答前述问题,所获研究结论对破裂立井井筒的长期安全管理有借鉴意义。

1 立井井筒周围地表长期沉降数值模型

1.1 孔庄主井基本地质条件

孔庄煤矿主井于1975年4月投入生产,井筒深度465.8 m,井筒净直径5.00 m。根据实际统计资料[14],贯穿主井井筒的第四系表土层厚度约为156 m,主要由黏土、砂质黏土及砂层组成。根据实际情况及研究需要,为简化计算,将相似相近地层进行合并,从而将实际地层简化为7个含水层及6 个弱透水层,其组成、厚度及埋深如表1所示。

表1 简化后的地层分布

1.2 数值模型采用的简化假设

由于该井筒建设年代较远,缺少必要的土层力学和水力试验结果,本文假设:①各弱透水层和含水层分别假设为同类均质黏土和砂土,除了先期固结压力和初始孔隙比不同外,其他力学参数均相同,渗透系数由孔隙比决定;②假定含水层及弱透水层均为正常固结状态的饱和土体;③假设模拟土层的力学性质和渗流性质均为各向同性,土中水的运动符合达西定律;④忽略地表建筑对地面沉降的影响;⑤根据该煤矿水文地质资料,仅考虑地层静水压力;⑥主要关注井筒周围地层的长期沉降变形,不考虑大气降水对浅部含水层孔隙水压的波动影响[1]。

1.3 数值模型的建立

计算基于商用有限元软件ABAQUS 进行,模型按常用的比奥固结理论,考虑渗流固结耦合作用;采用空间轴对称模型,选用二阶精度的减缩积分位移-孔压耦合单元;模型高度取为156 m(与实际表土层厚度相同),模型宽度取为50 m(根据最远沉降观测点距井筒的距离),因模型网格划分相对规则,故不再给出图示。参考该煤矿设计生产年限,模型总模拟时间为100 a。

需要确定表1所示土层的力学参数以及渗流参数。土体采用修正剑桥模型,其中黏土力学参数见表2,系由前期同地区深部黏土的实验结果确定[15]。其中,λ是临界状态线在v-lnp′平面(v为比容,p′为有效平均应力)上的斜率,κ是卸载曲线在v-lnp′平面上的斜率,M是临界状态线在p′-q平面(q为等效偏应力)上的斜率,e1是正常固结曲线在p′=1 kPa时的孔隙比,β是控制p′-q平面上屈服面形状的参数,K是控制偏平面上屈服面形状的参数,ν是泊松比。

表2 黏土力学参数

模型的砂土力学参数则参考文献中典型标准砂的临界土力学参数[16],具体见表3。

表3 砂土力学参数

由于不同深度土层的渗透系数不同,因此,考虑渗透系数随孔隙比的非线性变化,具体变化规律采用作者前期的实验结果[17]:

SPME获得的样品于安捷伦GC*GC-Q-TOF(7890B/7200型)进行分析。调制周期：3 s；热喷时间：350 ms。一维色谱柱为HP-5ms，30 m×0.25 mm×0.25 μm；二维色谱柱为DB-17HT，1.8 m×0.1 mm×0.1 μm。溶剂延迟时间：3 min，不分流进样，载气为He，流速：1 mL/min。进样口温度为250 ℃。升温程序为40 ℃保持2 min，5 ℃/min升至270 ℃，保持10 min，共58 min。质谱为EI源，接口温度：280 ℃，扫描范围：40～700 m/z。

黏土的渗透系数由公式(1)计算,而砂土则参照经验,取有相同孔隙比黏土渗透系数的3倍。

模型的定界条件包括初始条件和边界条件。其中模型初始条件包括各地层的初始竖向有效应力、初始孔隙比、初始孔压和初始干密度等,根据工程实际地质情况以及前期深部土压缩试验结果[17]确定,具体见表4。

表4 各地层的初始条件

模型边界条件包括位移边界条件以及水力边界条件,其中,位移边界条件较为常规,模型上边界竖向可自由移动,而水平方向固定;模型底部设定为刚性基座,水平与垂直方向均不可移动;模型左、右边界垂直方向可移动,水平方向固定。水力边界条件对地表长期沉降模拟相当关键,因此,下文着重介绍模型水力边界的考虑与确定。

模型上部水力边界较易确定,设为排水边界,而底部水力边界则需要根据工程实际区分不同情况。首先,工程实测资料表明,该煤矿底部含水层水位多年来一直持续下降,很可能是造成地层沉降的主要因素,因此,仅将底部设为孔压降低边界条件,记为边界条件1。其次,该煤矿主井井筒内实际淋水观测表明,靠近底部含水层的上临含水层很可能沿破裂井筒外壁与底部含水层形成水力联系,李定龙也曾分析指出,地层水可能在底部含水层附近产生向井筒微裂隙的渗流运动[6]。因此,本模型将最靠近底部的含水层⑥也设置为疏水地层,其水位下降量为底部含水层的1/4,将底部含水层和含水层⑥孔隙水压力同时下降,记为边界条件2。可见,上述底部水力边界条件1对应于“单纯底含失水”疏水条件,而上述水力边界条件2则对应于“底含失水+井筒破裂失水”的疏水条件。

目前,孔庄煤矿主井仅有近10 a的附近水文观测数据(见图1)。数据表明,2008年到2018年底含水层水位下降8 m,即其水头降低80 k Pa;此外,实际运营结果显示,随着井筒服役时间的不断增长,井壁损坏程度日益加大,井筒的涌水量不断增加,水位下降速度越来越快,有理由推断近10 a是水位下降最多的10 a。据此分析,水力边界条件1可由图2所示的底含上边界孔压变化曲线表示,即前40 a的水位每10 a分别下降1 m、2 m、4 m、8 m,相应底部含水层水头降低分别为:10 k Pa、20 kPa、40 k Pa和80 k Pa。水力边界条件2的具体数值可以根据前述假设由边界条件1确定。

图1 底含水位变化

图2 底含上边界孔压变化

2 模拟结果与实测数据的对比与分析

通过模型前40 a的模拟结果与已有实测数据的对比,以检验模型的合理性。

2.1 模拟结果与实测数据的对比

自2005年以来,矿方在孔庄主井井塔四周布置有4个沉降测点(即主1～主4),并在距井筒约30 m 处布置沉降测点N6,由于观测期间部分测点被破坏,目前,只有主2和N6测点近10 a的数据最为完整。由于主井井塔的基础和井筒本身相互分离,即井塔的沉降和井筒沉降无关,故监测数据可反映地层沉降。

模拟含水层疏水的第30～40 a可以大致对应为2008～2018年(目前仅有该时间段的数据),为了与监测数据进行对比,将井筒周围主2、N6测点的实际监测值和地层沉降模拟值分别绘于图3中。

图3 模拟和实测地表沉降对比

由图3可以看到,主2和N6两点的实际监测值尽管不完全重合,但其趋势基本相同,并且过去十年的总沉降量也基本相同。这意味着,实测地表沉降有相当比例源自于深部土层压缩变形的贡献,进而使得井筒-地层相互作用对上述两个测点的实际沉降变形影响差异较小,因此,两者沉降规律大致相当。

对比模拟结果与实际监测值可以看到,前述两类水力边界条件所对应的模拟沉降规律和实际沉降规律均较为接近,模拟结果与实际沉降的最大相对误差为-16.7%～+28.3%。考虑到前述所采用的模型参数均由类比法选取,而且水力边界条件是由实测水位下降推出,不可避免地会造成模拟结果和实际值存在偏差,但两者总体规律较为接近,因此,可以认为模型基本合理。由图3可以看到,近10 a内,实际监测沉降量约为60 mm,边界条件1的模拟沉降量约为50 mm,而边界条件2模拟沉降量约为87 mm,由此可见,考虑井筒破裂新增失水通道会明显增大总地表沉降量。

含水层疏水不仅会引起其自身的压缩,而且还导致相邻[5,9,13]甚至更远含水层和弱透水层的固结压缩。为了分析各土层特别是黏土层对地表沉降的贡献量,对比了这两种水力边界条件下前40 a各层的具体贡献值(如图4所示),图4中横轴h代表“含水层”,r代表“弱透水层”。可见,无论何种水力边界条件情况下,深部弱透水层即深部黏土地层对地表沉降的贡献均较大,具体而言,在边界条件1下,深部黏土地层对地表沉降的贡献值约为80%,而在边界条件2下,深部黏土地层对地表沉降的贡献值约为70%,而含水层的贡献值较低,其原因在于,含水层是由砂土组成,其埋深较深,砂土较为密实,砂土的孔隙总体较小,因此,沉降量相对较小。总而言之,砂土与黏土交替出现的深厚地层中,总的地表沉降中深部黏土层的贡献相当可观。

图4 过去40 a各地层沉降贡献值

此外,可以看到,如果只有底部含水层疏水,紧临底部含水层的弱透水层⑥会成为对地层沉降贡献最大的地层,并且与其他地层有着很大的差异,导致这种差异的主要原因是:由于弱透水层⑥的埋深相对较大,有效应力大、孔隙比小,相应的渗透系数非常低,孔隙水渗流速度极其缓慢,短期内不会影响上覆含水层和弱透水层中孔隙水的运动,更不会在短时间内进一步造成更远的上覆地层的大量沉降;而如果含水层⑥也疏水,即为边界条件2,则紧临含水层⑥的弱透水层⑤和⑥都将成为对地表沉降贡献最大的两个地层,并且与其他地层之间的差异也有逐渐缩小的趋势;可以预测,当井筒中存在多处漏水或下部多个含水层之间存在水力联系时,地层沉降速率将更快,并且每个地层的贡献将相对更加均匀。

2.2 服役周期内井筒未来地表沉降预测分析

上述模型前40 a的模拟结果与已有实测数据的对比表明了模型的合理性。以此为基础,模型后60 a的计算结果可以作为预测该井筒周围地层长期沉降规律的合理依据。

考虑一种理想的极限情况:第40 a之后,地层经过加固、井筒经过治理,底部含水层原有的以及井筒破裂新增的疏水通道均被关闭,不再疏水,即上述相应边界孔隙水压不再变化。图2 中40～100 a的水平线段给出了相应的底部含水层孔压。

模拟获得的后60 a及前40 a地层的平均沉降速率见表5,各个地层沉降贡献值见图5,总体沉降曲线见图6。可以看到,尽管含水层不再疏水,但其内部孔隙水压在重新平衡的过程中,其上覆地层的压缩潜力逐渐释放,沉降量持续增加。由于砂土渗透系数相对较大、变形较易稳定,因此,每个含水层的后期贡献值均较小,甚至之前疏水的含水层已经停止下沉,而黏土层由于渗透系数小、压缩性大,随着时间的增加,其贡献值逐渐增加。具体而言,在未来60 a内,弱透水层对地层沉降的总贡献率将达到90%。由表5 和图6可知,在水力边界条件1情况下,井筒周围地表的沉降速率在整个井筒工程服役期限内变化幅度很小;在水力边界条件2情况下,后60 a的沉降速率较前40 a的降低也不足30%,这归因于图5中所示的黏土地层对井筒周围地表沉降的贡献具有长期性。由此可见,含深厚黏土地层的长期沉降性质决定了井筒破裂治理也将是个长期、多次的过程[4]。

图5 未来60 a各地层沉降贡献值

图6 0～100 a模拟沉降曲线

表5 地层平均沉降速率

除上述计算模型外,笔者采用类似方法对徐庄煤矿主井周围地层的长期沉降变形也进行了模拟计算,所获结论与本文结果相近[18]。

3 结 论

笔者对孔庄煤矿主井周围地表长期沉降进行了数值模拟研究,围绕深立井周围地表沉降中深厚黏土层的贡献、井筒破裂新增失水通道对地表沉降的贡献等问题,开展了数值模拟,并与实测结果对比分析,获得了如下结论:

对孔庄煤矿主井的模拟结果表明,2008～2018年的井筒周围地层沉降模拟结果和实测结果趋势一致,量值相当,表明模型较为合理;模拟结果显示:深厚黏土地层在立井井筒周围地表沉降中贡献显著,前40 a内,与疏水含水层相邻的深部黏土层压缩变形显著,其压缩量占地层总沉降量的70～80%;即使是在考虑含水层不再疏水的理想极限情况下,后60 a内,地层沉降速率依然可观,且与前40 a内疏水含水层相距较远的弱透水层也产生了压缩变形,其中,黏性弱透水层的压缩变形对立井周围地表沉降的总贡献率达到了90%。此外,由于井筒破裂造成的新失水通道会加速总的地表沉降,对未来长期沉降有显著影响。