[工程勘察]承压含水层非稳定流有限层分析-修改稿

承压含水层非稳定流有限层分析

王旭东 宰金珉

(1. 南京工业大学 岩土工程研究所,江苏 南京 210009;2. 浙江省建筑科学设计研究院,浙江 杭州 310012)

摘 要:推导了以伽辽金法结合傅立叶函数为基础的各向异性承压含水层非稳定流有限层方程,编制了相应的计算程序。通过算例的有限层解与解析解的对比分析,在验证有限层法正确性的基础上,深入探讨了计算区域、计算项数等因素对计算结果的影响,得到了计算区域尺寸影响有限层解的收敛速度等结论,最后给出了一定精度条诸宏博

1,211

件下计算项数选用的近似估算公式。

关键词:有限层法;承压含水层;非稳定流;降深

0 引 言

有限层法是一种对某一方向进行数值离散,而在其余两方向采用连续函数的半数值半解析方法,这种方法能有效地将三维问题简化为一维问题求解,从而大大减少单元数,节省计算工作量和内存的需要量。基于这些特点,Cheung在1968年首先提出了有限条法,并成功地将有限条法应用于结构

工程问题的分析[1]

。Small和Booker等运用了有限层法对层状弹性地基和土的固结问题做了较为系统

的研究[2-4]

。A.M. oskoorouchi和 B. Novrouzian将有限条方法推广应用于岩土工程分析中,对基坑开

挖、隧道和土坝等问题进行了研究分析[5]

。Stanley,

Myron等将有限层理论应用于地下水分析[6]

。张问清、赵锡宏、宰金珉等采用有限层方法探讨了任意

力系作用下的层状弹性半空间问题[7]

。宰金珉、王旭东等将有限层法应用于土与结构物的共同作用分析[8-9]。

地下水计算方法主要有解析法和数值法。解析法可以将描述地下水运动的各种物理量反映在一个表达式中,使用方便,但自然界地下水运动状态复杂,往往与解析解的假设条件相差甚远,其应用范围也受到一定的限制。对于一些复杂的地下水问题,就必须借助于数值模拟方法求其近似解。目前,地下水数值分析中主要有有限差分法、有限元法、边

界元法等数值方法[10-12]

。数值法能较为方便地处理解析解难以解决的复杂边界等问题,因此,得到广泛的应用。但是,数值法计算工作量大,尤其是在考虑三维模型的计算分析时,占用内存多。国内外学者不断尝试寻找有效的数值方法用于地下水模拟计算,薛禹群,叶淑君,谢春红等将多尺度有限元法应用于地下水模拟,节省了计算量提高了计算精度[13]。

本文吸取有限层法的优点,以伽辽金法结合傅立叶函数为基础,推导了承压含水层非稳定流有限

层方程,编制了相应的计算程序,在验证方法的正确性基础上,探讨了计算区域、计算项数等因素对计算结果的影响。

1 承压含水层非稳定流模型的数学描述

假设均质、各向异性、等厚的有限区域承压含水层,在区域中心(x0 , y0)处有一个以恒定流量Q抽水的完整井,初始水位水平,定水头边界条件,如图1所示。根据质量守恒定律和达西定律,承压含水层非稳定流的数学模型用降深可表示为:

KssssxxxyKyyzKzzR(x,y,z,t)S

s

t

t>0, a≤x≤a,b≤y≤b

(1)



s(x,y,z,t)0 t0, a≤x≤a,b≤y≤bs(x,y,z,t)x-a0 t>0, b≤y≤bxa



s(x,y,z,t)y-b0 t>0, a≤x≤ayb式中,s(x,y,z,t)为降深,Kx、Ky、Kz分别为x,y,z

方向的渗透系数,Ss为单位贮水系数,R(x,y,z,t)为源汇项,补给源取正值,排泄汇取负值。

(1)式可用微分算子表示为:

LsR0 (2) 式中,Ls

s



s



s

s

x

Kxx

y

Kyy

z

Kzz

Sst,L为微

分算子。

───────

收稿日期:

基金项目:国家自然科学基金项目(50278042);江苏省高校自然科学研究计划项目(03KJB560044) 作者简介:诸宏博(1981-),男(汉族),浙江杭州人,硕士研究生,主要研究方向为环境岩土工程及数值分析。E-mail: [email protected]

图1 承压含水层物理模型与有限层层元 Fig.1 Configuration of confined aquifer physical model

and finite layer models

2 有限层格式推导

有限层计算简图如图1所示。x和y方向的计算尺寸分别取为2a和2b,沿z轴方向将承压含水层离散成L个层元。

2.1 试探函数

设降深在x和y方向分别用解析函数表示,而z方向采用有限层元离散,层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数: L1M

N

s

x,y,z,tΦmnjtAmn(x,y)Nj(z) (3) j1m1n1式中,~sx,y,z,t为降深试探函数;Φmnj (t)为待求系数;Amn(x,y)sinkmxasinknyb为满足定水头边界的已知函数,kmm/2a,knn/2b;M,N为级数项数;L为有限层层元数。

Ni(z)为线性插值基函数,其形式如下:

zzi1/zizi1zi1≤zNz≤zizz (4) ii1/zizi1zi≤z≤z

i1

0

z≤zi1或 z≥zi1

2.2 伽辽金方程

根据加权余量法原理,可得有限层法伽辽金方程:

Ls~

D

Rmni

dxdydz0,其中 i1,2,,L1 (5)

式中,mni为权函数,mnix,y,zNizAmn(x,y)。

将含水层离散为L个层元,则(5)式可改写为:

L

Ls~

Rmni

dxdydz0,其中i1,2,,L1 (6)

e1

De

将(3)式代入(6)式得:

LL1MN

e

e

e

LΦmnjtAmnNj

Re1j1m1n1

D

Ne

iAmndxdydz0

(7)

其中 i1,2,,L12.3 层元有限层方程

取一典型有限层层元l,其节面分别为l,l+1,层厚c,如图2所示。

图2 典型有限层层元

Fig.2 Configuration of finite layer element models

由Amn(x, y)和Am′n′(x, y) 的正交性,可得:

a

 0 asinkm

xsinkmxdxmm

a mm

(8)

b

 bsinkn

ysink0 nn

nydyb nn

(9)

故,当 mm和nn时,由(7)式得:

l1M

N

ΦZl1mnjKNeNedza

jlm1n1Zlxijab

2Amn

bx2Amndxdy

b

2 Φl1

eea

mnj

Z(10)

ZKl

yNNdzAmn

ij

aby2Amndxdy

ΦNej

b

mnjKzzNec

i0

a

a

b

A2mndxdy ΦZel1

KNNe

ij

a

b

mnjZl

z

zz

dz

a

b

A2mndxdy

Φmnj

Zl1

t

eZSA2l

sNiNejdz

a

a

b

b

mndxdy

Zl1a

b

Zl

a

b

RlNejAmndxdydz0, 其中 il,l1

由基函数的性质可得:

Zl1

ZNeNecl3 ijij

dzl

cij

(11) l 由(11)式可分别得到m、n阶的层元渗透矩阵系

数ge

、贮水矩阵系数beemnijmnij和水量矩阵系数qmnij:

g

e



Zl1

eea

b

mnij

ZKNiNj

dzl

xa2Amnb

x2Amndxdy



Zl1

Keea

b

2ZyNiNjdz

Amnb

y

2

Alamndxdy (12a)

e

KN

i

Zab

zz

Nj0

a

b

A2mndxdy 

Zl1

KNeNe

ija

b

Zl

zzz

dz

a

b

A2mndxdy

l1

b

eZl1

eemnij

sNiNjdz (12b)

il

ZSl

a

a

b

b

A2

mndxdyqe1

a

b

emnj

ZlZRlNjAmndxdydz

(12c) l

a

b

式中,i=l,l+1,j=l,l+1。

由(12a)式、(12b)式积分可得:

abK22xlkmKylkncl

Kzl ij

ge



3cl

mnij



ab

Kk2

xlmKylk2

n

c (13a)

lKz6lc ijl

abclSsij be

mnij

3abc (13b)

lSs

6 ij对于位于区域中心处的完整抽水井,由(12c)

式积分可得: l1b

qemnj

zzNyl

jdz

a

a

b

x0,0x,yAmndxdy

(13c)

cl

2

x0,y0sinkmasinknb, 其中jl,l1式中,x0,y0为抽水率x0,y0Qx0,y0/cl。

层元l的有限层方程为:

gel,lgebel,l1ll,lbel,lgeqel

el1,l

ge1

dll1,l1mnl1mnbel1,lbl1,l1mndtl1mn

qe0(14) l1mn

用矩阵可简化表示为:

Ge

e

de

mnΦmnBmndt

ΦmnQmn0 (15) 2.4 整体有限层方程

由层元有限层方程组装可得整体有限层方程:

GdmnΦmnBmndt

ΦmnQmn0 (16)

L

L

L

式中,GemnG,Bmne1

Be,Qmne1

Qe

mn

mn

mn,

e1

ΦmnT

12

jLL1mn。

整体矩阵Gmn、Bmn由层元矩阵叠加而成,叠加方法如图3所示。

图3 [G]mn和[B]mn矩阵叠加 Fig.3 Matrix assembly for [G]mn and [B]mn

3 有限层方程求解

对于非稳定流来说,降深是时间t的函数,则

Φmnj(t)也是时间t的函数,

采用差分法对时间进行离散得:

dΦk1k

mn

ΦmnΦmndt

t

(17)

假设在时间k和k+1之间,水头呈线性变化,

则有:

Φk

mn

Φk1k

mn1Φmn

(18) 式中, Φk1

mn和Φkmn分别为k时刻和k+1时刻的值,t为时间步长,为加权因子,取值范围在0~1之间。

将(18)式代入(16)式可得有限层方程的差分格式:

(G1k1

mn

t

Bmn)Φmn

(19) (1kkt

Bmn1Gmn)ΦmnQmn

在(19)式中,当0时,即为显式差分格式;1时,为隐式差分格式;0.5时,为Crank-Nicolson差分格式,其中,隐式差分格式和Crank-Nicolson差分格式无条件稳定。

(19)式有限层方程的差分格式为L+1阶线性代数方程组,求解方程组可得到未知系数Φmnj (t),并由(3)式计算k+1时刻的降深。

数值计算时初始降深为0,故Φ0

mn0,计算程序流程图如图4所示,并编制了相应的计算程序,采用追赶法

[14]

求解方程,以节省计算时间。

图4 有限层计算程序流程图

Fig.4 Program flow graph of finite layer method

4 数值算例及分析

4.1 有限承压含水层中抽水井

算例:承压含水层厚c=100m,为了便于对比分析,渗透系数Kx=Ky=Kz=5m/d,单位贮水系数Ss=0.2×10-4m-1;计算区域取为2a×2b=5000m×5000m。(1).在计算区域中心布置一口单井,单井抽水量Q =2000 m3/d;(2).在计算区域中200m×200m范围内布置8口群井,单井抽水量Q =800m3/d。分别采用有限层法模拟计算抽水引起的降深。

为了验证承压含水层非稳定流有限层分析方法的正确性,引入无限承压含水层非稳定流完整井解

析解[12]

(泰斯公式)进行对比分析。

图5和图6分别给出了单井和群井降深随距离的分布曲线。由图可知,无论是单井还是群井,有限层解与解析解都有较好的吻合程度,特别是在抽水较短时间内,降深不受边界影响,此时两种方法的计算结果尤为一致,表明了本文方法的正确性。 4.2 计算区域的影响

由于本算例为有限含水层问题,所以有限层分析中,随着抽水时间的延长,降深逐渐受到边界的影响,尤其是边界附近的降深值,在图6中,抽水10天时,边界附近的有限层解与解析解之间已存在

一定的差异。

图7 给出了不同计算区域情况下距抽水井100m处有限层和解析解的标准降深曲线,两者对比表明,在抽水初期,抽水引起水位下降的区域较小,基本不受边界的影响,不同计算区域的有限层解与解析解都有较好的一致性,且具有一致性的区域(时间) 随着计算区域的增大而增加,表明有限区域内的有限层解随着计算区域的增大逐渐逼近无限区域的解析解(泰斯公式)。

m

/ 深降距抽水井的距离 / m

图5 单井降深分布曲线 Fig.5 The drawdown curve of single well

距群井中心的距离/m

m

/深降

图6 群井降深分布曲线

Fig. 6 The drawdown curve of multiple wells

cKsQ

41/u = 4Kt/r2

SS

图7 标准化降深曲线(距井100m)

Fig.7 Normalized drawdown curves (100 m from well location)

4.3 计算项数的影响

为了避免边界对计算结果的影响,选取t =1d时的有限层和解析解的计算结果,以讨论计算项数对有限层计算结果的影响,由图7可知,t =1d时可忽略计算区域对有限层计算结果的影响。从图8给出的不同计算区域时有限层解的收敛性过程可知,计算结果的精度随计算项数的增加而呈周期性变化,并逐渐提高;随着计算区域的增大,计算结果的收敛速度减慢。根据有限层解与解析解的对比分析,当计算项数M(或N)取近似估算公式a/25或

b/25

解层解

限析有解M(或N)

图8 M、N与降深值收敛的关系曲线 Fig. 8 Convergence of drawdown value with M and N

5 结 论

(1) 以伽辽金法结合傅立叶函数为基础,推导了承压含水层非稳定流有限层方程,并编制了相应的计算程序。

(2) 算例分析验证了有限层方法的正确性;同时,表明了该方法对地下水非稳定流模拟具有较好的适用性,为有限区域承压含水层非稳定流计算提供一条新的有效途径。

(3) 计算区域和计算项数对计算结果影响的分析表明,计算结果的收敛速度随计算区域的增大而减慢,当计算项数M(或N)取a/25或b/25时能使有限层解的精度保持在±3%以内。 参考文献:

[1] Cheung Y. K., Finite strip method in structural mechanics [M]. Pergamon, New York,1976

[2] Small J C, Booker J R. Finite layer analysis of layered elastic materials using a flexibility approach [J]. Int. J. Numer. Methods Eng.,1984, 20(6):1025-1037.

[3]Booker J.R., Small J.C.,Finite layer analysis of consolidation I [J]. International Journal for Numerical Analysis Methods Geomechanics.1982a,6(2),151-171.

[4]Booker J.R., Small J.C.,Finite layer analysis of consolidation II [J]. International Journal for Numerical Analysis Methods Geomechanics.1982b,6(2),173-194.

[5] A.M. oskoorouchi, B. Novrouzian. Zoned finite strip method and its application in geomechanics [J]. Computers and Geotechnics, 1991,11(2):265-294.

[6] Stanley Smith S, Myron B. Allen, Jay Puckett. The finite layer method for groundwater flow models [J].Water Resources Research. 1992, 28(6): 1715-1722.

[7] 张问清, 赵锡宏, 宰金珉. 任意力系作用下的层状弹性半空间的有限层分析方法[J].岩土工程学报,1981, 3(2):27-42. (ZHANG Wen-qing, ZHAO Xi-hong, ZAI Jin-min. Finite layer analysis for layered and elastic half-space under arbitrary force systems [J]. Chinese Journal of Geotechnical Engineering, 1981, 3(2):27-42. (in Chinese))

[8] 宰金珉, 宰金璋. 高层建筑基础分析与设计[M].北京:中国建筑工业出版社,1993. (ZAI Jin-min, ZAI Jin-zhang. Foundation analysis and design of high-rise building[M]. Beijing: China Architecture and Building Press, 1993. (in Chinese))

[9] 王旭东,魏道垛,宰金珉.群桩-土-承台结果共同作

用的数值分析[J].岩土工程学报,1996,18(4),27-33. (WANG Xu-dong ,WEI Dao-duo, ZAI Jin-min. Numerical Analysis of Pile Groups-Soil-Pile Cap Interaction[J].Chinese Journal of Geotechnical Engineering, 1996, 18(4) :27-23. (in Chinese)) [10] 朱学愚,钱孝星,刘新仁.地下水资源评价[M],南京:南京大学出版社,1987 (ZHU Xue-yu, QIAN Xiao-xing, LIU Xin-ren. Groundwater Resource Evaluation[M]. Nanjing: Nanjing University Press, 1987. (in Chinese))

[11] 陈文华,王旭东.二维非均质体中非稳定流问题的边界单元法[J].岩土工程学报,1991,13(3),65-71. (CHEN Wen-hua, WANG Xu-dong . 2D Analysis of Unsteady groundwater in Non- homogeneous

aquifer by Boundary Element

method[J].Chinese Journal of Geotechnical Engineering,

1991,13(3),65-71. (in Chinese))

[12] 孙讷正. 地下水流的数学模型和数值方法[M].北京:地质出版社,1981. (SUN Ne-zheng. Mathematical Model and Numerical Method of Groundwater Flow[M]. Beijing: Geological Publishing House, 1981. (in Chinese))

[13] 薛禹群,叶淑君,谢春红,张云. 多尺度有限元法在地下水模拟中的应用[J].水利学报,2004(7),7-13. (XUE Yu-qun, YE Shu-jun, XIE Chun-hong, ZHANG Yun. Application of multi scale finite element method to simulation of groundwater flow[J]. Journal of Hydraulic Engineering, 2004(7),7-13. (in Chinese))

[14] 徐士良. FORTRAN常用算法程序集(第二版)[M].北京:清华大学出版社,1995. (XU Shi-liang. FORTRAN Common Algorithms Program (Edition Two) [M]. Beijing: Tsinghua University Press, 1995. (in Chinese))

承压含水层非稳定流有限层分析

王旭东 宰金珉

(1. 南京工业大学 岩土工程研究所,江苏 南京 210009;2. 浙江省建筑科学设计研究院,浙江 杭州 310012)

摘 要:推导了以伽辽金法结合傅立叶函数为基础的各向异性承压含水层非稳定流有限层方程,编制了相应的计算程序。通过算例的有限层解与解析解的对比分析,在验证有限层法正确性的基础上,深入探讨了计算区域、计算项数等因素对计算结果的影响,得到了计算区域尺寸影响有限层解的收敛速度等结论,最后给出了一定精度条诸宏博

1,211

件下计算项数选用的近似估算公式。

关键词:有限层法;承压含水层;非稳定流;降深

0 引 言

有限层法是一种对某一方向进行数值离散,而在其余两方向采用连续函数的半数值半解析方法,这种方法能有效地将三维问题简化为一维问题求解,从而大大减少单元数,节省计算工作量和内存的需要量。基于这些特点,Cheung在1968年首先提出了有限条法,并成功地将有限条法应用于结构

工程问题的分析[1]

。Small和Booker等运用了有限层法对层状弹性地基和土的固结问题做了较为系统

的研究[2-4]

。A.M. oskoorouchi和 B. Novrouzian将有限条方法推广应用于岩土工程分析中,对基坑开

挖、隧道和土坝等问题进行了研究分析[5]

。Stanley,

Myron等将有限层理论应用于地下水分析[6]

。张问清、赵锡宏、宰金珉等采用有限层方法探讨了任意

力系作用下的层状弹性半空间问题[7]

。宰金珉、王旭东等将有限层法应用于土与结构物的共同作用分析[8-9]。

地下水计算方法主要有解析法和数值法。解析法可以将描述地下水运动的各种物理量反映在一个表达式中,使用方便,但自然界地下水运动状态复杂,往往与解析解的假设条件相差甚远,其应用范围也受到一定的限制。对于一些复杂的地下水问题,就必须借助于数值模拟方法求其近似解。目前,地下水数值分析中主要有有限差分法、有限元法、边

界元法等数值方法[10-12]

。数值法能较为方便地处理解析解难以解决的复杂边界等问题,因此,得到广泛的应用。但是,数值法计算工作量大,尤其是在考虑三维模型的计算分析时,占用内存多。国内外学者不断尝试寻找有效的数值方法用于地下水模拟计算,薛禹群,叶淑君,谢春红等将多尺度有限元法应用于地下水模拟,节省了计算量提高了计算精度[13]。

本文吸取有限层法的优点,以伽辽金法结合傅立叶函数为基础,推导了承压含水层非稳定流有限

层方程,编制了相应的计算程序,在验证方法的正确性基础上,探讨了计算区域、计算项数等因素对计算结果的影响。

1 承压含水层非稳定流模型的数学描述

假设均质、各向异性、等厚的有限区域承压含水层,在区域中心(x0 , y0)处有一个以恒定流量Q抽水的完整井,初始水位水平,定水头边界条件,如图1所示。根据质量守恒定律和达西定律,承压含水层非稳定流的数学模型用降深可表示为:

KssssxxxyKyyzKzzR(x,y,z,t)S

s

t

t>0, a≤x≤a,b≤y≤b

(1)



s(x,y,z,t)0 t0, a≤x≤a,b≤y≤bs(x,y,z,t)x-a0 t>0, b≤y≤bxa



s(x,y,z,t)y-b0 t>0, a≤x≤ayb式中,s(x,y,z,t)为降深,Kx、Ky、Kz分别为x,y,z

方向的渗透系数,Ss为单位贮水系数,R(x,y,z,t)为源汇项,补给源取正值,排泄汇取负值。

(1)式可用微分算子表示为:

LsR0 (2) 式中,Ls

s



s



s

s

x

Kxx

y

Kyy

z

Kzz

Sst,L为微

分算子。

───────

收稿日期:

基金项目:国家自然科学基金项目(50278042);江苏省高校自然科学研究计划项目(03KJB560044) 作者简介:诸宏博(1981-),男(汉族),浙江杭州人,硕士研究生,主要研究方向为环境岩土工程及数值分析。E-mail: [email protected]

图1 承压含水层物理模型与有限层层元 Fig.1 Configuration of confined aquifer physical model

and finite layer models

2 有限层格式推导

有限层计算简图如图1所示。x和y方向的计算尺寸分别取为2a和2b,沿z轴方向将承压含水层离散成L个层元。

2.1 试探函数

设降深在x和y方向分别用解析函数表示,而z方向采用有限层元离散,层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数: L1M

N

s

x,y,z,tΦmnjtAmn(x,y)Nj(z) (3) j1m1n1式中,~sx,y,z,t为降深试探函数;Φmnj (t)为待求系数;Amn(x,y)sinkmxasinknyb为满足定水头边界的已知函数,kmm/2a,knn/2b;M,N为级数项数;L为有限层层元数。

Ni(z)为线性插值基函数,其形式如下:

zzi1/zizi1zi1≤zNz≤zizz (4) ii1/zizi1zi≤z≤z

i1

0

z≤zi1或 z≥zi1

2.2 伽辽金方程

根据加权余量法原理,可得有限层法伽辽金方程:

Ls~

D

Rmni

dxdydz0,其中 i1,2,,L1 (5)

式中,mni为权函数,mnix,y,zNizAmn(x,y)。

将含水层离散为L个层元,则(5)式可改写为:

L

Ls~

Rmni

dxdydz0,其中i1,2,,L1 (6)

e1

De

将(3)式代入(6)式得:

LL1MN

e

e

e

LΦmnjtAmnNj

Re1j1m1n1

D

Ne

iAmndxdydz0

(7)

其中 i1,2,,L12.3 层元有限层方程

取一典型有限层层元l,其节面分别为l,l+1,层厚c,如图2所示。

图2 典型有限层层元

Fig.2 Configuration of finite layer element models

由Amn(x, y)和Am′n′(x, y) 的正交性,可得:

a

 0 asinkm

xsinkmxdxmm

a mm

(8)

b

 bsinkn

ysink0 nn

nydyb nn

(9)

故,当 mm和nn时,由(7)式得:

l1M

N

ΦZl1mnjKNeNedza

jlm1n1Zlxijab

2Amn

bx2Amndxdy

b

2 Φl1

eea

mnj

Z(10)

ZKl

yNNdzAmn

ij

aby2Amndxdy

ΦNej

b

mnjKzzNec

i0

a

a

b

A2mndxdy ΦZel1

KNNe

ij

a

b

mnjZl

z

zz

dz

a

b

A2mndxdy

Φmnj

Zl1

t

eZSA2l

sNiNejdz

a

a

b

b

mndxdy

Zl1a

b

Zl

a

b

RlNejAmndxdydz0, 其中 il,l1

由基函数的性质可得:

Zl1

ZNeNecl3 ijij

dzl

cij

(11) l 由(11)式可分别得到m、n阶的层元渗透矩阵系

数ge

、贮水矩阵系数beemnijmnij和水量矩阵系数qmnij:

g

e



Zl1

eea

b

mnij

ZKNiNj

dzl

xa2Amnb

x2Amndxdy



Zl1

Keea

b

2ZyNiNjdz

Amnb

y

2

Alamndxdy (12a)

e

KN

i

Zab

zz

Nj0

a

b

A2mndxdy 

Zl1

KNeNe

ija

b

Zl

zzz

dz

a

b

A2mndxdy

l1

b

eZl1

eemnij

sNiNjdz (12b)

il

ZSl

a

a

b

b

A2

mndxdyqe1

a

b

emnj

ZlZRlNjAmndxdydz

(12c) l

a

b

式中,i=l,l+1,j=l,l+1。

由(12a)式、(12b)式积分可得:

abK22xlkmKylkncl

Kzl ij

ge



3cl

mnij



ab

Kk2

xlmKylk2

n

c (13a)

lKz6lc ijl

abclSsij be

mnij

3abc (13b)

lSs

6 ij对于位于区域中心处的完整抽水井,由(12c)

式积分可得: l1b

qemnj

zzNyl

jdz

a

a

b

x0,0x,yAmndxdy

(13c)

cl

2

x0,y0sinkmasinknb, 其中jl,l1式中,x0,y0为抽水率x0,y0Qx0,y0/cl。

层元l的有限层方程为:

gel,lgebel,l1ll,lbel,lgeqel

el1,l

ge1

dll1,l1mnl1mnbel1,lbl1,l1mndtl1mn

qe0(14) l1mn

用矩阵可简化表示为:

Ge

e

de

mnΦmnBmndt

ΦmnQmn0 (15) 2.4 整体有限层方程

由层元有限层方程组装可得整体有限层方程:

GdmnΦmnBmndt

ΦmnQmn0 (16)

L

L

L

式中,GemnG,Bmne1

Be,Qmne1

Qe

mn

mn

mn,

e1

ΦmnT

12

jLL1mn。

整体矩阵Gmn、Bmn由层元矩阵叠加而成,叠加方法如图3所示。

图3 [G]mn和[B]mn矩阵叠加 Fig.3 Matrix assembly for [G]mn and [B]mn

3 有限层方程求解

对于非稳定流来说,降深是时间t的函数,则

Φmnj(t)也是时间t的函数,

采用差分法对时间进行离散得:

dΦk1k

mn

ΦmnΦmndt

t

(17)

假设在时间k和k+1之间,水头呈线性变化,

则有:

Φk

mn

Φk1k

mn1Φmn

(18) 式中, Φk1

mn和Φkmn分别为k时刻和k+1时刻的值,t为时间步长,为加权因子,取值范围在0~1之间。

将(18)式代入(16)式可得有限层方程的差分格式:

(G1k1

mn

t

Bmn)Φmn

(19) (1kkt

Bmn1Gmn)ΦmnQmn

在(19)式中,当0时,即为显式差分格式;1时,为隐式差分格式;0.5时,为Crank-Nicolson差分格式,其中,隐式差分格式和Crank-Nicolson差分格式无条件稳定。

(19)式有限层方程的差分格式为L+1阶线性代数方程组,求解方程组可得到未知系数Φmnj (t),并由(3)式计算k+1时刻的降深。

数值计算时初始降深为0,故Φ0

mn0,计算程序流程图如图4所示,并编制了相应的计算程序,采用追赶法

[14]

求解方程,以节省计算时间。

图4 有限层计算程序流程图

Fig.4 Program flow graph of finite layer method

4 数值算例及分析

4.1 有限承压含水层中抽水井

算例:承压含水层厚c=100m,为了便于对比分析,渗透系数Kx=Ky=Kz=5m/d,单位贮水系数Ss=0.2×10-4m-1;计算区域取为2a×2b=5000m×5000m。(1).在计算区域中心布置一口单井,单井抽水量Q =2000 m3/d;(2).在计算区域中200m×200m范围内布置8口群井,单井抽水量Q =800m3/d。分别采用有限层法模拟计算抽水引起的降深。

为了验证承压含水层非稳定流有限层分析方法的正确性,引入无限承压含水层非稳定流完整井解

析解[12]

(泰斯公式)进行对比分析。

图5和图6分别给出了单井和群井降深随距离的分布曲线。由图可知,无论是单井还是群井,有限层解与解析解都有较好的吻合程度,特别是在抽水较短时间内,降深不受边界影响,此时两种方法的计算结果尤为一致,表明了本文方法的正确性。 4.2 计算区域的影响

由于本算例为有限含水层问题,所以有限层分析中,随着抽水时间的延长,降深逐渐受到边界的影响,尤其是边界附近的降深值,在图6中,抽水10天时,边界附近的有限层解与解析解之间已存在

一定的差异。

图7 给出了不同计算区域情况下距抽水井100m处有限层和解析解的标准降深曲线,两者对比表明,在抽水初期,抽水引起水位下降的区域较小,基本不受边界的影响,不同计算区域的有限层解与解析解都有较好的一致性,且具有一致性的区域(时间) 随着计算区域的增大而增加,表明有限区域内的有限层解随着计算区域的增大逐渐逼近无限区域的解析解(泰斯公式)。

m

/ 深降距抽水井的距离 / m

图5 单井降深分布曲线 Fig.5 The drawdown curve of single well

距群井中心的距离/m

m

/深降

图6 群井降深分布曲线

Fig. 6 The drawdown curve of multiple wells

cKsQ

41/u = 4Kt/r2

SS

图7 标准化降深曲线(距井100m)

Fig.7 Normalized drawdown curves (100 m from well location)

4.3 计算项数的影响

为了避免边界对计算结果的影响,选取t =1d时的有限层和解析解的计算结果,以讨论计算项数对有限层计算结果的影响,由图7可知,t =1d时可忽略计算区域对有限层计算结果的影响。从图8给出的不同计算区域时有限层解的收敛性过程可知,计算结果的精度随计算项数的增加而呈周期性变化,并逐渐提高;随着计算区域的增大,计算结果的收敛速度减慢。根据有限层解与解析解的对比分析,当计算项数M(或N)取近似估算公式a/25或

b/25

解层解

限析有解M(或N)

图8 M、N与降深值收敛的关系曲线 Fig. 8 Convergence of drawdown value with M and N

5 结 论

(1) 以伽辽金法结合傅立叶函数为基础,推导了承压含水层非稳定流有限层方程,并编制了相应的计算程序。

(2) 算例分析验证了有限层方法的正确性;同时,表明了该方法对地下水非稳定流模拟具有较好的适用性,为有限区域承压含水层非稳定流计算提供一条新的有效途径。

(3) 计算区域和计算项数对计算结果影响的分析表明,计算结果的收敛速度随计算区域的增大而减慢,当计算项数M(或N)取a/25或b/25时能使有限层解的精度保持在±3%以内。 参考文献:

[1] Cheung Y. K., Finite strip method in structural mechanics [M]. Pergamon, New York,1976

[2] Small J C, Booker J R. Finite layer analysis of layered elastic materials using a flexibility approach [J]. Int. J. Numer. Methods Eng.,1984, 20(6):1025-1037.

[3]Booker J.R., Small J.C.,Finite layer analysis of consolidation I [J]. International Journal for Numerical Analysis Methods Geomechanics.1982a,6(2),151-171.

[4]Booker J.R., Small J.C.,Finite layer analysis of consolidation II [J]. International Journal for Numerical Analysis Methods Geomechanics.1982b,6(2),173-194.

[5] A.M. oskoorouchi, B. Novrouzian. Zoned finite strip method and its application in geomechanics [J]. Computers and Geotechnics, 1991,11(2):265-294.

[6] Stanley Smith S, Myron B. Allen, Jay Puckett. The finite layer method for groundwater flow models [J].Water Resources Research. 1992, 28(6): 1715-1722.

[7] 张问清, 赵锡宏, 宰金珉. 任意力系作用下的层状弹性半空间的有限层分析方法[J].岩土工程学报,1981, 3(2):27-42. (ZHANG Wen-qing, ZHAO Xi-hong, ZAI Jin-min. Finite layer analysis for layered and elastic half-space under arbitrary force systems [J]. Chinese Journal of Geotechnical Engineering, 1981, 3(2):27-42. (in Chinese))

[8] 宰金珉, 宰金璋. 高层建筑基础分析与设计[M].北京:中国建筑工业出版社,1993. (ZAI Jin-min, ZAI Jin-zhang. Foundation analysis and design of high-rise building[M]. Beijing: China Architecture and Building Press, 1993. (in Chinese))

[9] 王旭东,魏道垛,宰金珉.群桩-土-承台结果共同作

用的数值分析[J].岩土工程学报,1996,18(4),27-33. (WANG Xu-dong ,WEI Dao-duo, ZAI Jin-min. Numerical Analysis of Pile Groups-Soil-Pile Cap Interaction[J].Chinese Journal of Geotechnical Engineering, 1996, 18(4) :27-23. (in Chinese)) [10] 朱学愚,钱孝星,刘新仁.地下水资源评价[M],南京:南京大学出版社,1987 (ZHU Xue-yu, QIAN Xiao-xing, LIU Xin-ren. Groundwater Resource Evaluation[M]. Nanjing: Nanjing University Press, 1987. (in Chinese))

[11] 陈文华,王旭东.二维非均质体中非稳定流问题的边界单元法[J].岩土工程学报,1991,13(3),65-71. (CHEN Wen-hua, WANG Xu-dong . 2D Analysis of Unsteady groundwater in Non- homogeneous

aquifer by Boundary Element

method[J].Chinese Journal of Geotechnical Engineering,

1991,13(3),65-71. (in Chinese))

[12] 孙讷正. 地下水流的数学模型和数值方法[M].北京:地质出版社,1981. (SUN Ne-zheng. Mathematical Model and Numerical Method of Groundwater Flow[M]. Beijing: Geological Publishing House, 1981. (in Chinese))

[13] 薛禹群,叶淑君,谢春红,张云. 多尺度有限元法在地下水模拟中的应用[J].水利学报,2004(7),7-13. (XUE Yu-qun, YE Shu-jun, XIE Chun-hong, ZHANG Yun. Application of multi scale finite element method to simulation of groundwater flow[J]. Journal of Hydraulic Engineering, 2004(7),7-13. (in Chinese))

[14] 徐士良. FORTRAN常用算法程序集(第二版)[M].北京:清华大学出版社,1995. (XU Shi-liang. FORTRAN Common Algorithms Program (Edition Two) [M]. Beijing: Tsinghua University Press, 1995. (in Chinese))


相关文章

  • 承压水地层基坑底部突涌及解决措施
  • 23第23卷 第5期隧道建设(5):25-27 2003年10月TunnelConstructionOct,2003 承压水地层基坑底部突涌及解决措施 郑剑升1 张克平2 章立峰1 (1 同济大学铁道建筑工程系 上海2003312 中铁二局 ...查看


  • 采掘场地下水涌水量计算
  • 采矿工程 露天采矿技术 2007年第4期・39・ 采掘场地下水涌水量计算 金祎民 (中煤国际工程集团沈阳设计研究院,辽宁沈阳110015) 摘要:地下水涌水量预测是确定采掘场疏干系统布置的前提.对胜利东二号露天区的水文地质条件进 行分析,确 ...查看


  • 承压含水层减压降水对既有盾构隧道影响研究
  • 第36卷第1期 岩 土 力 学 Vol.36 No. 1 2015年1月 Rock and Soil Mechanics Jan. 2015 DOI:10.16285/j.rsm.2015.01.025 承压含水层减压降水对既有盾构隧道影响 ...查看


  • 钱塘江潮汐作用下地层降水试验研究
  • 摘 要:本文以杭州市望江路过江隧道项目降水试验为依托,对钱塘江潮汐作用下,复杂地层承压水降水参数以及地层渗透系数进行研究,对水文地质参数进行参数反演,对下一步深基坑降水施工进行实际指导.该工程基坑降水井单井出水量大,周边环境复杂敏感,通过抽 ...查看


  • 基坑降水(承压含水层减压降水)施工组织方案及技术保证措施
  • 基坑降水施工组织方案及技术保证措施 1基坑降水分析 本工程降水分为坑内开挖范围内疏干性降水和深层承压含水层减压降水.针对这两类降水分别进行降水设计和现场施工运行. 1.1含水层降水风险分析 根据工程场地工程地质条件与水文地质条件分析,可能引 ...查看


  • 地下结构抗浮设防水位和浮力计算
  • 第3l卷第4期 2009年11月河北埋工大学学报(目然科学版)JournalofHebeiPolytechnicUniversity(NaturalScienceEdition)VoL31no.4Nov.2009文章编号:1674-0262 ...查看


  • (工程地质)答案
  • 1. 内外力地质作用定义及其包含具体内容有哪些,以及内外力作用的特征及其之间差别: 外力地质作用:指以太阳能以及日月引力能为能源并通过大气,水,生物等因素引起的地质作用,包括风化作用,剥蚀作用,搬运作用,沉积作用,固结成岩作用. 内力地质作 ...查看


  • 登封市向阳煤业防隔水煤柱设计doc
  • 河南金丰煤业集团有限公司 向阳煤业防隔水煤柱设计批复通知 向阳煤业: 河南金丰煤业集团有限公司总工办通过对<登封市向阳煤业有限公司技术改造初步设计说明书>会审,并根据<建筑.水体.铁路及主要井巷煤柱留设及压煤开采规程> ...查看


  • 地下水水力特征对土自重应力计算的影响分析
  • 刘振京等:地下水水力特征对土自重应力计算的影响分析 ・37・ 地下水水力特征对土自重应力计算的影响分析 刘振京孙刚 (河北工程技术高等专科学校河北沧州061001) 摘要利用太沙基的有效应力原理,根据地层中地下水的不同水力特征,论述了地下水 ...查看


热门内容