地下水建模方法和步骤

稳定性：如果在求解差分方程过程中，某时间步引入 某个误差，而在以后的各时段计算中，该误差不再扩 大，则称该差分格式是稳定的。
一维显示格式的收敛条件和稳定条件是： 0 1
2
(6)算例：显式有限差格式
应用实例：河间地块承压水流模型
设两条河流平行、完全切割含水层，含水层等厚、均质 各向同性。
步骤：
0
20
及其误差。
已知泰勒公式
f(x 0 x)f(x 0)f'(x 0) xf'2 '(!)( x)2 f(x 0 x)f(x 0)f'(x 0) xf'2 '(!)( x)2
方法一
A
B
① 由A得：
f(x0)f(x0 xx )f(x0)O ( x)
称
f(x0
x)f(x0) x
为f(x)在x0处的一阶前向差商，O(x)
KMt
e (x)2
将该式子代入得到：
K M h in 1 (2 h x in )2 h in 1(xi,tn)eh in 1 th in
h in 1h in 1 (1 2)h inh in 1 t (xi,tn) e
（i=1，2，3，.....N-1），（n=1，2，3，.....）
( x ) 2
t
H in 1 2 H in H in 1H in 1 H in
( x )2
t
i 2 ,3 ,,n x 1
隐式差分格式
H ( x 0 x , t 0 t ) 2 H ( x 0 , t 0 t ) H ( x 0 x , t 0 t ) H ( x 0 , t 0 t ) H ( x 0 , t 0 )
n 1 2
( x)2
H
n 1 3
H
n 1 2
t
H
n 2
H
n 1 2
2
H
n 1 3
( x)2
H
n 1 4
H
n 1 3
t
H
n 3
方
程
组
H
n 1 ni 3
2H
n 1 ni 2
( x )2
H n1 ni 1
H
n 1 ni 2
H
n ni 2
t
H
n 1 ni 2
2H
n 1 ni 1
1.1有限差分法
(1)有限差分法原理 (2)两种方法建立有限差分方程 (3)求解有限差分方程 (4)收敛性和稳定性概念 (5)算例
(1)有限差分法的基本原理
将连续的问题离散后求解：
➢ 方法一．以地下水流基本微分 方程及其定解条件为基础， 在 渗流区剖分基础上，用差商代 替微商，将地下水流微分方程 的求解转化为差分方程（代数 方程）求解。
v |x (x,y,z,t)
v |x (xx,y,z,t)
x方向流入 (vx)|(x,y,z,t) yzt
x方向流出 (vx)|(xx,y,z,t) yzt
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
基于达西定律，x，y，z方向流入—流出分别为:
A [ v x |x v x |x x ] y z t [ H i 1 ,j , k x H i ,j , k H i 1 ,j , k x H i ,j , k ] z x y t
D 源汇项 zxyt
t时段内，侧向流入与源汇项导致六面体水量变化量为： A+B+C+D
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
六面体内地下水储存量的变化为
s x y z H s x y z H ( x , y , z , t t ) H ( x , y , z , t )
③由A-B可以得：
f'(x 0 )f(x 0 x 2 ) x f(x 0 x ) O ( x )2
称 f(x0x)f(x0x)
2x
为f(x)在x0处的一阶中心差商，
O(x)2 为截断误差。
④由A+B可以得：
f''(x 0 ) f(x 0 x ) 2 ( f( x x )0 2 ) f(x 0 x ) O ( x )2
（1）基础资料的分析 （2）概念模型 （3）数学模型 （4）数值方法及计算机程序 （5）参数 （6）结果分析
建立数学模型
（1）模型概化
由所述水文地质条件，可以概化为一维承压水流问题。
（2）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内， 从而量化各变量的取值范围。本例，取x-轴原点位于左端河，右侧 为正向，设两河流间距为L.
称 f(x0x)2( f(xx )0 2)f(x0x)为f(x)在x0处的二阶中心差商，O(x)2 为截断误差。
方法一
(2)有限差分方程建立(续)
对于偏导数（偏微商），类似可以得到相应的差商：
H (x 0,t0)H (x 0,t0 t) H (x 0,t0)
t
t
H (x0,t0)H (x0 x ,t0) H (x0,t0)
2）建立差分方程：在网格系统中任意取一点 (xi ,tn )
记为（i，n） 设 H(x,t) 是问题的解，则在 (xi , tn ) 处有
e
Hn t i
KM2xH 2 n i
(xi,tn)
显式格式（续1）
用差商代替微商：eH t in
KM2xH 2 n i
(xi,tn)
Hn Hin1Hin O(t)
（3）数学模型
e Ht KM2xH2
0xL,t0
H t0
H0(x)
0xL
t 0 H x 01 (t)H , X L2 (t)
差分方程及其解法—显式格式
1)网格剖分：
①将（0—L）分成 N 等份， x l
记 xi ix，(i=0,1,2,3，4……N N)
②取时间步长 t，记 tn nt （n=0、1、2、3、4……）
显式格式（续3）
3）显示差分方程的求解
计算各结点初始时刻水头值
H t0
H0(x)
利用差分方程计算各结点t1时刻水头值
利用边界条件计算边界结点水头值
重复2、3步，直到计算出拟计算的各个时刻的水头 值
算例（续4）
e Ht KM2xH2
H t0
H0(x)
H x 01 (t)H , X L2 (t)
地下水建模方法和步骤
中国地质大学(武汉)环境学院 2012.8
地下水建模方法和步骤
➢ 1.求解地下水运动方程的数值方法 ➢ 2.地下水数值模型建模步骤 ➢ 3.建模所需要的基本资料
1.数值方法
绝大部分数学模型是无法用解析法求解的，数 值化就是将数学模型转化为可解的数值模型。
✓ 1.1有限差分法 ✓ 1.2有限单元法 ✓ 1.3分有限差分法 ✓ 1.4半解析半数值法 ✓ 1.5边界元法
x
x
2 H ( x x 2 0 ,t0 ) H (x 0 x ,t0 ) 2 H ( (x x 0 ) ,2 t0 ) H (x 0 x ,t0 )
一维控制方程差分格式
方法一
控制方程
T2h(xx2,t)
h(x,t)
t
网格剖分nx个
显式差分格式
H ( x 0 x ,t 0 ) 2 H ( x 0 ,t 0 ) H ( x 0 x ,t 0 ) H ( x 0 ,t 0 t ) H ( x 0 ,t 0 )
10mx,0 在上述模型中，设L=1000米 H0(x) 20mx,0
e 0.00K 4 8 m M 2 /m d 0 1 ( t ) ,, 2 m 2 0 ( t ) 1 ,m 0
取空间步长为200米，时间步长为0.25天，分别计算各节点各时刻的 水头值。
算例（续5）
h i n 1 h i n 1 ( 1 2 ) h i n h i n 1 t/ei,n
(2)有限差分方程建立
方法一：差商代替微商
导 导数的定义
数 的
f(x0) lx i0m f(x0 xx )f(x0)
有 当 x 非常小的时候，有
限 差
f(x0)f(x0 xx)f(x0)
商 上式右端项即为f(x)在x0处的差商。 近 这样定义的差商很容易理解，但不知道用差商代 似 替微商所产生的误差。下面利用泰勒公式导出差商
由水均衡原理得三维地下水流动方程的有限差分格式
[ H i1, j,k H i, j,k H i1, j,k H i, j,k ]zxyt
x
x
[ H i, j1,k H i, j,k H i, j1,k H i, j,k ]zxyt
y
y
[ H i, j,k 1 H i, j,k
H H i, j,k 1
为截断误差。
② 由B 得：
f(x0)f(x0) fx (x0 x)O ( x)
称
f(x0)f(x0x) x
为f(x)在x0处的一阶后向差商，O(x)
为截断误差。
方法一 f ( x 0 x ) f ( x 0 ) f '( x 0 ) x f ''2 ( ! x 0 ) ( x ) 2 f '''3 ( ! x 0 ) ( x ) 3 f ( 4 4 ) ! () ( x ) 4 A f ( x 0 x ) f ( x 0 ) f '( x 0 ) x f ''2 ( ! x 0 ) ( x ) 2 f '''3 ( ! x 0 ) ( x ) 3 f ( 4 4 ) ! () ( x ) 4 B
B [ v y |y v y |y y ] x z t [ H i,j 1 ,k y H i,j,k H i,j 1 ,k y H i,j,k ] z x y t C [ v z |z v z |z z ] x y t [ H i ,j ,k 1 z H i ,j ,k H i ,j ,k 1 z H i ,j ,k ] z x y t
➢ 方法二．在渗流区剖分的基础 上，直接由达西定律和水均衡 原理，建立各个均衡区的水均 衡方程，即差分方程。
矩形网格 多边形网格
网格划分的基本类型
（1）先划格线，格 点位于网格中心
均 衡 网
格
（2）先规定格点位
置，再垂直平分两相
节
邻结点的连线作格线， 点
形成的网格即为水均
网 格
衡区
MODFLOW网格系统
( x ) 2
t
H in 1 1 2 ( H x in ) 2 1 H in 1 1 H in 1 tH in i 2 ,3 ,,n 1 x
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
取右图所示得微小六面 体。设与x,y,z,方向对应得 主渗透系数分别为Kx, Ky,Kz；建立均衡期t时 段内，微小均衡六面体的 水量守恒方程。
t i
t
2xH 2 nHin1(2H x)in2Hin1O(x)2 i
将上述两式舍去余项，代入方程并记
H
n i
为
h
n i
得到
K M h in 1 (2 h x in )2 h in 1(xi,tn)eh in 1 th in
显然该式具有截断误差 O(t)O(x)2
显式格式（续2）
引入无量纲变量：
n 4
H
n 1 3
t
H
n 3
H
n ni
3
2
H
n ni 2
( x )2
Hn ni 1
H
n 1 ni 2
H
n ni 2
t
H
n ni 2
2
H
n ni 1
( x )2
H
n ni
H
n 1 ni 1
H
n ni 1
t
(3)差分方程求解
一维隐式差分格式 网格个数为ni
H
n 1 1
2H
( x)2
H n1 ni
H
n 1 ni 1
H
n ni 1
t
迭代求解
MODFLOW
PCG SIP SOR WHS SAMG GMG
(4)差分方程的收敛性和稳定性
截断误差：用差商代替微商时，地下水流动方程产生
的误差为截断误差。O(t)O(x)2
收敛性：当空间步长和时间步长趋于0时，有限差分方 程的精确解趋于地下水流动问题微分方程定解问题的 精确解。则称该差分格式是收敛的。
K e( x M )t208 .0 *20 ** (0 2 0 4 .20 )25 01/4
h i n 1 1 4 h i n 1 ( 1 2 1 4 ) h i n 1 4 h i n 1 1 4 ( h i n 1 2 h i n h i n 1 )
Time/day x=0 m x=200 m x=400 m x=600 m x=800 m x=1000 m
i, j,k
]zxyt
z
z
zxyt
s xyz
H n1 i, j,k
H
n i,
j,k
有限差分法:三维（MODFLOW）
差商代替微商
(3)差分方程求解
一维显式差分格式
网格个数为ni 直接求解
H
n 1
2H ( x
n 2
)2
HLeabharlann Baidu
n 3
H
n 1 2
t
H
n 2
H
n 2
2H ( x
n 3
)2
H