越流承压含水层不稳定流问题
有限差分法的稳定条件
唐仲华 向东进
(中国地质大学环境地质研究所，武汉 430074)
摘 要 证明了一维承压越流含水层系统不稳定流问题差分格式的稳定性条件.显式差分格式的稳定性条件是Δt｛4T/［μ(Δx)2］+b/μ｝≤2，隐式格式是无条件稳定的.以河间地块一维承压越流含水层为例，根据给出的稳定性条件进行了数值计算，结果表明，在满足稳定性条件时，其计算结果符合实际.
关键词 越流含水层，稳定性条件，有限差分方法.
中图法分类号 P614.2
第一作者简介 唐仲华，男，教授，1958年生，1989年毕业于中国地质大学研究生院，获硕士学位，主要从事应用数学及水文地质、工程地质方面的教学和研究工作.
0 引言
有限差分方法已广泛应用于求解地下水流动问题.对任何一种数值方法，都必须给出其稳定性条件，许多学者对此作过研究，如Remson等［1］对承压含水层一维有限差分法的稳定性进行了讨论.Rushton［2］研究了承压含水层二维有限差分格式的稳定性条件.Peaceman［3］讨论了对流-扩散型方程的有限差分格式的稳定性条件.最近，Chen等［4］研究了有限单元法求解反应输运方程的稳定性和收敛性.然而，对承压越流含水层系统地下水流动问题的有限差分方法的稳定性和收敛性条件的讨论还没有公开的报道.由Lax定理［5］知，收敛性条件与稳定性条件是一致的，即具有相容性.因此，我们只需讨论其中一个方面的问题.本文针对一维承压越流含水层不稳定流问题，证明了显式差分格式和隐式差分格式的稳定性条件.由地下水动力学理论容易推出，一维承压越流含水层不稳定流问题的数学模型为
其中：H是水头，T是水力传导系数，μ是储水系数，Kz，Mz分别是弱透水层渗透系数和厚度，Hz为相邻含水层的水头.
1 显式格式的稳定性条件
根据有限差分方法原理，容易推导方程(1)的显式差分格式为
我们知道，用差分方法求解不稳定流问题时，是从初始时刻的水头分布开始，以给定的时间步长Δt，逐个时段计算出以后各个时刻的水头值分布.所谓差分格式的稳定性是指，如果在求解差分方程时，某个时刻引进了误差ε，在以后的计算过程中，差分方程解的误差的传播不会增大，则称该差分格式是稳定的；否则是不稳定的［6］.不失一般性，我们设在初始时刻节点i处存在误差εi(i=1，…，N-1)，用 表示差分问题的近似解与差分问题精确解的差，则Vni满足
按稳定性定义，差分格式的稳定性也可表示为，对n＞0， 成立，其中 表示某种范数.
令 ，A为三对角矩阵，其主对角线上的元素为1-2γ-bΔt/μ，主对角线两边的元素为γ，则(3)式可写成
取2范数，由(4)式可得
由此可见，只要 A2 ≤1，则必有 Vn 2≤ V0 2
由计算数学理论［7］知道， A 2=A′A的特征值最大者的平方根.注意到A的对称性，由此即知 A 2为A的特征值中绝对值最大者的绝对值.因此，问题转化为求A的特征值.容易求得A的特征值为
因此，稳定性条件为
由此得到稳定条件为
2 隐式格式的稳定性条件
根据有限差分方法原理，容易推导方程(1)的隐式格式为
假设在时刻t=0时引入误差ε， 是差分方程的近似解与精确解之差，则Vni满足方程
定义B为三对角矩阵，其主对角线上的元素为1+2γ+bΔt/μ，主对角线两边的元素为-γ，则(10)式可写成
其中B-1是B的逆矩阵
由(11)式可以得到
在(12)式中两边取2范数，得
由(13)式可以看出，如果 B-1 2≤1，则 Vn 2≤ V0 2.
这说明，稳定性条件可以表示为 B-1 2≤1.
根据矩阵的特征值理论［7］知，如果λ是B的特征值，则λ-1是B-1的特征值，所以，问题归结为求B的特征值.
类似于求矩阵A特征值的方法，容易求得矩阵B的特征值为
于是，B-1的特征值为
所以，当
时，有 B-1 2≤1.
由于Δt≥0,b≥0,γ≥0，所以(16)式显然成立.这说明隐式格式是无条件稳定的，或者说其稳定性条件是Δt＞0.
图1 显式差分格式(a)和隐式差分格式(b)的计算结果
Fig.1The results caculated by explicit method (a) and implicit method (b) when time step is 0.45d and 0.50d
3 应用实例
某河间地块，有一均质、各向同性、等厚的承压含水层，底部有一相邻含水层，河水切割主含水层.取水头函数基准面为主含水层初始水位，则主含水层初始水头为H=0m，左端边界水头为1m，右端边界水头为0m，两河间距为1000m，主含水层导水系数为T=48m2／d，储水系数为μ=0.005,弱透水层系数为0.01m/d，弱透水层厚度为5m，相邻含水层水头为2m.
我们分别用显式和隐式格式求解了该问题.取Δx=100m，则显式格式的稳定性条件是Δt≤0.47d.图1a给出了取Δt=0.45d和Δt=0.50d时，显式格式的计算结果.图1b给出了取Δt=0.45d和Δt=0.50d时，隐式格式的计算结果.
由图1a可以看到，尽管时间步长只是略大于或小于稳定性条件的时间步长，但由两种时间步长得出的结果却完全不同.由图1b看到，由两种时间步长得出的结果基本相同.总之，如果所选参数符合稳定性条件，则其结果符合实际情况；否则，所得到的结果完全不符合实际.
4 结论
(1)一维越流承压含水层显示有限差分法的稳定性条件是Δt｛4T/［μ(Δx)2］+b/μ｝≤2.如果越流系数b=0，则稳定性条件简化为γ≤1/2，这就是无越流一维承压含水层显式差分法的稳定性条件.(2)隐式格式是无条件稳定，与显式格式相比，具有明显的优越性.
参考文献
1 Remson I, Hornberger C M, Molz F J. Numerical methods in subsurface hydrology. New York: Wiley_Interscience, 1971. 48～50
2 Rushton K R. Seepage and groundwater flow. New York: Wiley, 1997. 28～30
3 Peaceman D W. Fundamentals of numerical reservoir simulation. Amsterdam: Elsevier, 1997. 78～82
4 Chen Z X, Ewing R E. Stability and convergence of a finite element method for reactive transport in ground water. SIAM Journal on Numerical Analysis, 1997, 34(3): 881～904
5 Celia M A, Gray W G. Numerical methods for differential equations. New Jersey: Englewood Cliffs, 1992. 45～46
6 陈崇希，唐仲华.地下水流动问题数值方法.武汉：中国地质大学出版社，1990. 15～35
7 徐萃薇.计算方法引论.北京：高等教育出版社，1985. 142～143
STABILITY CONSTRAINTS OF FINITE DIFFERENCE METHOD
FOR SOLVING UNSTEADY CONFINED LEAKY AQUIFER SYSTEM
