地下水溶质运移方程有限差分格式的实证研究

焦　甜,王军霞,唐仲华，刘　耿

(1.中国地质大学(武汉)环境学院，湖北 武汉 430074；2.中国地质大学(武汉)数学与物理学院，湖北 武汉 430074；3.北京得力合环境治理有限公司，北京 100025)



地下水溶质运移方程有限差分格式的实证研究

焦甜1,王军霞2,唐仲华1，刘耿3

(1.中国地质大学(武汉)环境学院，湖北 武汉 430074；2.中国地质大学(武汉)数学与物理学院，湖北 武汉 430074；3.北京得力合环境治理有限公司，北京 100025)

摘要：在求解对流占优的地下水溶质运移方程时，为克服有限差分法和有限单元法产生的过量和值弥散等虚假数值现象，出现了许多数值方法，其中包括Crank-Nicholson格式、上游加权法和人工扩散量法等。以一维地下水溶质运移方程为例，分析三种差分格式的截断误差、稳定性和收敛性，并通过数值试验，在不同的对流强度下，将不同算法的数值解与解析解进行对比，研究不同算法的优缺点和适用条件。结果表明：在对流占优较弱时，三种方法都适用，但Crank-Nicholson差分格式误差最小；当对流强度增大时，适用时间步长较小的人工扩散量法；当对流强度较大时，适用时间步长较大的人工扩散量法，或者空间步长较小的上游加权法。该研究可为解决实际问题时算法的选择提供依据。

关键词：地下水溶质运移方程；对流占优强度；有限差分格式；解析解；实证对比

随着经济的快速发展，地下水污染问题日益严重，由于地下水污染具有滞后性、隐蔽性以及一旦污染便难以清除的特点[1]，因此地下水污染问题得到了广泛关注。卢晓华[2]通过建立地下水溶质运移模型，对废水处理车间发生泄漏情况下地下水重金属污染情景进行了数值模拟与环境影响预测评价，对环保部门从源头上防治和管理潜在的地下水重金属污染、保护地下水环境具有重要的意义；李合莲等[3]采用有限元法对尾矿中238U在地下水中的迁移情况进行了模拟；邓鼎兴[4]采用有限元法对矿山地下水污染的分区治理效果进行了数值模拟检验，结果表明数值模拟预测结果与实际监测值较为接近，可为类似矿山地下水污染治理提供借鉴。但地下水溶质运移方程的求解精度会直接影响对水环境的分析，例如当对流占优时，传统的数值方法(如有限单元法和有限差分法)会产生数值弥散和数值波动等虚假数值现象。目前，针对对流作用占优的地下水溶质运移方程求解已经有许多方法，包括欧拉法、拉格朗日法和混合-欧拉拉格朗日法[5]等，但这些方法仍然存在许多科学问题亟待解决。

有限差分法和有限单元法[6]等是常用的欧拉法。对于有限差分法，有学者提出了上游加权的方法,此后不少学者也对此进行了研究[7]，结果表明此方法消除了数值波动，但却增大了数值弥散。相比之下有限单元法精度高、较灵活，但当对流占优时，时间导数若采用中心差分会出现过量，向后差分则会出现数值弥散，对此很多学者提出了不同的上游加权方法[8]，但由于上游加权法结点所在的边线不一定是流线方向，就出现了“侧风现象”，为了消除这种现象，Brooks等[9]提出了沿流线上游加权的彼得罗夫-迦辽金方法，但实践表明这种方法也是在消除数值振荡的同时导致了更大的数值弥散。梅一等[10]在上游加权法的基础上引入一个数值弥散因子，并适当调节因子的大小来有效控制数值振荡，减少数值弥散,并且达到了满意的精度。动坐标法通过控制坐标系的运动速度有效减小了求解对流项引起的数值误差;网格变形法通过固定坐标系，使网格结点自动跟着锋面的移动而聚在锋面附近，来减小锋面附近的Peclet值，消除过量和数值弥散问题。虽然这两种方法处理对流占优问题时比较有效，并且保持了较高的精度，但实际问题中计算工作量比较大，难以推广到三维问题上。由于对流项占主导地位时，地下水溶质运移方程接近于一个双曲线方程，所以启发人们利用双曲线的特征线法来求解，Garder首先用此法求解可溶混溶质运移迁移，此种方法原理简单、程序冗长，但能有效减少数值弥散，也能处理各种复杂条件。Douglas等[11]提出了特征有限元法和有限差分法，其数值计算方面的有效性曾引起了许多学者的关注，并得到了广泛的应用研究[12]。随机步行法将溶质运移视为大量示踪剂特征点运动的平均结果，当特征点足够多时，基本上能刻画溶质运移过程，优点是由于特征点在研究区域内连续运动，所以计算浓度时只需要统计区域上特征点的数目，不用在每个时间步长内计算整个研究区的浓度。Kinzellbach将随机步行法与有限差分法、有限单元法和特征值法进行了对比，开发出与MODFLOW配合使用的MODWALK软件，但是其计算精度取决于示踪点数目，如果数目不够，统计出的浓度解会比较粗糙。Leismann等[13]提出了引入人工扩散量和加权法来求解，定义并引入人工扩散项进行校正。刘扬[14]提出一种加罚函数的Crank-Nicholson差分格式，形成的代数矩阵对角元严格占优，在保持数值解稳定性的基础上有效克服了数值振荡。

对于各种方法的优劣性，成建梅等[15]对溶质运移问题数值解从欧拉法、拉格朗日法、混合-欧拉拉格朗日法三大方面进行了现状综述；刘苑等[16]对比了MT3D模块中MOC、MMOC、HMOC和UFDM 4种数值方法；宁一鸣[17]对比了二维问题的弥散占优和对流占优的几种数值解差分格式，但没有讨论不同对流强度的情况，并探讨了一种基于三角形网格的上游加权有限差分格式。除此之外，许多学者也从各个方面研究了对流占优数值解差分格式的特点，但是对不同对流占优强度下不同数值解差分格式的适用性没有做详细的讨论。因此，本文在前人研究工作的基础上，对比了三种差分格式的截断误差、稳定性和收敛性，并通过数值试验在不同的对流强度下，将不同算法的数值解与解析解进行对比，研究了不同算法的优缺点和适用条件，讨论了时间步长对解的影响，以便于研究者在实际问题中选择适合的算法以及网格剖分大小，从而提高模拟的精度。

1三种差分格式

对于一维对流-弥散问题[18]，地下水溶质运移方程为

(1)

式中：C为浓度;DL为水动力弥散系数;u为地下水流速;x为空间距离;t为时间。

对于公式(1)分别采用Crank-Nicholson格式(C-N格式)、上游加权法、人工扩散量法，可得到以下三个差分方程：

(1) Crank-Nicholson格式的差分方程为

-ACi-1,n+1+BCi,n+1-DCi+1,n+1=ACi-1,n+

(4-B)Ci,n+DCi+1,n

(2)

式中：i为空间距离x的剖分离散结点;n为时间t的剖分离散结点;A、B、D为方程的系数，具体取值如下：

(2) 上游加权法的差分方程为

(3)

由上式可以发现，当其他参数保持不变，u越大时ω取值要求越高。

(3) 人工扩散量法的差分方程为

-ACi-1,n+1+BCi,n+1-DCi+1,n+1=ACi-1,n+

(4-B)Ci,n+DCi+1,n

(4)

式中：A、B、D为方程的系数，具体取值如下：

所谓过量就是由于数值解的振动，出现了超过最大浓度和小于最小浓度的计算值，违背了基本的物理意义。数值弥散则是由截断误差引起的，这里我们采用Pe=Δx/α来描述对流和弥散的相对大小，Δx为空间步长，α为纵向弥散率。Peclet数能够有效地判定给定的网格剖分是否会引起浓度振荡，对网格剖分具有指导意义[19]。按照Peclet数可以将对流-弥散方程分成两大类：当Peclet数小于或等于2时，为弥散占优方程；当Peclet大于2时，为对流占优方程。

2三种差分格式的误差及稳定性分析

2.1Crank-Nicholson格式

采用泰勒级数展开法分析可知，由于Crank-Nicholson格式空间项采用的是中心差商，所以其截断误差为O(Δt2)+O(Δx2)，下面来分析其稳定性。

首先扩充函数的定义域，使其在整个实轴上有定义，令

U(x,tn+1)=Ci,n+1

所以公式(2)可以改写为-AU(x-Δx,tn+1)+BU(x,tn+1)-DU(x+Δx,tn+1)=

AU(x-Δx,tn)+(4-B)U(x,tn)+DU(x+Δx,tn)

令

1-cos(kΔx)≥0,因为上式的分母为正，所以有|G|2-1≤0。

得出Crank-Nicholson格式是无条件稳定的。

Crank-Nicholson格式是隐式的，即

F2(Ci+1,n+1,Ci,n+1,Ci-1,n+1)=F1(Ci+1,n,Ci,n,Ci-1,n)

Ci+1,n,Ci,n,Ci-1,n和两端函数已知，只有Ci+1,n+1,Ci,n+1,Ci-1,n+1待求。因为高斯-塞德尔迭代法[20]的收敛速度比较快，所以本文采用高斯-塞德尔迭代法来求解方程组。在此种迭代方法下，为了满足迭代的收敛性，系数矩阵必须是严格对角占优的，即|B|>|A|+|D|。将系数代入,并化简得

(5)

因此求解方程组时采用高斯-塞德尔法迭代，要求时间和空间剖分Δt、Δx以及弥散系数DL和实际速度u满足式(5)Crank-Nicholson格式收敛。这里讨论的计算时迭代的收敛性分析和差分格式的数值解对于精确解的收敛性是不同的概念。由于在建立差分方程以后，给定初值和边界条件后，就必须通过解方程组来计算每个点的数值，本文在解方程时用高斯-塞德尔迭代法，所以在计算时就有一个迭代过程的收敛性问题，在文中“计算时迭代的收敛性”是指方程组的数值解是否收敛到方程组的精确解。另外两种方法的收敛性也表示的是同样的含义。

2.2上游加权法

由于公式(3)取的不是中心差分，所以截断误差是一阶的，即为O(Δt)+O(Δx)。上游加权法的稳定性分析与Crank-Nicholson格式类似，首先来扩充函数的定义域，使其在整个实轴上有定义，每一项用Fourier积分来表示，引入欧拉公式，得出过度因子为G=m(4)/{-m(1)cos(kΔx)+m(2)-m(3)cos(kΔx)+i[m(1)sin(kΔx)-m(3)sin(kΔx)]}其中:

将系数代入上式，并整理得

H2[1+cos2(kΔx)]-2H2cos(kΔx)=

H2[1+cos2(kΔx)-2cos(kΔx)]≥0;

当H≥0时，过度因子等式右端满足小于或等于1，故增长矩阵一致有界，式(3)是稳定的，因此上游加权法的稳定性条件是满足下式：

(6)

由于差分格式是隐式的，与Crank-Nicholson格式类似，所以方程组采用高斯-塞德尔法迭代时，为了满足迭代的收敛性，系数矩阵必须是严格对角占优的,即

当认为u>0时，化简为

(7)

因此求解方程组时采用高斯-塞德尔法迭代，要求时间和空间剖分Δt、Δx以及弥散系数DL和实际速度u满足式(7)时结果收敛。通过前面分析可知上游加权法的截断误差为O(Δt)+O(Δx)，说明此方法与Crank-Nicholson格式相比解的精度降低了，尤其在速度很大时由(7)式收敛条件可知上游加权法就达不到收敛的要求。

2.3人工扩散量法

(8)

用泰勒公式展开上式的C(x+Δx)与C(x-Δx)，并整理得

说明方程式(8)右端两项都很准确，再分析左端的时间项，对时间项进行泰勒展开并代入方程式(8)左端，有

所以式(1)真正离散后的方程变为

可见，为了校正这种误差，引入人工扩散系数，就此推出了人工扩散量方法的表达式。

(9)

同理，若求解方程组时采用高斯-赛德尔法迭代，要求时间和空间剖分Δt、Δx以及弥散系数DL和〗实际速度u满足公式(9)人工扩散量法收敛。

3数值试验与结果分析

本文考虑下述地下水溶质运移问题，即半无限长均质砂柱，地下水以实际平均流速u稳定流动，初始浓度为0，其一端为定浓度边界，假设无穷远处浓度分布为零，求其浓度分布。该问题的解析解为

(10)

为了用无限长的解析解近似有限长的实际问题，缩小误差，本文通过多次试验确保数值解与解析解在有限长的边界处附近浓度都为0，取时间为4 d，空间为60 m，并取参数C0=1 mg/L、u=6 m/d、Δt=0.5 d、Δt=0.1 d、Δx=2 m，这里取DL=uα[21]，同时按照Pe=Δx/α取不同的弥散率来代表不同的对流强度，α取值分别为2 m、0.5 m、0.125 m和0.062 5 m。通过计算验证表明其都满足三种差分格式的稳定性和收敛性条件，对应的对流占优强度分别为1、4、16和32。图1和图2为对流强度分别为1和4时数值解与解析解的浓度对比图，其中图1和图2中(a)和(b)分别为取时间段为1 d和3 d时三种差分格式随空间位置的浓度分布图。

图1　当对流强度为1时数值解与解析解的浓度对比图Fig.1　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 1

通过对比分析图1和图2可见，当对流占优较弱时，三种方法都很适用，其数值解与解析解几乎重合，差分格式误差都很小；当对流强度为4时C-N格式开始出现过量现象，上游加权法差分格式误差最大。通过将图2(b)与图2(c)进行对比分析可见，缩小时间步长，三种方法的误差都相对减小。

图3和图4为对流强度分别为16和32时数值解与解析解的浓度对比图，其中图3和图4中(a)和(b)为取时间为3 d、时间步长分别为0.5和0.1时三种差分格式随空间位置的浓度分布图。

通过对比分析图3(a)和图3(b)可见，人工扩散量法对时间步长比较敏感，时间步长缩小时，差分格式误差显著变小；相比之下上游加权法在缩小时间步长时，差分格式误差变化并不明显。

图3　当对流强度为16时数值解与解析解的浓度对比图Fig.3　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 16

图4　当对流强度为32时数值解与解析解的浓度对比图Fig.4　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 32

通过对比分析图3和图1及图2可见，随着对流强度的增大，C-N格式的过量现象越发明显，数值解与解析解的总体误差有所增大，这时人工扩散量法数值解与解析解的误差相比较小，尤其是在时间步长较小时，此种方法开始发挥优势，上游加权法精度不如人工扩散量法。

通过对比分析图4(a)和图4(b)可见，当对流强度为32时，人工扩散量法在时间步长较大时没有出现过量现象，但在时间步长过小时也出现了过量现象，说明了在减小数值弥散的同时增大了解的振荡；只有上游加权法没有出现过量现象，但是误差稍大，因为上游加权法对空间步长比较敏感，这时可以缩小空间步长，通过更小的网格剖分，来减弱对流强度，以减小差分格式误差。

4结论

采用数值法求解对流占优的地下水溶质运移问题时，难点是克服求解过程中的虚假数值现象。本文首先分析了差分格式的截断误差、稳定性和收敛性，再进行数值试验，结果表明:①在对流占优比较弱时，三种方法都适用，但C-N格式最接近理论值；②当对流强度增大时，例如数值试验中对流强度大于4时，建议使用时间步长较小的人工扩散量法，因为C-N格式开始出现过量现象，其误差开始增大，而上游加权法差分格式的误差比前两种方法稍大；③当对流强度较大时，例如数值试验中对流强度大于32时，建议使用时间步长较大的人工扩散量法，或者使用空间步长较小的上游加权法，通过调节空间步长来减弱对流强度，再进行求解。

对于一个实际问题，弥散系数和流速是一定的，只有网格剖分可以调节，从前面分析可知当弥散系数和流速已知时，可以适当调节网格剖分，在满足稳定性和收敛性条件下，根据实际问题在保证适当计算量的前提下提高算法的精度。

参考文献：

[1] 王焰新.地下水污染与防治[M].北京：高等教育出版社,2007.

[2] 卢晓华.基于数值模拟的企业地下水重金属污染的环境影响预测评价[J].安全与环境工程,2014,21(1):93-97.

[3] 李合莲,陈家军.铀尾矿库中238U运移数值模拟[J].安全与环境工程,2007,14(1):36-39.

[4] 邓鼎兴.某矿山地下水污染分区治理效果的数值模拟[J].安全与环境工程,2013,20(4):27-31.

[5] 魏恒,肖洪浪.地下水溶质运移模拟研究进展[J].冰川冻土,2013,35(6):1584-1586.

[6] Hossain M A,Barber M E.Optimized Petrov-Galerkin model for advective-dispersive transport[J].AppliedMathematicsandComputation,2000,115(1):1-10.

[7] Dehghan M.Weighted finite difference techniques for the one-dimensional advection-diffusion equation[J].AppliedMathematicsandComputation,2004,147(2):307-319.

[8] Sun N Z,Yeh W G.A proposed upstream weight numerical method for simulation pollutant transporting ground water[J].WaterResourcesResearch,1983,19(6):1489-1500.

[9] Brooks A N,Hughes T J R.Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J].ComputerMethodsinAppliedMechanicsandEngineering,1982,32(1/2/3):199-259.

[10]梅一,吴吉春.地下水溶质运移数值模拟中减少误差的新方法[J].水科学进展,2009,20(5):640-645.

[11]Douglas J,Russell T F.Numerical method for convection-dominated diffusion problems based on combing the method of characteristics with finite element or finite difference procedures[J].SIAMJournalonNumericalAnalysis,1982,19(5):871-885.

[12]Tsai T L,Yang J C,Huang L H.Characteristics method using cubic-spline interpolation for advection-diffusion equation[J].JournalofHydraulicEngineering,2004,130(6):580-585.

[13]Leismann H M,Frind E O.A symmetric-matrix time integration scheme for the efficient solution of advection-dispersion problems[J].WaterResourcesResearch,1989,25(6):1133-1139.

[14]刘扬.对流-扩散方程的新型Crank-Nicholson差分格式[J].数学杂志,2005,25(4):464-467.

[15]成建梅,胡进武.饱和水流溶质运移问题数值解法综述[J].水文地质工程地质,2003,30(2):99-101.

[16]刘苑,武晓峰.地下水中污染物运移过程数值模拟算法的比较[J].环境工程学报,2008,2(2):230-234.

[17]宁一鸣.地下水溶质运移对流-弥散方程数值解方法对比研究[D].武汉:中国地质大学,2008.

[18]陈崇希,李国敏.地下水溶质运移理论及模型[M].武汉:中国地质大学出版社,1996.

[19]马倩,李海龙.变异Henry问题检验网格Peclet数在数值计算中的重要性[J].水文地质工程地质,2014,41(3):1-5.

[20]李庆扬,王能超,易大义.数值分析[M].北京:清华大学出版社,2008.

[21]Fetter C W.ContaminantHydrology[M].New Jersey:Prentice-Hal,1999:45-186.

Empirical Study of the Finite Difference Schemes of Solute Transport Equation in Groundwater

JIAO Tian1,WANG Junxia2,TANG Zhonghua1,LIU Geng3

(1.SchoolofEnvironmentalStudies,ChinaUniversityofGeosciences,Wuhan430074,China;2.SchoolofMathematicsandPhysics,ChinaUniversityofGeosciences,Wuhan430074,China;3.BeijingDeliheEnvironmentalGovernanceLtd.,Beijing100025,China)

Abstract:In solving the advection-dominated solute transport equation in groundwater,in order to overcome the excess and numerical dispersion produced by finite-difference and finite element methods,many numerical methods are adopted,such as Crank-Nicholson scheme,upstream weighting method and artificial diffusion method.Taking the one-dimensional groundwater solute transport equation as an example,this paper analyses the truncation error,stability and convergence of the three difference schemes,and compares the numerical solutions of different algorithms with the analytical solutions in different advection-dominated intensities through numerical experiments in order to explore the advantages and disadvantages and applicable conditions of different algorithms.The results indicate that when the advection-dominated intensity is weak,all the methods are applicable,while the error of CN scheme is the minimum;with the increasing advection-dominated intensity,the artificial diffusion method with smaller time step is applicable;once the advection-dominated intensity is relatively strong,the artificial diffusion method with larger time step or the upstream weighting method with smaller space step is applicable.This research provides a basis for the selection of algorithm for solving the practical problems.

Key words:groundwater solute transport equation;advection-dominated intensity;finite difference scheme;analytical solution;empirical comparison

文章编号：1671-1556(2016)03-0017-07

收稿日期：2015-10-11修回日期：2016-03-20

基金项目：中国博士后科学基金项目(201404406064);中国地质调查局项目(1212011121142)

作者简介：焦甜(1993—)，女，硕士研究生，主要研究方向为地下水及溶质运移数值模拟。E-mail:jiaotianedu@sina.cn

中图分类号：X143;P641.2

文献标识码：A

DOI：10.13578/j.cnki.issn.1671-1556.2016.03.003

通讯作者：唐仲华(1958—)，男，教授，主要从事渗流理论及其数值模拟方面的研究。E-mail:zhhtang@cug.edu.cn