地下水数值模拟03

差分公式对比
名称
前差
公式
df f ( x x) f ( x) dx x
截断误差
O(x) O(x)
一 阶 导 数
后差 中心差
df f ( x) f ( x x) dx x
df f ( x x) f ( x x) dx 2x
O(x )
2
二阶导数
该问题属于一维承压非稳定流的定解问题。求解该问题，需要对空间 和时间进行离散，形成的差分网格称为时空网格。
网格剖分
• • •
以等距剖分为例 将研究区域[0，L]用直线等分为n份，把时间段[0,T]用直线等分成m份 以 H ik 表示结点（i,k）处的水头
2H * H T 2 x t H ( x, t ) H 0 ( x ) t 0 H ( x, t ) (t ) 0 x 0 H ( x, t ) x L L (t )
开始
输入初始参数 t=0 对i=0,1,…,n循环执行
hi0 H0 ix
t=t+Δt 对i=1,…,n-1循环执行
hik 1 hik1 (1 2 )hik hik1
• 用差分方程的解作为微分方程的近似解。
一、有限差分法的基本思想
• 2、求解步骤
选取网格； 对微分方程及定解条件选择差分近似，列出
差分格式；
求解差分格式； 讨论差分格式解对于微分方程解的收敛性及 误差估计。
差分网格
结点 格点
y
x
Δx，Δy——空间步长 Δt——时间步长
一、有限差分法的基本思想
显式差分方程的求解
差分方程与离散后的定界条件构成了数学模型 的显式差分方程问题。其求解步骤如下： 由初始条件给出t0时刻各结点的水头值h00， h10，…hl0；再根据差分方程，在k=0时， 分别取i=1，i=2,…i=n-1,便可求得t1时刻各 内结点的水头值h11，…hn-10。 由以上计算的h11， h21 …hn-11值及由边界 条件计算的h01和hl1,再次利用差分公式 （取k=1， i=1，i=2,…i=n-1，便可计算得 t2时刻各结点的水头值。如此重复，便可 计算出t3 ，t4 ，…各时刻的水头分布值。
二、导数的有限差分近似表示
• 2、导数的有限差分形式
用泰勒级数展开可以推导出导数的有限差分形式。
df (x) 2 d 2 f (x)3 d 3 f f ( x x) f ( x) x 2 3 dx 2! dx 3! dx
df (x) 2 d 2 f (x)3 d 3 f f ( x x) f ( x) x 2 3 dx 2! dx 3! dx
2 1 h1 Wx 2 / T 0 1 2 1 h 2 Wx / T 2 2 h 1 2 1 W x / T n2 2 1 2 hn 1 Wx / x f x f x 2x f x x f x x 4x 2 2x 2x x f x 2x f x x f x x f x 2x 1 4 f? x 3 x 2 2 f x 2x 2 f x x 2 f x x f x 2x 2x 3
整理得：
Tt k k k k 1 k h 2 h h h h i 1 i i 1 i i * (x) 2


Tt 定义 * (x) 2
hik 1 hik1 (1 2 )hik hik1
只要知道了k时段初始时刻tk各结点的hik值，便可计算出k时段末了时刻tk+1的 hik+1值，各方程可独立求解，因此，这种方程称为显式有限差分方程。
思考与练习
• 高阶导数的有限差分形式
2 f ? xy
2 f xy 0
f f y y 1 3 2h 1 f5 f 6 f8 f 7 2h 2h 2h 1 2 f5 f 7 f 6 f8 4h
10 8 11 3 7 2 0 4 12 5 1 6 9
4 f x 4
0
2 f x 2
2 f 2 f x 2 2 2 1 3 x h2
0
f1 f 3 2 f 0 1 f 9 f 0 2 f1 f 0 f11 2 f 3 2 h2 h2 h2 h2 1 4 6 f 0 4 f1 4 f 3 f 9 f11 h
三、简单水文地质模型的有限差分方程组
• 3、有限差分方程
（2）建立有限差分方程
先任取一结点i进行分析。
2H T W 0 x 2
hi 1 2hi hi 1 T W 0 2 (x)
W x 2 T
n-1个线性代 数方程，未知 量共n-1个 故方程可解。
移项整理，得： hi 1 2hi hi 1
( 0 x L ,0 t T ) (0 x L ) (0 t T ) (0 t T )
导数可以利用一阶、二阶导数的差商代替，由于一阶导数 可以有三种差商表示，因此分别对水头关于时间的导数项分 别运用前差、后差、中心差将得到三种差分格式。
显式有限差分
隐式有限差分
←
←
前差
后差 中心差
中心式有限差分 ←
一、一维显式有限差分格式
2H H T * x 2 t
向前差分
k 1 hik hik1 2hik hik1 * hi = T 2 t (x)
(i-1,k) (i,k+1) (i,k) (i+1,k)
截断误差为：O([Δx]2)+O(Δt)
d2 f f ( x x) 2 f ( x) f ( x x) 2 O ( x ) 2 2 dx (x)
思考与练习
• 高阶导数的有限差分形式
d 3 f d 3 fd d 2 f d f x x f x x 2 f x 3 2 2 ? dx dx dx x 3 dx dx 1 df x x df x x df x f 2 2 2 x dx dx dx ?
线性代数方程组！
最有效的求 解方法 ——追赶法
系数矩阵中的元素仅位于三条对角线上，系数矩阵对称且 正定，故称为三对角线性代数方程组。
第二节
几种主要的差分格式
水文地质模型描述
以一维河间地块承压含水层 中的水流问题为例。 含水层均质各向同性，不考 虑垂向补给，两河流边界的 水位随时间变化，分别φ0(t)、 φL(t) 。 试研究含水层的水头分布。
三、简单水文地质模型的有限差分方程组
• 1、水文地质模型
以一维河间地块承压含水层中的水流问题为例。 考虑一个以通过x＝0和x＝L处的长且直的河流为界的承压含水层，该含水 层均质各向同性，顶底板水平，上覆弱透水层，垂向补给强度为W，两河流 边界的水位分别φ0、 φL为且不随时间变化。 试研究含水层的水头分布。
——dy、dx分别是函数及自变量的微分, ——
dy 是函数对自变量的导数，又称微商。 dx
—— y 、x 分别称为函数及其自变量的差分 y —— x为函数对自变量的差商。 由导数（微商）和差商的定义可知， 当自变量的差分（增量）趋近于零时， 就可以由差商得到导数。因此在数值 计算中常用差商近似代替导数。
2H * H T 2 x t H ( x, t ) H 0 ( x ) t 0 H ( x, t ) (t ) 0 x 0 H ( x, t ) x L L (t ) ( 0 x L ,0 t T ) (0 x L ) (0 t T ) (0 t T )
对于所有内结点1、 2、…、n-1，建立 结点相应的差分方 程组
Wx 2 0 2h1 h2 T Wx 2 h1 2h2 h3 T 2 W x hn 3 2hn 2 hn 1 T Wx 2 hn 2 2hn 1 L T
三、简单水文地质模型的有限差分方程组
• 3、有限差分方程
（3）求解
将方程表示成矩阵形式，则有：
Wx 2 2 h h 0 2 1 T Wx 2 h 2 h h 2 3 1 T Wx 2 hn 3 2hn 2 hn 1 T Wx 2 hn 2 2hn 1 L T
第三章
有限差分法
主要内容
• • • • 有限差分法的基本原理 几种主要的差分格式 二维渗流问题的差分方程 一般差分方程组的解法
第一节
有限差分法的
基本原理
一、有限差分法的基本思想
• 1、基本原理
• 从物理现象引出相应微分方程（方程+边界条件）； • 用差分网格离散求解域；
• 用差分公式将基本方程转化为差分方程（代数方程）；
三、简单水文地质模型的有限差分方程组
• 3、有限差分方程
（1）网格剖分
2H W 0 T 2 x H ( x) x 0 0 H ( x) x L L
(0 x L )
沿河流的方向取单宽[0,L]作为研究区域， L 将L等分为n份，空间步长 x n 对每个结点进行编号，结点编号由左向右 依次为0，1，2，…，i，…，n。共有n+1个 结点，其中2个边界结点，n-1个内结点。
• 3、优缺点
优 点 缺 点
1.简单问题的数学表达 式和计算的执行过程 比较直观、易懂； 2.算法效率比较高，易 编程；
1、对自然边界处理的灵 活性较差。 2、对溶质运移等问题， 精度受限
二、导数的有限差分近似表示
• 1、差分的概念
设有x的解析函数y(x)，函数y对x的导数为：
dy y f ( x x) f ( x) lim lim x 0 dx x x 0 x
当x 0时，C2 0
当x L时， W 2 L C1 L 0 L 2T
L 0 W L C1 L 2T C2 0
H x
W 2 L 0 WL x x 0 2T L 2 T
三、简单水文地质模型的有限差分方程组
• 2、数学模型
2H (0 x L ) T 2 W 0 x H ( x) x 0 0 H ( x) x L L 2H W 解析解 x 2 T H W xC x T W 2 H x C1 x C2 2T