土体自重固结压缩对地面沉降影响研究

，，

(1.江苏南京地质工程勘察院，江苏 南京 210041; 2.江苏省地质矿产勘查局，江苏 南京 210018)

0 引 言

地下水开采引起地下水渗流场变化和土体自重固结压缩引起的应变场变化相互叠加，致使区域地面沉降愈来愈严重。国内外学者对区域地面沉降的研究起步很早，但更多的是对基于抽水引发的地面沉降影响研究(陈希崇等，2001；骆祖江等，2008；Luo et al,2011)，对土体自重固结压缩引起的地面沉降研究尚不多见。

为了准确评价土体自重固结压缩对地面沉降的影响，以比奥固结理论为基础，结合土体非线性流变理论，将比奥固结理论中的本构关系推广到黏弹塑性，同时考虑水力学参数及土力学参数随有效应力的动态变化过程。以江苏南通为例，建立南通市土体自重固结压缩、地下水开采与地面沉降三维全耦合模型，模拟预测了南通在地下水停止开采、仅在土体自重固结压缩影响下的地面沉降发生、发展趋势，为南通市地面沉降科学防控提供决策依据。

1 研究区概况

研究区域包括南通市区(崇川、港闸)和通州区。全区第四纪松散沉积物发育，分布范围广，最大厚度达300 m左右，含水层以第四纪松散层孔隙水为主。垂向上从上往下分别为潜水含水层、第I承压含水层、第II承压含水层、第III承压含水层以及各含水层间的黏性土弱含水层。其中，潜水含水层由第四系全新统地层组成，岩性以粉砂、亚砂土为主，潜水位埋深1～3 m，厚度20～30 m；第I承压含水层由上更新统冲积、冲海积松散地层组成，岩性以中细砂为主，顶板埋深50～60 m，承压水头埋深1～3 m；第Ⅱ承压含水层由中更新统河流、河口相沉积的松散层组成，岩性以粉细砂、中粗砂及砂砾为主，顶板埋深约140 m，承压水头埋深3～5 m；第III承压含水层由下更新统长江古河道沉积的松散层组成，岩性以中细砂、中粗砂为主，局部有含砾卵石，顶板埋深187～270 m，承压水头埋深30～40 m。

2 地质概念模型

研究区平面上以行政区域为界，面积为1 586 km2，垂向深度在230～350 m。垂向上，从上往下将潜水含水层、第I承压含水层、第II承压含水层、第III承压含水层及各含水层之间的黏性土弱含水层剖分成独立的7个层位参与计算，各含水层在水力学上均概化为非均质各向异性，在土力学上均概化为非均质地层。研究区南部是长江边界，其从上往下切割潜水含水层、第I承压含水层及之间的黏性土弱含水层，故在水力学上可概化上述含水层的长江边界为定水头边界，其余在水力学上均概化为第二类流量边界，系统的顶部一方面接受大气降水及农业灌溉入渗的补给，是补给边界；另一方面，地下水又通过其蒸发，是排泄边界，系统的底部为隔水边界。研究区所有侧向边界和底部边界在土力学上均概化为零位移边界。研究区地下水开采是唯一的源汇项，将其概化为大井处理。土体在自重应力作用下发生固结压缩。

3 数学模型的建立与求解

3.1 比奥固结理论

饱和土体中假定土骨架变形为线弹性、微小变形、渗流符合达西定律、水不可压缩或微压缩的三维比奥固结方程如下(钱家欢等，1996)。

(1)

(2)

式(1)、式(2)中，G为剪切模量；ν为泊松比；wx、wy、wz分别为x、y、z方向上的位移分量；u为孔隙水压力；kx、ky、kz分别为主轴x、y、z方向上的渗透系数；γ为土的重度；γw为水的重度。

3.2 土体本构模型

土的本构关系是土的力学特性即应力-应变-强度-时间等关系的数学表达式。对于考虑流变特性的土体来说，其变形特征主要表现为变形的时间与应力水平有关，所显示的是具有弹性、塑性和黏滞性的黏弹塑性体，若将此类土体的总应变增量dε分为弹塑性应变增量dεep、黏弹性应变增量dεve、黏塑性应变增量dεvp，则具有流变特性的土体中任意点在任意时刻的应变增量为(维亚诺夫，1987；钱家欢等，1996；陈晓平等，2001；姚仰平等，2009)：

dε=dεep+dεve+dεvp

(3)

式(3)中各部分的应变增量可以由下述方法确定。

3.2.1 弹塑性应变增量 由土体的弹塑性本构模型可得：

dεep=[C]dσ′

(4)

3.2.2 黏弹性应变增量 由Kelvin流变模型E1ε+Keε=σ可得应力不变时复杂应力状态下的黏弹性应变增量为：

(5)

式(5)中，ηe=E1/Ke，E1为Kelvin体黏弹性模量；Ke为黏滞系数；[A]为应力矩阵，

3.2.3 黏塑性应变增量 利用黏塑性法确定黏塑性应变增量，在该种方法中允许材料在有限“期间”内超越破坏准则(以破坏准则函数F的值>0来表示)。在讨论土体的黏塑性应变而非塑性应变时，应变的变化率与超越量有关，即有以下关系式：

(6)

式(6)中，Qs为塑性势函数；F为破坏准则函数，对于摩尔-库伦材料来说，

其中：φ为摩擦角，c为黏聚力。σ=(σx,σy,σz,τxy,τyz,τzx)T，

E,ν分别为弹性模量和泊松比。

如果将黏塑性应变率与一个伪时间步相乘，就可以得到累加到下一个时间步的黏塑性应变增量，于是有：

(δεvp)i=Δt(εvp)i

(7)

(Δεvp)i=(Δεvp)i-1+(δεvp)i

(8)

数值计算绝对稳定的时间步与假定的破坏准则有关。对于摩尔-库伦材料有：

(9)

塑性势函数对应力的偏导数可以表示为：

(10)

(11)

将式(4)、(5)、(11)代入式(3)，可求得弹塑-黏弹-黏塑性体的应变增量为(殷宗泽等，2006；缪林昌，2007；姚仰平等，2009)：

(12)

式(12)即为土体的黏弹塑性本构方程。

3.3 比奥固结有限元方程

利用伽辽金加权余量法离散方程，考虑到土体的非线性特性，取Δt时间内的位移增量来代替位移，将式(1)、(2)离散成增量形式(李医民等，2004)：

(13)

因为渗流取决于孔隙压力全量的分布，而不是取决于时间内孔隙压力增量，所以孔压要用全量的形式表示，记时刻tn和tn+1时单元节点i的孔压全量分别为ui(n)和ui(n+1)，且Δui=ui(n+1)-ui(n)，则式(4)可变换为：

(14)

式(14)即为三维比奥固结有限元方程。

3.4 参数的动态变化模型

3.4.1 孔隙度与渗透系数的非线性 流固耦合问题实际上是孔隙应力的消散引起土体骨架的变形，孔隙系数的变化，从而影响土体的渗透性，宏观上表现为土体的固结变形。在比奥固结的假定条件下，根据孔隙度的相关定义和渗流力学Kozeny-Carman方程推得孔隙度n和渗透系数k的动态表达式(田杰等，2005)：

(15)

3.4.2 土体参数的非线性 采用邓肯-张非线性模型，将土体的本构关系推广到非线性，则本构关系{Δσ}=[D]{Δε}中矩阵[D]中的弹性常数E、ν不再视为常量，而是随着应力状态改变而改变，其切线弹性模量和切线泊松比的表达式如下(罗刚等，2004)：

(16)

(17)

式(16)、式(17)中，Rf为破坏比；c为黏聚力；φ为内摩擦角；σ1为第一主应力；σ3为第三主应力；n为弹性模量与固结压力曲线的斜率(lgα)；G为土体常规三轴压缩实验结果所绘曲线截距；F=0.04，D=3为土体实验参数；pa为大气压强。

3.5 定解条件

3.5.1 初始条件 初始条件是指在初始时刻(t=0)所研究对象各个求解变量的空间分布情况。初始条件越接近真实值，迭代收敛所用的时间越短，模型计算结果越符合实际情况。

(1) 地应力初始条件：

采用土体的自重应力估算土体的初始应力：

(18)

(2) 位移初始条件：

w(x,y,z,t)|t=0=0

(19)

(3) 孔隙水压力初始条件：

u(x,y,z,t)|t=0=u0(x,y,z)

(20)

式(20)中，u0(x,y,z)为研究区域内已知初始孔隙水压力。

3.5.2 边界条件 (1) 孔隙水压力边界条件Γ1：

u(x,y,z)|Γ1=us

(21)

式(21)中，us为水头边界Γ1上的已知孔隙水压力。

(2) 流量边界条件Γ2：

(22)

(3) 自由面边界条件Γ3：

(23)

式(23)中，μ为土体给水度；θ为自由面外法线方向与垂线的交角；q为通过自由面边界Γ3的单位面积流量；Z为自由面所在的高程。

(4) 位移边界条件Γ4：

(24)

式(24)中，wx,wy,wz为位移边界Γ4上3个方向的已知位移。

上述比奥固结有限元数学模型即可运用Fortran语言编制相应的有限元程序进行求解(钱家欢等，1996；毛昶熙等，1999)。

4 模型识别与验证

对整个研究区域用八节点六面体单元进行离散化，在平面上剖分成6 489个矩形网格单元，剖面上按地层岩性(含水沙层和黏性土相对隔水层)剖分成7层。单元总数为45 423个，节点总数为13 514个。研究区三维网格剖分详见图1。

图1 三维网格剖分图

利用前述开发的有限元计算机模型，选取2009年12月31日—2010年12月31日作为模型识别、验证的时段，每个月作为1个应力期，共分12个应力期，每个应力期为1个时间步长。各含水层的初始孔隙水压力由初始水位值换算得出，初始位移值及边界上的位移均为0。对观测井的水位进行拟合、反演，获得了各含水层的水文地质参数，整个模型参数分区共58个，以第III承压含水层组为例，参数分区如图2所示，其参数分区参数值见表1。

图2 第III承压含水层组参数分区

图3 地下水位拟合图(南通市烟滤嘴厂)

图4 地下水位拟合图(通州二甲水厂)

图5 2010年地面沉降量计算值与观测值对比(mm)

图3、图4、图5举例说明了模型的地下水位值和地面沉降值拟合精度较高，从识别的结果看，各参数分区参数值的级别大小均符合常规，识别过程中所涉及的物理量都是在实际调查的基础上确定的，因此模型的识别具有一定的精度和可信度。

5 模型预测

利用上述经过识别和验证的土体自重固结压缩、地下水开采与土体变形三维全耦合数值模型，模拟预测了研究区在地下水停采、仅考虑土体自重固结压缩的影响下，自2010年12月31日—2025年12月31日各含水层地下水流场变化特征和固结压缩沉降的发展趋势。模型计算得出：南通市各层最大固结压缩量均位于通州区，累计由土体自重固结压缩导致的最大地面沉降量为2.42 mm，最大地面沉降速率为0.16 mm/a。图6、图7分别举例说明了2025年12月31日第III承压含水层的地下水预测流场变化和预测固结压缩量分布特征，图8举例说明了2025年12月31日的预测地面沉降量分布特征。

从计算结果可以看出，南通市由土体自重固结压缩引发的地面沉降，在平面上远离市区的地方地面沉降值逐渐增大，符合远离市区的土体固结程度偏低，弹性模量较低的规律。

表1 第Ⅲ承压含水层参数分区表

图6 2025年12月31日第III承压含水层预测流场(m)

图7 2025年12月31日第III承压含水层压缩量预测等值线图(mm)

图8 2025年12月31日预测地面沉降量等值线图(mm)

6 结 论

(1) 以比奥固结理论为基础，并结合土体流变理论，将土体水力学参数以及土力学参数随有效应力的动态变化引入到模型中，建立的土体自重固结压缩、地下水开采与地面沉降三维全耦合模型，可准确刻画土体自重应力、地下水开采与地面沉降之间的相互作用机制，提高模型计算的置信度。

(2) 南通市在地下水停采、仅在土体自重固结压缩影响下，自2010年12月31日—2025年12月31日的最大地面沉降量为2.42 mm，最大地面沉降速率为0.16 mm/a，土体自重固结压缩对南通市地面沉降的影响极为有限。

(3) 南通市由土体自重固结压缩引发的地面沉降，在平面上远离市区的地方，地面沉降值逐渐增大，符合土体固结程度偏低，弹性模量较低的规律。

陈崇希，裴顺平.2001.地下水开采-地面沉降模型研究[J].水文地质工程地质，28(2)：5-8.

陈晓平,白世伟.2001.软黏土地基黏弹塑性比奥固结的数值分析[J].岩土工程学报,23(4):481-484.

李医民,周凤燕.2004.一类三维偏微分方程边值问题的解法[J].江苏大学学报:自然科学版,25(4):328-331.

罗刚,张建民.2004.邓肯-张模型和沈珠江双屈服面模型的改进[J].岩土力学,25(6):887-890.

骆祖江，刘金宝，李朗.2008.第四纪松散沉积层地下水疏降与地面沉降三维全耦合数值模拟[J].岩土工程学报，30(2):193-198.

毛昶熙,段祥宝,李祖贻,等.1999.渗流数值计算与程序应用[M].南京:河海大学出版社.

缪林昌.2007.非饱和土的本构模型研究[J].岩土力学,28(5):855-860.

钱家欢,殷宗泽.1996.土工原理与计算[M].北京:水利水电出版社.

SMITHI M, GRIFFITHS D V. 2003.有限元方法编程[M].王菘,周坚鑫,王来,等译.北京:电子工业出版社.

田杰,刘先贵,尚根华.2005.基于流固耦合理论的套损力学机理分析[J].水动力学研究与进展:A辑,20(2):221-225.

维亚诺夫C C.1987.土力学的流变学原理[M].杜余培，译.北京:科学出版社.

殷宗泽,周建,赵仲辉,等.2006.非饱和土本构关系及变形计算[J].岩土工程学报,28(2):137-146.

姚仰平,侯伟.2009.土的基本力学特性及其弹塑性描述[J].岩土力学,30(10):2881-2902.

LUO ZUJIANG, ZENG FENG. 2011. Finite element numerical simulation of land subsidence and groundwater exploitation based on visco-elastic-plastic biot’s consolidation theory [J]. Journal of Hydrodynamics, 23(5): 615-624.