地下水开采与地面沉降控制三维全耦合模型研究骆祖江1,曾　峰1,李　颖21.河海大学水文地质研究所,南京　2100982.长春工程学院勘查与测绘学院,长春　130021

摘要:针对第四纪松散沉积层中地下水开采所引起的地面沉降问题,以比奥固结理论为基础,考虑到土体的非线性特征及土体的渗透性随应力状态的动态变化,引入邓肯2张非线性模型和渗透率动态模型,将地下水渗流场和土体应力场耦合起来,并以吴江市浅层地下水开采为例,建立了浅层地下水开采与地面沉降三维全耦合数值模型。在对模型进行校正、识别的基础上,以微承压含水层未来10a内地下水位埋深不低于含水层顶板、地面沉降量不超过50mm为约束条件,预测了吴江市各镇的地下水可采资源量合计为1

203.59×104m3/a。关键词:比奥固结理论;地面沉降;三维全耦合模型;非线性;地下水;开采中图分类号:P641 文献标识码:A 文章编号:167125888(2009)0621080207

收稿日期:2009202219

基金项目:江苏省与国土资源部合作重大项目(200312300009)作者简介:骆祖江(1964—),男,江苏吴江人,教授,博士生导师,主要从事水文地质与工程地质研究,E2mail:luozujiang

@sina.com。

StudyonThree2dimensionalFullCouplingModelofGroundwaterExploitationandLand2SubsidenceControl

LUOZu2jiang1,ZENGFeng1,LIYing21.ResearchInstituteofHydrogeology,HohaiUniversity,Nanjing　210098,China

2.SchoolofProspecting&SurveyingEngineering,ChangchunInstituteofTechnology,Changchun　130021,China

Abstract:Aimedattheland2subsidencecasedbygroundwaterexploitationinQuaternaryloosese2

diments,basedonBiot’sconsolidationtheory,tooktheshallowgroundwaterexploitationinWujiangasthecasestudy,athreedimensionalfullcouplingnumericalmodelbetweenthehydraulicactionoftheseepagefieldandthestressfieldwasestablished,inwhichDuncan2Changnon2linearmodelandthedy2namicmodelsofpermeabilitywereintroducedwithconsideringthenon2linearcharacteristicsandthedy2namicvariationofhydraulicconductivityofsoilswithstress.Basedonidentificationandverificationofthemodel,withtheconstraintconditionthatthegroundwaterlevelwouldnotbelowerthanhalfoftheroofdepthofthefirstconfinedaquiferandtheland2subsidencewoundnotexceed50mm,theallowableexploitablequantityofgroundwaterineverytownwaspredicted,theexploitablequantityofgroundwa2teris1203.59×104m3/ainthewholearea.Keywords:Biotconsolidationtheory;land2subsidence;threedimensionalfullcouplingmodel;non2linear;groundwater;exploitation

0引言地面沉降是一种由多种因素引起的地面高程缓慢降低的地质现象,其中过量开采地下水是造成地面沉降的最主要原因。江苏省吴江市地处长江三角洲,经济发达,城市供水对地下水需求量大且开采历

　第39卷　第6期2009年11月 吉林大学学报(地球科学版)JournalofJilinUniversity(EarthScienceEdition) Vol.39　No.6Nov.2009 史较早,已形成了区域性的地下水位降落漏斗,进而导致了区域性的地面沉降地质灾害问题;因此建立区域性的地下水开采与地面沉降模拟模型,准确预测基于地面沉降控制的地下水可采资源量已成为当务之急。目前,国内外关于地下水渗流与地面沉降问题的计算模型主要有水土分算模型[1]、部分耦合模型[2]和全耦合模型[3]。但大多数模型并未考虑到土体的变形参数以及渗透性随着土体中应力场改变的动态变化。笔者根据比奥固结理论,并考虑土体的非线性特征及土的渗透性随应力状态的动态变化,结合吴江市地下水系统的水文地质机制,建立了吴江市浅层地下水开采与地面沉降的三维全耦合模型;同时计算出地下水开采过程中地下水位和地面沉降的动态变化,为地面沉降的科学防控提供决策依据。1土水耦合模型1.1比奥(Biot)理论假定饱和土体中,土体骨架是线弹性、变形微小、地下水渗流符合达西定律且不可压缩,取z坐标向上为正,应力以压为正,则三维比奥固结方程如下[4]:-G󰂈2wxG1-2v・55x(5wx5x+5wy5y+5wz5z)+5u5x=0,-G󰂈2wyG1-2v・55y(5wx5x+5wy5y+5wz5z)+5u5y=0,-G󰂈2wzG1-2v・55z(5wx5x+5wy5y+5wz5z)+5u5z=-γ。(1)55t(5wx5x+5wy5y+wz5z)1γw[55x(kx5u5x)+55y(ky5u5y)+55z(kz(5u5z+γw))]=0。(2)式中:G为剪切模量(Pa);v为泊松比;wx、wy、wz分别为x,y,z方向上的位移分量(m);u为孔隙水压力(Pa);kx,ky,kz分别为x、y、z方向上的渗透系数(m/d);γ为土的重度(N/m3);γw为水重度(N/m3)。利用Galerkin加权余量法离散方程,考虑到土体的非线性特性,取Δt时间内的位移增量来代替位移,将式(1)、(2)离散成增量形式[526]:KCCTΔtH+BΔrΔu=ΔfΔtQ-ΔtHu。(3)式中:K,H分别是固体刚度矩阵和流体“刚度”矩阵;C是耦合矩阵;Q是流量矩阵;B是关于自由面的积分矩阵;Δr为结点位移增量(m);Δf为等效节点力增量(Pa)。上述方程结合一定的定解条件即可运用伽辽金有限单元法进行求解。1.2方程的定解条件1.2.1初始条件地应力初始条件:采用土体的自重应力估算土体的初始应力,

σz0=γz,

σx0=σy0=K0γ

z。

(4)

式中:

σ

x0,σy0为土体的初始水平向应力(Pa);σz0为土

体的初始垂向应力;z为计算点深度(m);K0为静止

侧压力系数,K0=

1-sinφ′砂土

0.95-sinφ′粘土

,φ′为有效

内摩擦角。位移初始条件:w(x,y,z,t)|t=0=0。(5)

孔隙水压力初始条件:

u(x,y,z,t)|t=0=u0(x,y,z)。(6)式中:u

0(x,y,z)为区域内已知初始孔隙水压力。

1.2.2边界条件孔隙水压力边界条件Γ1:

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

。(7)

式中:us为水头边界Γ1上的已知孔隙水压力。

流量边界条件Γ2:

K5H5n|Γ2=q󰁫L

。(8)

式中:q󰁫L为边界Γ2上的已知单位面积流量。

自由面边界条件Γ4:

u=0;q=-μ5u5tcosθ。(9)

式中:μ为土体给水度;θ为自由面外法线方向与垂线的交角;q为通过自由面边界Γ4的单位面积流

量。1.2.3自由面边界条件的处理在地下水开采的计算过程中,自由面的位置是未知的,需要通过迭代才能求出[728]。应用改进的截止负压法[9210]将式(2)拓展到整个区域,使待定的自由面边界固化,迭代求解。并引入罚函数M(P),得到式(2)的弱表达式:

∫Ω(󰂈sk󰂈(pγ)+M(p)󰂈sk󰂈z+

s5εv5t)dΩ-∫Γ4sqdΓ=0。(10)式中:

1801　第6期 骆祖江,曾　峰,李　颖:地下水开采与地面沉降控制三维全耦合模型研究M(p)=0,pG

(1-θ)H

。(14)

式中:当θ=0时,为显示差分;θ=1时为隐式差分;

θ=1/2时为中心差分。其中隐式差分从理论上来

说,时间步长的选择是没有任何限制条件的,也就是说隐式差分是无条件稳定的,这也是本文基于隐式差分编制相应有限元方程的主要原因。

2应力场中的本构关系虽然比奥固结理论在方程推导时假定土体骨架是线弹性的,但采用了有限单元法,就可以将其推广应用于非线性问题[11]。本文采用邓肯2张非线性E2v模型来模拟土体的非线性本构关系。其中E、v不再视为常数,而是随着应力状态而改变。切线泊松比vt与切线模量Et表达式如下:

vt=G-Flg[σ3pa][1-D(σ1-σ3)kpa(σ3pa)n[1-Rf(1-sinφ)(σ1-σ3)2ccosφ+2σ3sinφ]]2,

(15)Et=[1-Rt(1-sinφ)(σ1-σ3)2ccosφ+2σ3sinφ]2kpa[σ3pa]n。

(16)式中:pa为大气压(Pa);参数c,φ,k、n、Rf、G、F、D

可由三轴固结排水剪切实验得到。

3孔隙度n和渗透率k的动态模型流固耦合问题实际上是孔隙应力的消散引起土体骨架的变形、孔隙率的变化,从而影响土体的渗透性,宏观上表现为土体的固结变形。在比奥固结的假定条件下,根据孔隙度的定义和渗流力学Kozeny2Carman方程推得孔隙度n和渗透率k的动态表达式[12213]:

n=n0+εv1+εv,(17)

k=k01+εv[1+εvn0]3。(18)

式中:n0为初始孔隙度(无量纲);k0为初始渗透率;

εv=εx+εy+εz为体应变。

根据以上的推导和论述,本文提出的n动态模型是有效的而且是较为严格的。对于此公式,体应变εv可以是弹性,可以是塑性。通常情况下,它是

弹塑性变形或者粘弹性变形[14]。

2801 吉林大学学报(地球科学版) 第39卷