地下水抽取对溶洞稳定性影响的流-固全耦合分析

梁文鹏, 钱静, 王明龙

(1.中国科学院深圳先进技术研究院, 深圳 518000; 2.深圳市自然资源和不动产评估发展研究中心, 深圳 518000;3.深圳市地质环境监测中心, 深圳 518000)

岩溶在中国分布广泛,约占国土面积的1/10[1-2]。随着“一带一路”“西部大开发”等一系列工程实施,在岩溶区进行的交通隧道、城市化建设和地下空间开发等地下工程数量急剧增加。然而在建设过程中,地下水位受迫变动将直接影响溶洞稳定,进而诱发岩溶塌陷、破坏上层建(构)筑物[3-4]。这将对人口密集区人民群众的财产和生命安全造成巨大威胁,并逐渐成为制约中国某些区域经济快速发展的关键问题。

目前,不少学者已经利用现场试验[5-6]和室内物理模型[7-8]分析了水位变动对溶洞稳定性的影响,结果均表明抽取地下水所导致的水位变动是岩溶地区塌陷形成的主要人为因素之一[9-10]。为从理论出发定量给出溶洞稳定与水力联系之间的关系,Chen等[11]推导了降水诱发椭球形洞体岩溶塌陷的解析解,分析了地下水参数与盖层厚度的关系,表明稳定系数与地下水位下降负相关,与盖层厚度正相关。Jia等[12]基于改进的Terzaghi土压力理论,给出了塌陷时降水深度临界值。Wei等[13]建立了洞体形成、柱体垂直塌落和表面漏斗形倒塌3个阶段的力学解析解。靳红华等[14]基于极限平衡法计算分析的结果表明长时间降水是形成塌陷的主要原因。以上解析解因操作简便且参数易于获取被广泛应用于工程分析。然而其在建立过程中存在大量假设且在机理分析、时空预测和物理场演化规律分析等方面略有不足。

应用数值模拟方法研究岩溶塌陷可较为直观地计算得出各物理场的时空演化,已经被应用于岩溶区的各类研究中[15-16]。Drumm等[17-18]先后利用有限元方法计算了覆盖层厚度、不同土性参数对溶洞的影响,所提出的“稳定图”概念可较为方便地应用于实践,然而其并未考虑流场影响。同样地,Soliman等[19]忽略流场影响后,仅从力学角度通过一系列有限元模拟,评价了洞体几何形态(深度、直径)和盖层土体力学性质对洞体稳定性的影响。陈峰等[20]利用有限差分软件FLAC3D分析了稳态条件下洞径及埋深对高水位岩溶区内深基坑稳定性的影响。薛明明等[21]同样采用有限差分法模拟不同降雨量对不同洞径溶洞稳定性的影响。陈冬琴等[22]基于太沙基固结理论和GMS软件,假设土体线弹性,采用单向解耦法建立了有限元耦合分析模型。熊启华等[23]利用弹塑性理论和渗流理论建立了流-固单向耦合有限元模型。综上所述,目前部分研究中直接忽视了流场,重点从力学角度分析溶洞稳定性,部分研究虽采用有限差分法或有限元法进行解耦计算分析了溶洞稳定性,但很少研究采用全耦合计算,相关研究表明在某些情况下流固耦合项的影响不可忽视[24]。

为解决上述问题,精确量化地下水抽取过程中地面沉降及洞周应力分布,分析评价该过程中不同影响因素对溶洞稳定性的影响,现基于连续介质力学推导耦合控制方程,结合渗透率关系、应力-应变关系等本构方程给出闭合数学描述,通过经典算例验证模型有效性后,以深圳市某场地地质条件为背景,研究不同抽水速率、洞径、洞形等多个因素对沉降和溶洞稳定性的影响,分析不同工况下孔压、应力和塑性区等关键物理场的时空分布,阐述多物理场作用下岩溶塌陷及地面沉降的力学破坏机制,并利用灰色理论确定主要影响因素,以期使相关研究成果为岩溶塌陷防治提供技术参考。

1 研究区地质条件概述

研究区域位于深圳市,地处北回归线以南,珠江三角洲南端,雨量充沛,地下水位变幅1～2 m[25]。其位于在建建筑基坑东侧,在建地铁站的东北侧为隐伏岩溶区,依据收集资料,第四系厚度5～30 m,下伏基岩为石炭系石磴子组灰岩、部分灰岩大理岩化,其区域地质图如图1所示。

图1 研究区区域地质图

研究区内基岩面上分布填土、软塑状粉质黏土,基岩溶洞发育强烈,溶洞顶板灰岩较薄,其地质剖面图如图2所示。地下水主要赋存于灰岩裂隙和溶洞中,区内岩溶水稳定水位埋深1.80～5.20 m,抽水试验结果表明该区地下水富水性及联通性较好。

图2 研究区地质剖面图

从2020年初,地铁站基坑、建筑基坑开挖抽水,基坑外地下水水位不断下降,引起地下水急剧变化或持续降低。降水漏斗作用在溶洞(隙)上覆土体,使得基岩面顶部的软-流塑土层发生流失破坏,从而导致研究区内局部建筑物因地表不均匀沉降产生裂缝。根据相关资料,目前由地表不均匀沉降已经出现,施工抽水对溶洞稳定性及地表沉降的影响亟待深入分析。

2 模型建立

2.1 数学模型

根据连续介质力学理论[26],固相s的质量守恒方程为d(nsρs)/dt+nsρs∇·vs=0。其中n为体积分数,ρ为密度,vs分别为固相的速度。代入小变形假设∇·vs≈dεv/dt后可得

dn/dt=(1-n)dεv/dt

(1)

式(1)即为岩土体发生体变εv过程中孔隙率随时间t的演化方程。将流相与固相物质导数间的关系dl()/dt=d()/dt+(vl-vs)·∇()代入流相的质量守恒方程后可得流相质量守恒方程为

(2)

式(2)中:Sl为流相饱和度;ql=nl(vl-vs)为流体相对于固相的达西流速,由压力pl的梯度、孔隙介质中渗透性Kl及流体黏滞性μl确定,ql=-Kl(∇pl-ρlg)/μl,vl为流相的速度。忽略温度效应,与流体密度随压力演化方程dρl/dt=ρldpl/Bldt代入式(2)后得

(3)

式(3)中:Bl为液相的体积压缩模量,通常取值4×10-10;Kl通常与孔隙率、土体孔隙尺寸分布数αm和土体饱和状态相关[23],即

(4)

根据动量平衡方程∇·(σ)=(nsρs+nlρl)g、有效应力原理σ′=σ-pl1,结合弹塑性本构模型可求解土体的应力应变关系为

σ′=D:(L:u-εp)

(5)

式(5)中:σ、σ′分别为总应力及有效应力;1为单位矩阵;D为弹塑性刚度矩阵;u为位移;L为应变位移矩阵。

式(1)、式(3)和式(5)即为以孔压、位移为变量的控制方程,与前人求解思路一致,利用Comsol软件中偏微分方程求解模块可对控制方程进行全耦合求解[24]。

事实上,式(3)反映了力场和流场的相互影响,表示了水-土之间是复杂的流-固全耦合过程,即孔隙流体压力改变导致土体有效应力发生变化[式(5)],进而产生土体压缩等行为。同样地,土体压缩后,在式(3)中耦合项dεv/dt的影响下,流体流动受阻,孔压变化率减缓。两者同时发生,相互影响。

2.2 模型验证

为验证模型有效性,用一维土体经典固结案例进行分析,边界条件与土体参数均与文献[27]中一致。一维土柱高10 m,两侧轴支撑无流动,底部位移约束无流动,上边界为透水位移自由面且施加200 kPa瞬时荷载。土体弹性模量10 MPa,泊松比0.3,渗透系数1×10-7m/d,网格为10个矩形单元。

图3表示施加荷载2000、1.8×104和9.4×104s与文献[27]所用模型对孔压计算结果的对比。与其结论一致,加载初期(2 000 s时刻)孔压瞬时升高,随着时间推移,压力逐渐消散且上层土体孔压消散更快。然而与原模型不同的是,加载初期,本文模型的孔压略低于文献[27]中单向耦合的计算结果。这是由于本文模型在计算过程中考虑了体积变化对压力的影响,加载初期上层土柱体变急剧增加,土体压缩后微观孔隙体积的减小将阻碍流体流动,进而使得孔压略低。但整体而言两者计算结果吻合较好,可验证本文模型的有效性。

图3 本文模型与其他模型计算结果对比

2.3 研究区基准工况分析

根据勘察报告,构建概化地质模型和网格剖分如图4所示,填土厚度、粉质黏土和分化灰岩层厚度分别取3、10和22 m。采用三角形网格,洞周进行局部加密处理,共计1 364个单元格。赋予各层土体物理力学参数后,如表1所示,设置力学边界条件为:左右边界轴承约束,下边界固定约束,上边界为自由边界。考虑到右侧潜水补给,左侧施工抽水,因此流场的边界条件为:右端为定水头边界,左端除长为2 m的抽水井按照指定速率抽水外,其余部分无流动,上层为自由边界。

表1 地层物理参数

图4 概化地质模型示意图

根据实际工况,同时考虑到不同影响因素的对比分析,基准工况按照下述条件设置:抽水速率Vs为水头每天降低0.3 m,即Vs=0.3 m/d;抽水口上端距离地表的高度(抽水位置)H=25 m,溶洞顶端距灰岩顶层距离(顶板厚度)D=1.5 m,溶洞的轴a与轴b相等且a=b=2.5 m,即洞形为圆形。

图5给出了不同时刻水压、位移、有效应力和塑性区的空间分布。各物理场随时间变化规律表明:抽水过程中,渗流力作用于土体,引起流场、力场变化,整个上覆层应力重分布,洞周产生应力集中。由于两侧无位移,导致洞顶竖向位移最大。恒定抽水速率下,抽水边界处的水压随着时间的推移逐渐降低[图5(a)],由于洞内渗流速度更快,从而在洞顶产生孔压变化较大的区域。

图5 不同时刻物理场的空间分布

根据有效应力原理,各地层所受有效应力随着水压的降低而增大。因工况中土体压缩仅由压力变化引发,因此地表竖向位移与水压的空间分布形式一致,均呈现碟状分布[图5(b)],这与目前勘察所得的结果一致。同时,孔压降低导致洞周出现应力集中[图5(c)],随着抽水的继续,导致洞周进一步产生塑性破坏[图5(d)]。

前人多基于有限差分或有限元法建立半耦合模型[19-23]。在这些模型中,首先通过计算流场方程得到压力场,然后通过有效应力原理、连续方程和本构方程求解位移,流场和力场并未同时求解。这就造成了这类模型在一个计算步中并未考虑由于体积变化所导致的对流场的影响。

为比较两种计算方法的差异,图6给出了灰岩顶板位置处不同时刻的水压空间分布。由于地层压实减小了流体的流动空间,因此全耦合模型孔压下降速度比半耦合模型慢,如抽水1 d后全耦计算的孔压比解耦的结算结果最大高10 kPa,即二者有效应力亦相差10 kPa。因此在给定时间内,半耦合模型所计算的有效应力将更大,进而计算得到安全系数也会更大。

图6 全耦和半耦对不同时刻顶层水压分布计算结果对比

为更直观地分析耦合项影响,图7展示了耦合项的空间分布。可以看出,由于岩层的弹性模量明显高于土层,因此其体变相对较小,耦合项的影响更小。而对于土层而言,水位下降过程中体变较大,变形对流场的影响也将更大。计算结果表明土层和岩层在一个时步下的耦合项数值可相差6～10倍,因此在重点分析土体渗流时,这种作用不可忽视。

图7 不同时刻耦合项dεv/dt的空间分布

3 溶洞稳定性影响因素分析

为进一步研究地下水抽取对溶洞稳定和地面沉降的影响,从工程抽水角度而言选取抽水位置H,抽水速率Vs为研究因素,从溶洞自身性质而言选取a=b时(圆形洞体)的洞径d,a≠b时溶洞形状因子a/b(b=2.5恒定)以及顶板厚度D为变量,共计25种工况,具体模拟设计方案见表2。

表2 不同影响因素下的计算工况

3.1 抽水速率和抽水位置对地面沉降的影响分析

随着人为基坑开挖抽水,各物理场之间的关联性和规律已经给出,考虑到地面沉降是岩溶区施工过程中主要监测内容,仅以抽水5 d后的地表垂向位移分析各因素的影响。

地表垂向沉降如图8所示,可以看出不同抽水速率下地面沉降整体呈现碟状分布,抽水速率较小时,5 d后地面沉降最大值出现在溶洞正对地表位置,Vs=0.1 m/d时,最大沉降约为3.2 mm,Vs=0.9 m/d时,洞口处最大沉降已经将近11.3 mm。随着抽水速率增加,抽水井一侧水压下降更为明显,因此井筒附近沉降更大。应注意在计算中未考虑抽水井一侧的施工支护,但在实际施工中,在支护结构的作用下,井筒处位移将被限制,此时洞顶地面沉降值仍将最大。计算结果表明在实际工程中应注意设置合理的抽水速率,尽可能地避免大范围水力压降的出现,进而较少溶洞顶部地面沉降。

图8 不同抽水速率下5 d的地表垂向位移分布

同样地,不同抽水位置的计算结果如图9所示。抽水位置表示抽水井顶部距离地面的距离,抽水速率相同的条件下,可以看出随着抽水位置的下降,抽水侧的地面垂向位移逐渐降低。这是由于随着抽水位置的下降,相同时刻抽水井上方空间同一位置土体的孔隙水压更高,有效应力更低所致。而非抽水侧由于边界条件的补给作用整体孔压下降不明显,因此在非抽水侧的地面位移差异不大。这意味着随着施工初期,抽水井位置一般较浅,其对地面沉降的影响相对更大,应作为监测重点期。

图9 不同抽水位置下5 d时地表垂向位移分布

3.2 洞径、顶板厚度和洞形对地面沉降的影响分析

1.5～3.5 m内5个不同洞径工况的计算结果如图10所示。可以看出洞径较小时洞顶上方的沉降也较小,为-3.3 mm。随着洞周附近水力交换的影响随洞径的增大逐渐增强,洞周压降更为明显,导致溶洞顶部地面沉降更大,洞径为3.5 m时沉降可达8.2 mm。因此实际施工时,应进行充分地前期勘察工作,探明溶洞大小,同时在施工过程中应重点对洞顶地层进行监测。

图10 不同洞径下5 d的地表垂向位移分布

图11展示了不同顶板厚度时的地面沉降计算结果。由于洞体形状和抽水速率完全相同,相同时刻下整个流场在空间的分布相差不大,因此地面沉降的变化并不明显。但从溶洞稳定性的角度分析,顶板越薄其力学性质越差,所对应的塑性区面积将会越大,洞体的稳定性越差。若以顶板塑性区贯通为破坏标准,则顶板坍塌风险与顶板厚度则成正相关。

图11 不同顶板厚度下5 d的地表垂向位移分布

图12 不同洞形下5 d的地表垂向位移分布

使图4中的轴b恒等于2.5 m,通过改变轴a的长短设置5个比例因子进行计算,a/b<1表示上下尖两侧扁的细长椭圆形溶洞,a/b=1表示圆形,a/b>1表示左右尖上下扁的椭圆形溶洞。从计算结果可以看出,随着椭圆形溶洞长短轴的变化,a/b越大,洞体上下的水力交换接触面积越大,压降越明显,地面沉降值也将随之增加,由3 mm增至将近10 mm。

3.3 不同影响因素对溶洞稳定性的影响分析

为分析不同因素对溶洞稳定性的影响,以抽水5 d时最大塑性应变表示溶洞稳定性。塑性形变越大表示溶洞越不稳定。为便于对比不同单位的影响因子ai,将其进行归一化处理,得到归一系数为ai/(amax-amin),其与最大塑性应变的关系如图13所示。可以看出溶洞的稳定性与抽水速率、洞径和a/b值呈负相关而与抽水位置呈正相关,即抽水速率、洞径和a/b值越大,抽水位置相对越靠近地面溶洞稳定性越差。应注意随着顶板厚度的增加,稳定性由差向好但随后又逐渐变差。这是由于一方面顶板越薄其力学性质越差,所对应的塑性应变将越大,另一方面溶洞还受其形成时初始地应力的影响,顶板越厚所对应的溶洞埋深越深,受到的初始地应力随之增大,塑性区也会更大,进而造成顶板厚度与塑性应变之间呈现抛物线形分布。

图13 不同影响因素对溶洞稳定性的影响

4 影响因素敏感性分析

为分析各影响因素对地面沉降和溶洞稳定性的影响程度,为工程防治提供依据,采用灰色关联法进行参数敏感性分析。灰色关联法作为因素敏感性分析的方法,因其克服了传统回归分析和方差分析需大量初始数据的不足,在岩土边坡分析[28]、基坑稳定性分析[29]中应用广泛。

首先以影响岩溶区地面沉降的各因素(抽水速率、抽水位置、洞径、顶板厚度、形状因子a/b)为比较列X(X=[X1X2X3X4X5]T),将其影响的指标,即以洞顶地面垂向位移表示的地面沉降Y和以最大塑性应变表示的溶洞稳定性为参考列Y1,表达式分别为

(6)

(7)

(8)

(9)

式(9)中:ξ∈[0,1],通常取0.5。

将无量纲化处理后的比较列和参考列矩阵代入式(9)计算得到地面沉降和溶洞稳定性的关联系数矩阵分别为

(10)

(11)

图14 敏感性分析雷达图地面沉降溶洞稳定性

可以看出,整体而言,相较于溶洞稳定性,地面沉降对各因子更为敏感。同时可以看出,各个因子对地面沉降和溶洞性定性的影响规律几乎一致,即抽水速率>洞径>形状因子>顶板厚度>抽水位置。但就地面沉降而言,前三者几乎绝对占优。因此为预防沉降,应重点控制抽水速率,关注溶洞大小和形状。而对于溶洞稳定性而言,抽水速率和洞径(溶洞大小)占主导地位。综上,为防止岩溶塌陷及地面沉降的发生,在岩溶区抽水过程中首先应充分控制水头下降速率,使其保持在合适区间,必要时增加回灌井;其次在勘察中应查明溶洞规模和分布,尽量避开较大溶洞区域范围抽水,或抽水前对溶洞进行相应的工程措施处理。

5 结论

岩溶区抽取地下水所导致塌陷的实质是流场和力场相互作用的结果。以深圳市岩溶塌陷易发区工程地质条件为背景,建立流-固耦合分析模型,对不同抽水速率、洞径和洞形等多个关键因素进行分析,基于灰色理论确定影响因素的敏感性大小,得到结论如下。

(1)基于连续介质力学建立了流-固耦合分析模型,采用一维固结案例验证模型有效性后进行基准工况分析。与半耦合模型的计算结果对比表明:由于地层压实减小了流体流动空间,全耦模型孔压下降速度更慢,使得土层体变更小。因此在重点分析松软土体渗流时耦合项的影响不可忽视;对力场的计算结果表明解耦计算结果所得到的安全系数将会偏大。

(2)抽水过程中,渗流力作用于土体,引起流场、力场变化,整个上覆层应力重分布,洞周产生应力集中;地表竖向位移与水压的空间分布形式一致,均呈现碟状分布,这与目前勘察结果一致。所分析的洞径、顶板厚度和抽水速率等5个影响因素中,对压降有明显影响的因素所计算得到的地面沉降和塑性区面积通常更大,即水力压降是形成地面沉降和溶洞坍塌的主要因素。

(3)基于灰色关联理论进行了敏感性分析。其中,地面沉降对抽水速率、洞径和形状因子更为敏感,溶洞稳定性对抽水速率和洞径更为敏感。因此为防止地面沉降的发生及保持溶洞稳定,在岩溶区抽水过程中首先应控制水头下降速率,使其保持在合适区间;其次在前期勘察中对溶洞大小进场查明,尽量避开溶洞较大区域抽水,或抽水前对溶洞进行相应的工程处理。

相关成果可为岩溶区灾害防治提供支持。