潜水含水层非稳定井流降水的有限差分分析
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Key words: non-steady; unconfined aquifer; hydro-geological parameters; precipitation; finite-difference
采用井点降水预压进行软基处理，其关键问 题在于降水设计。Dupuit 认为潜水井中进行定流 量（或定降深）抽水，经过一定时间 t，渗流将趋向 稳定，潜水面 由原来的水 平状态变为 漏斗状 [1]。 Muskat M 给 出 达 到 稳 定 的 时 间 t ≥0.5R2Sy/（kH） （Sy 为潜水含水层的给水度，R 为影响半径，k 为 土层渗透系数，H 为潜水含水层厚度）[2]。因此， 通过稳定井流抽水试验估计 k 值耗时过长。Theis C V 提出地下水径向非稳定承压井流计算公式[3]， 已知抽水量 Q，降深 s，经过抽水时间 t，可求得 土层压力传导系数 a，水力传导系数 T（=kM）和承
图 1 插设塑料排水板后的 s~t 曲线 图 2 插设塑料排水板后的降水漏斗时程曲线
·72·
水运工程
2010 年
图 3 插设塑料排水板后的 Neuman 拟合曲线
限差分公式：
CR i,j-
1
，k
(hi,j-1,k-hi,j,k)+CR i,j+
1
(hi,j,k-hi,j,k)+
，k
2
2
CC i,j-
计算过程分 8 个时间阶段：第 1~7 个阶段陆 续完成 24 个泵的安装，每安装完 1 个泵便开始抽 水，第 8 个阶段所有泵都开始抽水，整个过程共 历时 25 d。
第2期
刘 嘉，等：潜水含水层非稳定井流降水的有限差分分析
·73·
表 1 有限差分法计算参数一览
土层
kxx /
kyy /
kzz /
（cm·s-）1 （cm·s-）1 （cm·s-）1
（3）
式中符号意义见文献[19]。式（3）以矩阵形式可表
示为：
Ah=q
（4）
式中：A 为水头Biblioteka Baidu系数矩阵；h 为所求得水头矩
阵；q 表示各方程中所包含的所有常数项和已知
项，求解时根据 A 和 q 通过迭代获得 h。
2.4 有限差分法模型
如图 4 和图 5 所示，三维地下水渗流模型尺
寸为 300 m×245 m×20 m，加固区位于模型中部
2.3 有限差分法基本原理
在不考虑水的密度变化的条件下，孔隙介质
中地下水在三维空间的流动可以用下面的偏微分
方程来表示[19]：
! " ! " ! " 坠
坠x
kxx
坠h 坠x
+坠 坠y
kyy
坠h 坠y
+坠 坠z
kzz
坠h 坠z
-W =Ss
坠h 坠t
（2）
式中：kxx，kyy 和 kzz 分别为渗透系数在 x，y 和z 方
1 2
,j,k
(h-Ii,j,k)+CCj+
1 2
，j,k(hi+1,j,k-hi,j,k)+
CVi,j,k-
1 2
(hi,j,k-1-hi,j,k)+CV i,j,k+
1 2
(hi,j,k+1-hi,j,k)+
Pi,j,khi,j,k+Qi,j,k=Ssi,j,k(△rj△ci△υk)
△hi,j,k △t
Abstract: In the coastal area, an unconfined aquifer will be formed on the upper part of the foundation by
dredging and filling silts or sands to reclaim land, but the most existing methods are relying on the Theis, a well flow theory of non-steady confined aquifer to solve the hydro-geological parameters and conduct precipitation design that is lack of rationality. By applying the Neuman unsteady well flow theory to derive the hydro geological parameters of unconfined aquifer and then using the finite-difference method to analyze all derived parameters, the resulted values are consistent with the water level monitored results. The engineering examples show that this method is reasonable and reliable, and provides a way for the similar applications.
量 q 取实测值 48 m3/d，井内最大水位降深 15 m。
根据工程地质条件，将试验区含水层简化为自上
而下的 5 个土层：①0~6 m 第 1 个软软弱黏土层；
②6~8 m 砂层；③8~13 m 第 2 个超软弱黏土层；
④13~15 m 砂层；⑤15~20 m 硬黏土层，并将区内
视为单一潜水含水层，采用 Neuman 法求参结果，
Sy
加固区内等效土层 4.36×10-4 4.36×10-4 4.36×10-5 0.15
区外土层①，淤泥
1×10-7
1×10-7
1×10-8
0.02
区外土层②，砂层
1×10-2
1×10-2
1×10-3
0.26
区外土层③，淤泥
1×10-7
1×10-7
1×10-8
0.02
区外土层④，砂层
1×10-2
1×10-2
中图分类号：U 655.54+3
文献标志码：A
文章编号：1002-4972（2010）02-0070-05
Finite-difference analysis of unsteady well flow at unconfined aquifer
LIU Jia, ZHANG Gong-xin, LUO Yan, DONG Zhi-liang (CCCC Sihang Engineering Research Institute Co., Ltd., Guangzhou 510230, China)
收稿日期：2009-08-20 作者简介：刘嘉（1981—），男， 硕士，主要从事岩土工程数值分析方面的研究工作。
第2期
刘 嘉，等：潜水含水层非稳定井流降水的有限差分分析
·71·
目标值拟合井函数 W（u），形成了非稳定井流抽水 试验确定水文地质参数的解析方法[8]，较之图解法 和配线法更为简便。齐学斌按 Theis 基本方程，运 用 Gauss-Newton 迭代 法进行非线 性回归求参 [9]。 刘立才等借助梯度数学模型对 Theis 公式求解得出 T 和 S[10]。刘玉珍等基于 Theis 方程论述了灵敏度分 析法确定水文地质参数的基本原理，即把 s 视为 以 T 和 S 为变量的函数，利用最小二乘法求出理 论降深与观测降深最佳拟合时的 T 和 S[11]。巩彦文 通过对 W（u）的优化拟合，建立一个二元线性回归 方程求解回归方程中的待定常数，从而反算含水 层的水文地质参数[12]。类似的方法还见于 Swamee P K，Ojha C S P 和 Ajayi O 等人的研究中[13-17]。
压含水层给水度 S。其中，W（u）为 Theis 井函数，M 为承压含水层厚度，r 为井半径。Jacob C E 建议 用修正值 sc=s-s（/ 2H）描述潜水井流情况[4]，并要求 满足最大降深大于含水层厚度 5%[5]。
根据 Theis 解，由非稳定井流抽水试验所得 s~t 曲线 （s 为水位降深） 确定水文地质参数的主 要方法有图解法、配线法和解析法。前两种方法 人为误差大，求参过程复杂，适用性不大。梁文 康通过编写图解查表程序[6]，邹正盛[7]依据最小二 乘法建立自动配线法，都丰富了传统求解手段。 滕凯采用逐次逼近法，以标准剩余差最小为最优
1×10-8
0.26
区外土层⑤，黏土层 1×10-8
1×10-8
1×10-9
0.02
3 模型验证 如图 6 所示，两个中心点水位观测井的实测
值与有限差分法所得计算值趋势基本一致，在抽 水期的 0~6 d 计算水位高于实测水位，之后的时 间段则是预测值水位低于实测值。实测曲线在 0~ 5 d 之间没有水位观测点，所得的结果之间采用线 性插值，而计算曲线和实测曲线也在抽水后的第 6 d 相交，说明与观测结果相符，仅是实测数据偏 少未能用线性变化反映期间的过程。这一阶段潜 水泵正在逐个安装，区内抽水量较小，水位降低 速率应低于后期，预测曲线恰能更好地反映这点。 预测曲线的第 2 个阶段，地下水位下降速率增加， 而后慢慢变缓。这与非稳定井流理论相符，即认 为潜水面是可移动的物质界面，随着抽水时间的 增加，s~t 曲线斜率减小，有时甚至获得短时间的 “稳定”，这反映了潜水疏干排水的作用，如同含 水层得到“补给”，属于滞后给水现象。实测降深 曲线则是有高有低，一方面是测量过程中人为因 素的影响，另一方面潜水泵在抽水过程中需要定时 维护，也将使抽水量变小，水位出现短时间上升。 就总体趋势而言，预测水位降深略大于实测值。
向上的分量；h 为水头；W 为单位体积流量。离
散化后得到计算单元 （i，j，k） 的地下水渗流有
1 降水设计 1.1 工程地质条件
广州港南沙港区粮食及通用码头工程软基处 理试验区，加固面积 22 500 m（2 125 m×180 m），地 下水位离地表 0.5 m。根据地质勘察资料揭示，原 场地自上而下可分为：①淤泥混砂为主的近期人 工回填土，平均厚度 2 m；②淤泥或淤泥质黏土， 层①+②平均厚度 6 m，为第 1 个软土层；③混砂 层，平均厚度 2 m；④淤泥或淤泥质黏土，平均 厚度 5 m，为第 2 个软土层；⑤中粗或中细砂层， 平均厚度 2 m；⑥分布不均匀的粉质黏土、黏土、 砂层及硬黏土混合层。软基加固前在场地表层铺 设了 1.5 m 厚的砂垫层作为水平排水体及后续施 工垫层，并按 1.1 m 间距插设塑料排水板作为竖 向排水体，以及在加固区四周打设双排桩泥浆搅 拌墙作为止水帷幕。
2.2 水文地质求参
广义 Neuman 方程定义潜水含水层中的水位
降深 s 随时间 t 变化规律如下[18]：
s=
Q 4πt
W(uA，uB，β)
（1）
式中：W (uA，uB，β) 为潜水含水层井函数；uA=r2S/ （4Tt）（用于降水初 期 A 类型曲线），uB=r2Sy/(4Tt)
（用于降水后期 B 类型曲线），T=khH。抽水试验按
现有非稳定井流解析方法大都基于 Theis 系 统，其实质是承压井流，属二维情况，没有考虑 潜水井流存在垂直分流速，属于三维流动的实质。 并且从承压井中抽水是弹性贮量的释放，接近瞬 时完成，无压井则为含水层重力水疏干的释水机 制。针对软基处理工程中潜水井流降水的特点， 本文通过 Neuman 非稳定潜水完整井流求参，并 进行有限差分分析，研究数值计算和实测地下水 位降深规律，为今后相关设计和研究工作提供有 益的参考。
含水层井流理论求解水文地质参数，进行降水设计，缺乏合理性。结合广州港南沙港区粮食及通用码头软基试验区科研项
目，利用 Neuman 非稳定井流理论进行潜水含水层水文地质求参，将所得参数进行有限差分分析，数值结果与水位监测成果
一致。工程实例表明，该方法合理、可靠，为类似应用提供了一种新途径。
关键词：非稳定；潜水含水层；水文地质参数；降水；有限差分
照距离抽水井点 PW 一侧 5 m，10 m 和 15 m 远的
地方分别布设观测井 OW1，OW2 和 OW3。
根据图 1 和图 2 中的试验数据，借助商用软
件包 AquiferTest Pro 求参，可得插设塑料排水板
之后加固区内等效水平渗透系数 kh=4.36×10-4 cm/s，
Sy=0.15，拟合曲线见图 3。
（灰色区域），尺寸为 180 m×125 m×20 m，离散后
的计算单元 （i，j，k）尺寸为 2 m×2 m×2 m。加固区
内共打设 24 个抽水井 （PW1~PW24，间距 30 m×
30 m） 和 2 个水位观测孔 （OWZX1 和 OWZX2），
区外布置 8 个水位观测孔（OW1~OW8）。单井抽水
2010 年 2 月 第 2 期 总第 438 期
水运工程 Port & Waterway Engineering
Feb. 2010 No. 2 Serial No. 438
潜水含水层非稳定井流降水的有限差分分析
刘 嘉，张功新，罗 彦，董志良
（中交四航工程研究院有限公司，广东 广州 510230）
摘要：滨海地区利用吹填淤泥或砂围海造陆将在地基上部形成一潜水含水层，但现有方法大都依赖于 Theis 非稳定承压
区外则仍按原含水层建立，计算所用参数见表 1， 其中假定各含水层为水平各向同性，竖向渗透系 数 kzz 取水平渗透系数 kxx 的 10%。模型边界条件 为：初始水头高度 h0=19.5 m，加固区四周设不排 水边界等效泥浆搅拌密封墙。
图 4 三维地下水渗流模型
图 5 地下水渗流模型平面图 （单位: m）
采用井点降水预压进行软基处理，其关键问 题在于降水设计。Dupuit 认为潜水井中进行定流 量（或定降深）抽水，经过一定时间 t，渗流将趋向 稳定，潜水面 由原来的水 平状态变为 漏斗状 [1]。 Muskat M 给 出 达 到 稳 定 的 时 间 t ≥0.5R2Sy/（kH） （Sy 为潜水含水层的给水度，R 为影响半径，k 为 土层渗透系数，H 为潜水含水层厚度）[2]。因此， 通过稳定井流抽水试验估计 k 值耗时过长。Theis C V 提出地下水径向非稳定承压井流计算公式[3]， 已知抽水量 Q，降深 s，经过抽水时间 t，可求得 土层压力传导系数 a，水力传导系数 T（=kM）和承
图 1 插设塑料排水板后的 s~t 曲线 图 2 插设塑料排水板后的降水漏斗时程曲线
·72·
水运工程
2010 年
图 3 插设塑料排水板后的 Neuman 拟合曲线
限差分公式：
CR i,j-
1
，k
(hi,j-1,k-hi,j,k)+CR i,j+
1
(hi,j,k-hi,j,k)+
，k
2
2
CC i,j-
计算过程分 8 个时间阶段：第 1~7 个阶段陆 续完成 24 个泵的安装，每安装完 1 个泵便开始抽 水，第 8 个阶段所有泵都开始抽水，整个过程共 历时 25 d。
第2期
刘 嘉，等：潜水含水层非稳定井流降水的有限差分分析
·73·
表 1 有限差分法计算参数一览
土层
kxx /
kyy /
kzz /
（cm·s-）1 （cm·s-）1 （cm·s-）1
（3）
式中符号意义见文献[19]。式（3）以矩阵形式可表
示为：
Ah=q
（4）
式中：A 为水头Biblioteka Baidu系数矩阵；h 为所求得水头矩
阵；q 表示各方程中所包含的所有常数项和已知
项，求解时根据 A 和 q 通过迭代获得 h。
2.4 有限差分法模型
如图 4 和图 5 所示，三维地下水渗流模型尺
寸为 300 m×245 m×20 m，加固区位于模型中部
2.3 有限差分法基本原理
在不考虑水的密度变化的条件下，孔隙介质
中地下水在三维空间的流动可以用下面的偏微分
方程来表示[19]：
! " ! " ! " 坠
坠x
kxx
坠h 坠x
+坠 坠y
kyy
坠h 坠y
+坠 坠z
kzz
坠h 坠z
-W =Ss
坠h 坠t
（2）
式中：kxx，kyy 和 kzz 分别为渗透系数在 x，y 和z 方
1 2
,j,k
(h-Ii,j,k)+CCj+
1 2
，j,k(hi+1,j,k-hi,j,k)+
CVi,j,k-
1 2
(hi,j,k-1-hi,j,k)+CV i,j,k+
1 2
(hi,j,k+1-hi,j,k)+
Pi,j,khi,j,k+Qi,j,k=Ssi,j,k(△rj△ci△υk)
△hi,j,k △t
Abstract: In the coastal area, an unconfined aquifer will be formed on the upper part of the foundation by
dredging and filling silts or sands to reclaim land, but the most existing methods are relying on the Theis, a well flow theory of non-steady confined aquifer to solve the hydro-geological parameters and conduct precipitation design that is lack of rationality. By applying the Neuman unsteady well flow theory to derive the hydro geological parameters of unconfined aquifer and then using the finite-difference method to analyze all derived parameters, the resulted values are consistent with the water level monitored results. The engineering examples show that this method is reasonable and reliable, and provides a way for the similar applications.
量 q 取实测值 48 m3/d，井内最大水位降深 15 m。
根据工程地质条件，将试验区含水层简化为自上
而下的 5 个土层：①0~6 m 第 1 个软软弱黏土层；
②6~8 m 砂层；③8~13 m 第 2 个超软弱黏土层；
④13~15 m 砂层；⑤15~20 m 硬黏土层，并将区内
视为单一潜水含水层，采用 Neuman 法求参结果，
Sy
加固区内等效土层 4.36×10-4 4.36×10-4 4.36×10-5 0.15
区外土层①，淤泥
1×10-7
1×10-7
1×10-8
0.02
区外土层②，砂层
1×10-2
1×10-2
1×10-3
0.26
区外土层③，淤泥
1×10-7
1×10-7
1×10-8
0.02
区外土层④，砂层
1×10-2
1×10-2
中图分类号：U 655.54+3
文献标志码：A
文章编号：1002-4972（2010）02-0070-05
Finite-difference analysis of unsteady well flow at unconfined aquifer
LIU Jia, ZHANG Gong-xin, LUO Yan, DONG Zhi-liang (CCCC Sihang Engineering Research Institute Co., Ltd., Guangzhou 510230, China)
收稿日期：2009-08-20 作者简介：刘嘉（1981—），男， 硕士，主要从事岩土工程数值分析方面的研究工作。
第2期
刘 嘉，等：潜水含水层非稳定井流降水的有限差分分析
·71·
目标值拟合井函数 W（u），形成了非稳定井流抽水 试验确定水文地质参数的解析方法[8]，较之图解法 和配线法更为简便。齐学斌按 Theis 基本方程，运 用 Gauss-Newton 迭代 法进行非线 性回归求参 [9]。 刘立才等借助梯度数学模型对 Theis 公式求解得出 T 和 S[10]。刘玉珍等基于 Theis 方程论述了灵敏度分 析法确定水文地质参数的基本原理，即把 s 视为 以 T 和 S 为变量的函数，利用最小二乘法求出理 论降深与观测降深最佳拟合时的 T 和 S[11]。巩彦文 通过对 W（u）的优化拟合，建立一个二元线性回归 方程求解回归方程中的待定常数，从而反算含水 层的水文地质参数[12]。类似的方法还见于 Swamee P K，Ojha C S P 和 Ajayi O 等人的研究中[13-17]。
压含水层给水度 S。其中，W（u）为 Theis 井函数，M 为承压含水层厚度，r 为井半径。Jacob C E 建议 用修正值 sc=s-s（/ 2H）描述潜水井流情况[4]，并要求 满足最大降深大于含水层厚度 5%[5]。
根据 Theis 解，由非稳定井流抽水试验所得 s~t 曲线 （s 为水位降深） 确定水文地质参数的主 要方法有图解法、配线法和解析法。前两种方法 人为误差大，求参过程复杂，适用性不大。梁文 康通过编写图解查表程序[6]，邹正盛[7]依据最小二 乘法建立自动配线法，都丰富了传统求解手段。 滕凯采用逐次逼近法，以标准剩余差最小为最优
1×10-8
0.26
区外土层⑤，黏土层 1×10-8
1×10-8
1×10-9
0.02
3 模型验证 如图 6 所示，两个中心点水位观测井的实测
值与有限差分法所得计算值趋势基本一致，在抽 水期的 0~6 d 计算水位高于实测水位，之后的时 间段则是预测值水位低于实测值。实测曲线在 0~ 5 d 之间没有水位观测点，所得的结果之间采用线 性插值，而计算曲线和实测曲线也在抽水后的第 6 d 相交，说明与观测结果相符，仅是实测数据偏 少未能用线性变化反映期间的过程。这一阶段潜 水泵正在逐个安装，区内抽水量较小，水位降低 速率应低于后期，预测曲线恰能更好地反映这点。 预测曲线的第 2 个阶段，地下水位下降速率增加， 而后慢慢变缓。这与非稳定井流理论相符，即认 为潜水面是可移动的物质界面，随着抽水时间的 增加，s~t 曲线斜率减小，有时甚至获得短时间的 “稳定”，这反映了潜水疏干排水的作用，如同含 水层得到“补给”，属于滞后给水现象。实测降深 曲线则是有高有低，一方面是测量过程中人为因 素的影响，另一方面潜水泵在抽水过程中需要定时 维护，也将使抽水量变小，水位出现短时间上升。 就总体趋势而言，预测水位降深略大于实测值。
向上的分量；h 为水头；W 为单位体积流量。离
散化后得到计算单元 （i，j，k） 的地下水渗流有
1 降水设计 1.1 工程地质条件
广州港南沙港区粮食及通用码头工程软基处 理试验区，加固面积 22 500 m（2 125 m×180 m），地 下水位离地表 0.5 m。根据地质勘察资料揭示，原 场地自上而下可分为：①淤泥混砂为主的近期人 工回填土，平均厚度 2 m；②淤泥或淤泥质黏土， 层①+②平均厚度 6 m，为第 1 个软土层；③混砂 层，平均厚度 2 m；④淤泥或淤泥质黏土，平均 厚度 5 m，为第 2 个软土层；⑤中粗或中细砂层， 平均厚度 2 m；⑥分布不均匀的粉质黏土、黏土、 砂层及硬黏土混合层。软基加固前在场地表层铺 设了 1.5 m 厚的砂垫层作为水平排水体及后续施 工垫层，并按 1.1 m 间距插设塑料排水板作为竖 向排水体，以及在加固区四周打设双排桩泥浆搅 拌墙作为止水帷幕。
2.2 水文地质求参
广义 Neuman 方程定义潜水含水层中的水位
降深 s 随时间 t 变化规律如下[18]：
s=
Q 4πt
W(uA，uB，β)
（1）
式中：W (uA，uB，β) 为潜水含水层井函数；uA=r2S/ （4Tt）（用于降水初 期 A 类型曲线），uB=r2Sy/(4Tt)
（用于降水后期 B 类型曲线），T=khH。抽水试验按
现有非稳定井流解析方法大都基于 Theis 系 统，其实质是承压井流，属二维情况，没有考虑 潜水井流存在垂直分流速，属于三维流动的实质。 并且从承压井中抽水是弹性贮量的释放，接近瞬 时完成，无压井则为含水层重力水疏干的释水机 制。针对软基处理工程中潜水井流降水的特点， 本文通过 Neuman 非稳定潜水完整井流求参，并 进行有限差分分析，研究数值计算和实测地下水 位降深规律，为今后相关设计和研究工作提供有 益的参考。
含水层井流理论求解水文地质参数，进行降水设计，缺乏合理性。结合广州港南沙港区粮食及通用码头软基试验区科研项
目，利用 Neuman 非稳定井流理论进行潜水含水层水文地质求参，将所得参数进行有限差分分析，数值结果与水位监测成果
一致。工程实例表明，该方法合理、可靠，为类似应用提供了一种新途径。
关键词：非稳定；潜水含水层；水文地质参数；降水；有限差分
照距离抽水井点 PW 一侧 5 m，10 m 和 15 m 远的
地方分别布设观测井 OW1，OW2 和 OW3。
根据图 1 和图 2 中的试验数据，借助商用软
件包 AquiferTest Pro 求参，可得插设塑料排水板
之后加固区内等效水平渗透系数 kh=4.36×10-4 cm/s，
Sy=0.15，拟合曲线见图 3。
（灰色区域），尺寸为 180 m×125 m×20 m，离散后
的计算单元 （i，j，k）尺寸为 2 m×2 m×2 m。加固区
内共打设 24 个抽水井 （PW1~PW24，间距 30 m×
30 m） 和 2 个水位观测孔 （OWZX1 和 OWZX2），
区外布置 8 个水位观测孔（OW1~OW8）。单井抽水
2010 年 2 月 第 2 期 总第 438 期
水运工程 Port & Waterway Engineering
Feb. 2010 No. 2 Serial No. 438
潜水含水层非稳定井流降水的有限差分分析
刘 嘉，张功新，罗 彦，董志良
（中交四航工程研究院有限公司，广东 广州 510230）
摘要：滨海地区利用吹填淤泥或砂围海造陆将在地基上部形成一潜水含水层，但现有方法大都依赖于 Theis 非稳定承压
区外则仍按原含水层建立，计算所用参数见表 1， 其中假定各含水层为水平各向同性，竖向渗透系 数 kzz 取水平渗透系数 kxx 的 10%。模型边界条件 为：初始水头高度 h0=19.5 m，加固区四周设不排 水边界等效泥浆搅拌密封墙。
图 4 三维地下水渗流模型
图 5 地下水渗流模型平面图 （单位: m）