数值模拟

题目：一承压含水层，平面分布如图，地表高程为20m，隔水顶板底面高程为-20m，厚度M为20m，渗透系数K为10m/d。2008年9月18日测得水位埋深为10m。南北边界为隔水边界，西部边界水位埋深保持15m，东部边界水位埋深为10m。含水层平面区域东西长2000m，南北长1500m。请（1）建立该含水层的数学模型；（2）预测含水层中部以2000方/天开采地下水，含水层水头分布动态。

解：1）取坐标系如图，水头基准面为海平面，基准面高程为0。水头即为高程。

H2H2HeTxxTyy （1） 22txyH(x,y,0)10mH|x015m,H|x200010m （2）

HH|y00,|y15000 （3） yy

2）取x=400m，y=300m，剖分网格如上图。设开采井孔在含水层中心位置如图。设e=0.02，时间步长t=1d。

对于每个剖分网格，根据（1）式，使用二维显式差分法变换得

1nnnnnhin,jxhi1,jxhi1,jyhi,j1yhi,j1(12x2y)hi,jte （4）

（i，j = 1，2，… ，4）

其中xKMt10*20*10.063， 22e(x)0.02*400 yKMt10*20*10.111 22e(y)0.02*300 对于一般网格，有=0。对于有井孔的中心网格， 即

20000.017，

400*300te1*0.0170.833

0.02 将上述值带入公式（4）得

hi,j0.063hi1,j0.063hi1,j0.111hi,j10.111hi,j10.652hi,j（一般网格） hi,j0.063hi1,j0.063hi1,j0.111hi,j10.111hi,j10.652hi,j0.833（中

心网格）

根据以上公式，结合题目所给边界条件，利用EXCEL软件算得 t=0d 时， 各格点水头情况如下：

n1nnnnnn1nnnnn15 15 15 15 15

t=2d

10.315 10.315 10.315 10.315 10.315 10.59031 10.59031 10.59031 10.59031 10.59031 11.23482 11.23133 11.21657 11.23133 11.23482 11.90082 11.88008 11.84133 11.88008 11.90082

10 10 10 10 10 10 10 9.167 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 9.999516 9.996033 9.981271 9.996033 9.999516 9.974235 9.953496 9.914748 9.953496 9.974235

10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

15 15 15 15 15

t=5d

10.01985 10 10 10.01985 9.907537 10 9.967366 8.623884 9.947521 10.01985 9.907537 10 10.01985 10 10 10.13972 10.08814 9.891443 10.08814 10.13972 10.33392 10.19636 9.905021 10.19636 10.33392

9.895251 9.552422 7.747118 9.552422 9.895251 9.617304 9.126146 7.137881 9.126146 9.617304

9.986357 9.934774 9.738081 9.934774 9.986357 9.856042 9.731295 9.439958 9.731295 9.856042

15 15 15 15 15

t=10d

15 15 15 15 15

t=15d

15 15 15

12.26578 10.43874 9.40274 9.675463 9.909505 12.23528 10.28958 8.878008 9.527234 9.879356 12.18643 9.973078 6.855023 9.211438 9.830784

15 15

t=20d

15 15 15 15 15

12.23528 10.28958 8.878008 9.527234 9.879356 12.26578 10.43874 9.40274 9.675463 9.909505 12.46889 12.43571 12.384 12.43571 12.46889

10.50869 10.3542 10.03204 10.3542 10.50869

9.256081 8.723901 6.693265 8.723901 9.256081

9.508604 9.354232 9.032199 9.354232 9.508604

9.827826 9.794731 9.743108 9.794731 9.827826

10 10 10 10 10 10 10

从以上计算结果可以看出各节点水头动态，并可以得出以下结论：

1、 从左向右，第一二列网格水头值渐渐上升，第三四五列网格水头值则逐渐下降。 2、 以布有抽水井孔的中心网格为中心，形成一个水位下降漏斗，逐渐向四周辐射。 3、 t≥15d 后，各节点水头值渐渐趋于稳定。

注：采用反映法处理隔水边界。