![地下水数值模拟基础](https://wfqqreader-1252317822.image.myqcloud.com/cover/292/40936292/b_40936292.jpg)
二、几种主要的差分格式
(一)显式差分格式
为了便于说明,以下列一维河间地块均质各向同性承压含水层中的地下水流问题
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_10.jpg?sign=1739373091-6dLFc6mMaDhMaV9YH9cXAXgihRwgjUNm-0-67c4d7f1fc08a921c46ae4c6e9d721c2)
为例来加以说明。
首先将研究区域[0,L]用直线等分为l份,步长Δx=L/l,把时间段[0,Tsum]用直线等分成m份,时间步长Δt=Tsum/m,构成如图2-1所示的网格,结点坐标xi=iΔx,tκ=κΔt(i=0,1,…,l;κ=0,1,…,m),简记为(i,κ),并以表示H(iΔx,κΔt),以
表示原方程的差分方程解(即H的近似值)。
式(2-5)中的导数,用差商代替,在典型结点(i,κ)处表示为
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_13.jpg?sign=1739373091-AHpjPJIPIJ2e1M7JjlekWbDJ3DZKITgl-0-4850fa6ae39eac45f969bb71bcc38920)
图2-1 研究区域网格示意图
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_14.jpg?sign=1739373091-wbWVFFpM9ixW35UMhFmSVyOFRtWth8oT-0-76246fcead8cf26ec735b6647713327b)
略去O(Δt)和O([Δx]2),可得和式(2-5)以及图2-1中x,t平面上的网格对应的差分方程为
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_15.jpg?sign=1739373091-cRRTk4MsVFILHJqKTyeX9stSaYQuFxL4-0-ee9b740a0a2943ea7450623f9bffa647)
截断误差为O(Δt)+O([Δx]2)。即当Δt和Δx很小时,这一误差是Δt阶的一个量与(Δx)2阶的一个量之和。若定义
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_16.jpg?sign=1739373091-hfMOjnZpTdctiDDrBF6PGxrSs7sa6UXs-0-1b04086dc69a95be207e4221385f2853)
则式(2-10)可改写为
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_17.jpg?sign=1739373091-tuTaC9UTk3Cephh5aWhWwTKzwMyMJJ27-0-de5588976e892e2b0722389746c53c91)
由此可知,只要知道某一时段κ开始时刻tκ各结点的值,利用上式便能算出tκ+1时刻,即κ时段终了时刻的
值(1≤i≤l-1,1≤κ≤m)。所以称这一方法为显式方法。边界结点的水头值则由边界条件给出,即
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_20.jpg?sign=1739373091-IvnT6kaqRv6tpStgTMtfG0AnPCzChBX6-0-c329aa2824a9c3672a0d50d13250be2a)
这样利用t=0时刻时各结点的值已由初始条件式(2-6)
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_22.jpg?sign=1739373091-FkZF9Gcy8IaFCtl3K2IsmIWBGklouLHq-0-d36e861bb8b1764a2d3aa0803a0214f1)
给出的情况,直接计算t1时刻各内部结点的水头值,并利用边界条件补充边界结点上的水头值
。再把求得的t1时刻各结点的水头作为初值,重复上述过程求t2时刻各结点的水头值。如此一个时间水平、一个时间水平地做下去,就能求得计算区Ω上所有时刻的水头值。
前面我们给出了求的方法,但必须回答一个问题,即差分方程的解
是不是很逼近原微分方程的解在相应结点上的值
?为此,需要从两个方面,即差分方程的收敛性和稳定性来回答上述问题。
如果差分方程的解在步长Δx、Δt取得充分小时,和微分方程的解析解在某种意义上很接近的话,便说这种差分格式是收敛的。研究收敛性就是讨论当Δx→0、Δt→0时,差分方程的解和微分方程解的差(在一维条件下为)的绝对值在什么条件下趋近于零。其次,实际计算中由于只能用有限位计算,每一步都会有舍入误差,而且它还影响以后的计算结果。于是要考虑一个问题,当某一步结果本身有误差时,利用它去计算,若Δx和Δt固定,随着计算时间或计算次数的增加,误差是逐渐消除?还是逐步积累,愈变愈大?如是后者,则当t→∞时(或计算次数无限增多时),尽管某一步的误差很小,但其影响最终有可能达到十分可观的程度,使所得解面目全非。这时所考虑的差分格式便是不稳定的。显然,不收敛和不稳定的差分格式是没有实用价值的。
显式格式经证明,只有当0<λ≤1/2时才收敛和稳定。因此,Δt不能取的太大。
(二)隐式差分格式
如式(2-5)左端二阶导数取在tκ+1时间水平上[即用κ时段末t=(κ+1)Δt时刻的水头值],便得隐式差分方程为
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_29.jpg?sign=1739373091-lmQt1hVulTbEe5jSY0KSPIZ4unYkh6Hf-0-57c8160732b25f8d0e8b7caad6a44c8a)
截断误差为O(Δt)+O([Δx]2)。仍用式(2-11)定义的λ,则式(2-12)化为
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_30.jpg?sign=1739373091-tcUDw47G9thIMdV9bSomabJMD2ImJRej-0-6d5c852db58edc5012f801ef08ceebd5)
式(2-13)左端包含3个未知数,不能直接解出,所以称为隐式方法。必须对所有内结点(本例共有l-1个)都列出与式(2-13)相应的方程,并把边界条件
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_32.jpg?sign=1739373091-cvlyqS0XRJFNMn9Rl4KO9IrOoPAoWDBM-0-731b1d27a8845c05816a007fe1024111)
代入第一个和最后一个方程,形成由l-1个方程组成的方程组。联立求解,可得l-1个内结点上tκ+1时刻的水头值。所得代数方程组的系数矩阵只在三条对角线上有值,其余均为零,故称三对角线方阵。可用追赶法求解。
可以证明,隐式方法对任何λ都是收敛、稳定的,也就是说它的收敛、稳定是无条件的,Δt的取值不受Δx的严格限制。
(三)中心式差分格式
如果对在tκ和tκ+1的中点t=tκ+Δt/2处取中心差
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_34.jpg?sign=1739373091-V5IIouaIzFYgTmSNBcnwDU0qhULjnj16-0-a9cd77630d4e090125149c43efc3b6a7)
把沿t的正向在tκ+1/2处用Talyor级数展开,
沿t的负向在tκ+1/2处用Talyor级数展开,各取两项,两式相加,得
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_37.jpg?sign=1739373091-m59Z7eZHi9Lxh3KTCGGoSwDkXAsq42o9-0-cb610daf956ada3556cfab4333e6b418)
于是,式(2-5)可写成下列差分格式
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_38.jpg?sign=1739373091-IPe2xy2VZHWNz1W9cIj8BA9RszeRM8gL-0-6d4e685ab9d268301cc7d0f347cfcb0d)
截断误差为O([Δt]2)+O([Δx]2)称为Crank-Nicolson中心式差分格式或六点格式,形成的代数方程组的系数矩阵也是三对角线方阵。它和隐式方法一样也是无条件收敛、稳定的。
(四)加权显式-隐式格式
前面介绍的几种差分格式可以统一到下列一般形式中
![img](https://epubservercos.yuewen.com/424EE0/21277064901840706/epubprivate/OEBPS/Images/txt002_39.jpg?sign=1739373091-yzZjMF9HGZOVfemy56FMniSBrbum7ggG-0-d7e095eb876e59be905133f32f185669)
θ为权因子。上述格式称为加权显式-隐式格式。当θ=0时即为显式方法;θ=1时即为隐式方法;θ=1/2时即为中心式差分方法。不难证明,如取则式(2-15)是无条件稳定的,截断误差为O[(Δt)2+(Δx)4](Lapidus和Pinder,1982)。