承压含水层非稳定流有限层分析
王旭东 宰金珉
(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所示。根据质量守恒定律和达西定律,承压含水层非稳定流的数学模型用降深可表示为:
KssssxxxyKyyzKzzR(x,y,z,t)S
s
t
t>0, a≤x≤a,b≤y≤b
(1)
s(x,y,z,t)0 t0, a≤x≤a,b≤y≤bs(x,y,z,t)x-a0 t>0, b≤y≤bxa
s(x,y,z,t)y-b0 t>0, a≤x≤ayb式中,s(x,y,z,t)为降深,Kx、Ky、Kz分别为x,y,z
方向的渗透系数,Ss为单位贮水系数,R(x,y,z,t)为源汇项,补给源取正值,排泄汇取负值。
(1)式可用微分算子表示为:
LsR0 (2) 式中,Ls
s
s
s
s
x
Kxx
y
Kyy
z
Kzz
Sst,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方向采用有限层元离散,层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数: L1M
N
s
x,y,z,tΦmnjtAmn(x,y)Nj(z) (3) j1m1n1式中,~sx,y,z,t为降深试探函数;Φmnj (t)为待求系数;Amn(x,y)sinkmxasinknyb为满足定水头边界的已知函数,kmm/2a,knn/2b;M,N为级数项数;L为有限层层元数。
Ni(z)为线性插值基函数,其形式如下:
zzi1/zizi1zi1≤zNz≤zizz (4) ii1/zizi1zi≤z≤z
i1
0
z≤zi1或 z≥zi1
2.2 伽辽金方程
根据加权余量法原理,可得有限层法伽辽金方程:
Ls~
D
Rmni
dxdydz0,其中 i1,2,,L1 (5)
式中,mni为权函数,mnix,y,zNizAmn(x,y)。
将含水层离散为L个层元,则(5)式可改写为:
L
Ls~
Rmni
dxdydz0,其中i1,2,,L1 (6)
e1
De
将(3)式代入(6)式得:
LL1MN
e
e
e
LΦmnjtAmnNj
Re1j1m1n1
D
Ne
iAmndxdydz0
(7)
其中 i1,2,,L12.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
xsinkmxdxmm
a mm
(8)
b
bsinkn
ysink0 nn
nydyb nn
(9)
故,当 mm和nn时,由(7)式得:
l1M
N
ΦZl1mnjKNeNedza
jlm1n1Zlxijab
2Amn
bx2Amndxdy
b
2 Φl1
eea
mnj
Z(10)
ZKl
yNNdzAmn
ij
aby2Amndxdy
ΦNej
b
mnjKzzNec
i0
a
a
b
A2mndxdy ΦZel1
KNNe
ij
a
b
mnjZl
z
zz
dz
a
b
A2mndxdy
Φmnj
Zl1
t
eZSA2l
sNiNejdz
a
a
b
b
mndxdy
Zl1a
b
Zl
a
b
RlNejAmndxdydz0, 其中 il,l1
由基函数的性质可得:
Zl1
ZNeNecl3 ijij
dzl
cij
(11) l 由(11)式可分别得到m、n阶的层元渗透矩阵系
数ge
、贮水矩阵系数beemnijmnij和水量矩阵系数qmnij:
g
e
Zl1
eea
b
mnij
ZKNiNj
dzl
xa2Amnb
x2Amndxdy
Zl1
Keea
b
2ZyNiNjdz
Amnb
y
2
Alamndxdy (12a)
e
KN
i
Zab
zz
Nj0
a
b
A2mndxdy
Zl1
KNeNe
ija
b
Zl
zzz
dz
a
b
A2mndxdy
l1
b
eZl1
eemnij
sNiNjdz (12b)
il
ZSl
a
a
b
b
A2
mndxdyqe1
a
b
emnj
ZlZRlNjAmndxdydz
(12c) l
a
b
式中,i=l,l+1,j=l,l+1。
由(12a)式、(12b)式积分可得:
abK22xlkmKylkncl
Kzl ij
ge
3cl
mnij
ab
Kk2
xlmKylk2
n
c (13a)
lKz6lc ijl
abclSsij be
mnij
3abc (13b)
lSs
6 ij对于位于区域中心处的完整抽水井,由(12c)
式积分可得: l1b
qemnj
zzNyl
jdz
a
a
b
x0,0x,yAmndxdy
(13c)
cl
2
x0,y0sinkmasinknb, 其中jl,l1式中,x0,y0为抽水率x0,y0Qx0,y0/cl。
层元l的有限层方程为:
gel,lgebel,l1ll,lbel,lgeqel
el1,l
ge1
dll1,l1mnl1mnbel1,lbl1,l1mndtl1mn
qe0(14) l1mn
用矩阵可简化表示为:
Ge
e
de
mnΦmnBmndt
ΦmnQmn0 (15) 2.4 整体有限层方程
由层元有限层方程组装可得整体有限层方程:
GdmnΦmnBmndt
ΦmnQmn0 (16)
L
L
L
式中,GemnG,Bmne1
Be,Qmne1
Qe
mn
mn
mn,
e1
ΦmnT
12
jLL1mn。
整体矩阵Gmn、Bmn由层元矩阵叠加而成,叠加方法如图3所示。
图3 [G]mn和[B]mn矩阵叠加 Fig.3 Matrix assembly for [G]mn and [B]mn
3 有限层方程求解
对于非稳定流来说,降深是时间t的函数,则
Φmnj(t)也是时间t的函数,
采用差分法对时间进行离散得:
dΦk1k
mn
ΦmnΦmndt
t
(17)
假设在时间k和k+1之间,水头呈线性变化,
则有:
Φk
mn
Φk1k
mn1Φmn
(18) 式中, Φk1
mn和Φkmn分别为k时刻和k+1时刻的值,t为时间步长,为加权因子,取值范围在0~1之间。
将(18)式代入(16)式可得有限层方程的差分格式:
(G1k1
mn
t
Bmn)Φmn
(19) (1kkt
Bmn1Gmn)ΦmnQmn
在(19)式中,当0时,即为显式差分格式;1时,为隐式差分格式;0.5时,为Crank-Nicolson差分格式,其中,隐式差分格式和Crank-Nicolson差分格式无条件稳定。
(19)式有限层方程的差分格式为L+1阶线性代数方程组,求解方程组可得到未知系数Φmnj (t),并由(3)式计算k+1时刻的降深。
数值计算时初始降深为0,故Φ0
mn0,计算程序流程图如图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
cKsQ
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所示。根据质量守恒定律和达西定律,承压含水层非稳定流的数学模型用降深可表示为:
KssssxxxyKyyzKzzR(x,y,z,t)S
s
t
t>0, a≤x≤a,b≤y≤b
(1)
s(x,y,z,t)0 t0, a≤x≤a,b≤y≤bs(x,y,z,t)x-a0 t>0, b≤y≤bxa
s(x,y,z,t)y-b0 t>0, a≤x≤ayb式中,s(x,y,z,t)为降深,Kx、Ky、Kz分别为x,y,z
方向的渗透系数,Ss为单位贮水系数,R(x,y,z,t)为源汇项,补给源取正值,排泄汇取负值。
(1)式可用微分算子表示为:
LsR0 (2) 式中,Ls
s
s
s
s
x
Kxx
y
Kyy
z
Kzz
Sst,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方向采用有限层元离散,层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数: L1M
N
s
x,y,z,tΦmnjtAmn(x,y)Nj(z) (3) j1m1n1式中,~sx,y,z,t为降深试探函数;Φmnj (t)为待求系数;Amn(x,y)sinkmxasinknyb为满足定水头边界的已知函数,kmm/2a,knn/2b;M,N为级数项数;L为有限层层元数。
Ni(z)为线性插值基函数,其形式如下:
zzi1/zizi1zi1≤zNz≤zizz (4) ii1/zizi1zi≤z≤z
i1
0
z≤zi1或 z≥zi1
2.2 伽辽金方程
根据加权余量法原理,可得有限层法伽辽金方程:
Ls~
D
Rmni
dxdydz0,其中 i1,2,,L1 (5)
式中,mni为权函数,mnix,y,zNizAmn(x,y)。
将含水层离散为L个层元,则(5)式可改写为:
L
Ls~
Rmni
dxdydz0,其中i1,2,,L1 (6)
e1
De
将(3)式代入(6)式得:
LL1MN
e
e
e
LΦmnjtAmnNj
Re1j1m1n1
D
Ne
iAmndxdydz0
(7)
其中 i1,2,,L12.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
xsinkmxdxmm
a mm
(8)
b
bsinkn
ysink0 nn
nydyb nn
(9)
故,当 mm和nn时,由(7)式得:
l1M
N
ΦZl1mnjKNeNedza
jlm1n1Zlxijab
2Amn
bx2Amndxdy
b
2 Φl1
eea
mnj
Z(10)
ZKl
yNNdzAmn
ij
aby2Amndxdy
ΦNej
b
mnjKzzNec
i0
a
a
b
A2mndxdy ΦZel1
KNNe
ij
a
b
mnjZl
z
zz
dz
a
b
A2mndxdy
Φmnj
Zl1
t
eZSA2l
sNiNejdz
a
a
b
b
mndxdy
Zl1a
b
Zl
a
b
RlNejAmndxdydz0, 其中 il,l1
由基函数的性质可得:
Zl1
ZNeNecl3 ijij
dzl
cij
(11) l 由(11)式可分别得到m、n阶的层元渗透矩阵系
数ge
、贮水矩阵系数beemnijmnij和水量矩阵系数qmnij:
g
e
Zl1
eea
b
mnij
ZKNiNj
dzl
xa2Amnb
x2Amndxdy
Zl1
Keea
b
2ZyNiNjdz
Amnb
y
2
Alamndxdy (12a)
e
KN
i
Zab
zz
Nj0
a
b
A2mndxdy
Zl1
KNeNe
ija
b
Zl
zzz
dz
a
b
A2mndxdy
l1
b
eZl1
eemnij
sNiNjdz (12b)
il
ZSl
a
a
b
b
A2
mndxdyqe1
a
b
emnj
ZlZRlNjAmndxdydz
(12c) l
a
b
式中,i=l,l+1,j=l,l+1。
由(12a)式、(12b)式积分可得:
abK22xlkmKylkncl
Kzl ij
ge
3cl
mnij
ab
Kk2
xlmKylk2
n
c (13a)
lKz6lc ijl
abclSsij be
mnij
3abc (13b)
lSs
6 ij对于位于区域中心处的完整抽水井,由(12c)
式积分可得: l1b
qemnj
zzNyl
jdz
a
a
b
x0,0x,yAmndxdy
(13c)
cl
2
x0,y0sinkmasinknb, 其中jl,l1式中,x0,y0为抽水率x0,y0Qx0,y0/cl。
层元l的有限层方程为:
gel,lgebel,l1ll,lbel,lgeqel
el1,l
ge1
dll1,l1mnl1mnbel1,lbl1,l1mndtl1mn
qe0(14) l1mn
用矩阵可简化表示为:
Ge
e
de
mnΦmnBmndt
ΦmnQmn0 (15) 2.4 整体有限层方程
由层元有限层方程组装可得整体有限层方程:
GdmnΦmnBmndt
ΦmnQmn0 (16)
L
L
L
式中,GemnG,Bmne1
Be,Qmne1
Qe
mn
mn
mn,
e1
ΦmnT
12
jLL1mn。
整体矩阵Gmn、Bmn由层元矩阵叠加而成,叠加方法如图3所示。
图3 [G]mn和[B]mn矩阵叠加 Fig.3 Matrix assembly for [G]mn and [B]mn
3 有限层方程求解
对于非稳定流来说,降深是时间t的函数,则
Φmnj(t)也是时间t的函数,
采用差分法对时间进行离散得:
dΦk1k
mn
ΦmnΦmndt
t
(17)
假设在时间k和k+1之间,水头呈线性变化,
则有:
Φk
mn
Φk1k
mn1Φmn
(18) 式中, Φk1
mn和Φkmn分别为k时刻和k+1时刻的值,t为时间步长,为加权因子,取值范围在0~1之间。
将(18)式代入(16)式可得有限层方程的差分格式:
(G1k1
mn
t
Bmn)Φmn
(19) (1kkt
Bmn1Gmn)ΦmnQmn
在(19)式中,当0时,即为显式差分格式;1时,为隐式差分格式;0.5时,为Crank-Nicolson差分格式,其中,隐式差分格式和Crank-Nicolson差分格式无条件稳定。
(19)式有限层方程的差分格式为L+1阶线性代数方程组,求解方程组可得到未知系数Φmnj (t),并由(3)式计算k+1时刻的降深。
数值计算时初始降深为0,故Φ0
mn0,计算程序流程图如图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
cKsQ
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))