地下水建模方法和步骤

1.1有限差分法
(1)有限差分法原理
(2)两种方法建立有限差分方程
(3)求解有限差分方程
(4)收敛性和稳定性概念 (5)算例
(1)有限差分法的基本原理
将连续的问题离散后求解：
方法一．以地下水流基本微分 方程及其定解条件为基础， 在 渗流区剖分基础上，用差商代 替微商，将地下水流微分方程 的求解转化为差分方程（代数 方程）求解。 方法二．在渗流区剖分的基础 上，直接由达西定律和水均衡 原理，建立各个均衡区的水均 衡方程，即差分方程。
网格个数为ni
直接求解
(3)差分方程求解

一维隐式差分格式
网格个数为ni
n 1 n 1 n H1n 1 2 H 2 H 3n 1 H 2 H2 2 ( x ) t n 1 n 1 H2 2 H 3n 1 H 4 H 3n 1 H 3n 2 ( x ) t 方 程 组 n 1 n 1 n 1 n 1 n H ni 2 H H H H 3 ni 2 ni 1 ni 2 ni 2 2 (x) t n 1 n 1 n 1 n 1 n H ni 2 H H H H ni 1 ni ni 1 ni 1 2 2 (x) t
n i
2H KM ( xi , tn ) x 2 i
n
n
H Hin1 Hin O(t ) t i t
''
f ( x0 x) 2 f ( x0 ) f ( x0 x) 2 2 称 为 f ( x ) 在 x 处的二阶 中心 差商， 为截断误差。 O ( x ) 0 (x)
方法一
(2)有限差分方程建立(续)

对于偏导数（偏微商），类似可以得到相应的差商：
H ( x0 , t0 ) H ( x0 , t0 t ) H ( x0 , t0 ) t t
1 一维显示格式的收敛条件和稳定条件是： 0 2
(6)算例：显式有限差格式
应用实例：河间地块承压水流模型
设两条河流平行、完全切割含水层，含水层等厚、均质 各向同性。
步骤：
（1）基础资料的分析
（2）概念模型 （3）数学模型 （4）数值方法及计算机程序 （5）参数 （6）结果分析
建立数学模型
由水均衡原理得三维地下水流动方程的有限差分格式
[ H i 1, j ,k H i , j ,k x H i 1, j ,k H i , j ,k x ]zxyt ]zxyt ]zxyt
[ [
H i , j 1,k H i , j ,k y H i , j ,k 1 H i , j ,k
控制方程
2 h ( x, t ) h( x, t ) T x 2 t
方法一
网格剖分nx个
显式差分格式
H ( x0 x, t 0 ) 2H ( x0 , t 0 ) H ( x0 x, t 0 ) H ( x0 , t 0 t ) H ( x0 , t 0 ) 2 t (x)
②取时间步长
2）建立差分方程：在网格系统中任意取一点 ( xi , t n ) 记为（i，n） 设
H ( x, t ) 是问题的解，则在 ( xi , t n ) 处有
n
H e t
i
H KM ( xi , tn ) 2 x i
2
n
显式格式（续1）
用差商代替微商：
H e t
1 n 1 n 1 n 1 n H in 2 H H H H 1 i i 1 i i (x) 2 t
i 2,3, , nx 1
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
取右图所示得微小六面 体。设与x,y,z,方向对应得 主渗透系数分别为Kx, Ky,Kz；建立均衡期t时 段内，微小均衡六面体的 水量守恒方程。
H i 1, j ,k H i , j ,k x H i 1, j ,k H i , j ,k x ]zxyt
H i , j 1,k H i , j ,k y

H i , j 1,k H i , j ,k y
]zxyt
H i , j , k 1 H i , j , k z
H in1 2H in H in1 H in1 H in (x) 2 t i 2,3, , nx 1
隐式差分格式
H ( x0 x, t 0 t ) 2H ( x0 , t 0 t ) H ( x0 x, t 0 t ) H ( x0 , t 0 t ) H ( x0 , t 0 ) t (x) 2
地下水建模方法和步骤
中国地质大学(武汉)环境学院 2012.8
地下水建模方法和步骤
1.求解地下水运动方程的数值方法
2.地下水数值模型建模步骤 3.建模所需要的基本资料
1.数值方法
绝大部分数学模型是无法用解析法求解的，数 值化就是将数学模型转化为可解的数值模型。 1.1有限差分法 1.2有限单元法 1.3分有限差分法 1.4半解析半数值法 1.5边界元法
H i , j 1,k H i , j ,k y H i , j ,k 1 H i , j ,k
z z 1 n zxyt s xyz H in H , j ,k i , j ,k


有限差分法:三维（MODFLOW）
差商代替微商
(3)差分方程求解
MODFLOW
迭代求解
PCG SIP SOR WHS SAMG GMG
(4)差分方程的收敛性和稳定性



截断误差：用差商代替微商时，地下水流动方程产生 2 的误差为截断误差。O(t ) O(x) 收敛性：当空间步长和时间步长趋于0时，有限差分方 程的精确解趋于地下水流动问题微分方程定解问题的 精确解。则称该差分格式是收敛的。 稳定性：如果在求解差分方程过程中，某时间步引入 某个误差，而在以后的各时段计算中，该误差不再扩 大，则称该差分格式是稳定的。

一维显式差分格式
n n 1 n H1n 2 H 2 H 3n H 2 H2 2 (x) t n n H2 2 H 3n H 4 H 3n 1 H 3n 2 ( x ) t n n n n 1 n H ni H ni 3 2 H ni 2 H ni 1 2 H ni 2 2 (x) t n n n n 1 n H ni 2 H H H H ni 1 ni ni 1 ni 1 2 2 ( x ) t
已知泰勒公式
f '' ( ) f ( x0 x ) f ( x0 ) f ( x0 )x ( x ) 2 2! f '' ( ) ' f ( x0 x ) f ( x0 ) f ( x0 )x ( x ) 2 2!
'
方法一
A B
① 由A得：
f ( x0 x ) f ( x0 ) f ( x0 ) O ( x ) x
'
③由A-B可以得：
f ( x0 x ) f ( x0 x ) f ( x0 ) O ( x ) 2 2x
'
称
f ( x0 x) f ( x0 x) 为f(x)在x0处的一阶中心差商， O(x)2 为截断误差。 2x
④由A+B可以得：
f ( x0 x ) 2 f ( x0 ) f ( x0 x ) 2 f ( x0 ) O ( x ) ( x )2
矩形网格
多边形网格
网格划分的基本类型
（1）先划格线，格 点位于网格中心 （2）先规定格点位 置，再垂直平分两相 邻结点的连线作格线， 形成的网格即为水均 衡区
均 衡 网 格
节 点 网 格
MODFLOW网格系统
(2)有限差分方程建立
方法一：差商代替微商
导 数 的 有 限 差 商 近 似
导数的定义

H i , j , k 1 H i , j , k z
]zxyt
源汇项 zxyt
t时段内，侧向流入与源汇项导致六面体水量变化量为：
A+B+C+D
Leabharlann Baidu
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
六面体内地下水储存量的变化为
s xyzH s xyzH ( x, y, z, t t ) H ( x, y, z, t )
0 x L
H t 0 H 0 ( x)
H
x0
1 (t ), H
X L
2 (t )
t0
差分方程及其解法—显式格式
1)网格剖分：
①将（0—L）分成 N 等份，
记
xi ix，(i=0,1,2,3，4……N)
t，记 t n nt （n=0、1、2、3、4……）
l x N
vx |( x, y, z ,t )
vx |( xx, y,z,t )
x方向流入 (vx ) |( x , y , z ,t ) yzt
x方向流出 (vx ) |( x x , y , z ,t ) yzt
(2)有限差分方程建立(续)
方法二：达西定律和水均衡原理
基于达西定律，x，y，z方向流入—流出分别为: A [vx | x vx | x x ]yzt [ B [v y | y v y | y y ]xzt [ C [vz |z vz |z z ]xyt [ D
f ( x0 ) lim
x 0
f ( x0 x) f ( x0 ) x
当 x 非常小的时候，有
f ( x0 ) f ( x0 x ) f ( x0 ) x
上式右端项即为f(x)在x0处的差商。
这样定义的差商很容易理解，但不知道用差商代 替微商所产生的误差。下面利用泰勒公式导出差商 及其误差。
（1）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。 （2）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内， 从而量化各变量的取值范围。本例，取x-轴原点位于左端河，右侧 为正向，设两河流间距为L. （3）数学模型
H 2 H e KM t x 2
0 x L, t 0
f '' ( x0 ) f ''' ( x0 ) f (4) ( ) 2 3 f ( x0 x ) f ( x0 ) f ( x0 )x ( x ) ( x ) ( x ) 4 2! 3! 4!
'
A方法一
f '' ( x0 ) f ''' ( x0 ) f (4) ( ) 2 3 f ( x0 x ) f ( x0 ) f ( x0 )x ( x ) ( x ) ( x ) 4B 2! 3! 4!
称
f ( x0 x) f ( x0 ) 为f(x)在x0处的一阶前向差商，O(x) 为截断误差。 x
② 由B 得：
f ( x0 )
称
f ( x0 ) f ( x0 x ) O ( x ) x
f ( x0 ) f ( x0 x) 为f(x)在x0处的一阶后向差商，O(x) 为截断误差。 x
H ( x0 , t 0 ) H ( x0 x, t 0 ) H ( x0 , t 0 ) x x
2 H ( x0 , t0 ) H ( x0 x, t0 ) 2H ( x0 , t0 ) H ( x0 x, t0 ) 2 x (x) 2
一维控制方程差分格式