越流承压含水层不稳定流问题有限差分法的稳定条件*唐仲华 向东进(中国地质大学环境地质研究所，武汉 430074) {SwQ[$k=_
摘 要 证明了一维承压越流含水层系统不稳定流问题差分格式的稳定性条件.显式差分格式的稳定性条件是Δ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.
