变异Henry问题检验网格Peclet数在数值计算中的重要性

马 倩，李海龙

(1.生物地质与环境地质国家重点实验室/中国地质大学(北京)水资源与环境学院，北京 100083;2.地下水循环与环境演化教育部重点实验室/中国地质大学(北京)，北京 100083)

有限元法作为求解地下水流和溶质运移对流-弥散方程的常用数值方法，在水流模拟中有很多优点。它有固定的网格，通常满足质量守恒定律;可以精确、高效地处理以弥散为主的问题;便于编写程序并得以实现等［1］。但是，对于许多野外常见的以对流为主的问题，求解易引起显著的数值振荡。为此，许多学者进行了大量的计算实践和理论分析指出［2］:对于以对流占优的溶质运移问题，有限差分法和有限元法计算时皆不可避免地出现数值振荡。因此，掌握消除地下水流和溶质运移问题中浓度振荡的方法和理论依据是很有意义的。

Henry问题研究的是承压含水层中咸淡水静力平衡问题，在地下水流和溶质运移领域具有极高的代表性［3］。其常用的概念模型为:承压含水层均质各向同性，长2m高1m;上下顶板、底板为隔水边界;内陆边界为恒定的淡水流量边界;海岸边界为定水头定浓度的海水边界，海水水位1m［4］。它是一个理想的概念模型，为提高其实际适用性及科研价值(如提高模型参数敏感性)，很多学者在Henry问题的基础上，研究水文地质条件改变后的海水入侵问题，即变异Henry问题。Frind［5］研究的海水入侵问题，水动力弥散系数与原Henry问题不同;Huyakorn等［6］将右边界分为定浓度和零浓度通量边界条件两部分;Simpon和Clement［3］发现减小Henry问题中左边界的内陆淡水补给可以增加密度效应;Abarca等［7］综合考虑含水层渗透系数的各向异性及右边界水流方向对浓度的影响，发现变异后的Henry问题对密度变化更敏感。本研究对内陆和海岸边界条件改变后的Henry问题进行模拟求解，具体为:左边界为定水头边界、右边界水位为2m，发现该变异Henry问题对网格剖分要求更高，对网格Peclet数更为敏感(详见2.3.3)，为消除浓度振荡提供了实例，同时，验证和说明了网格Peclet数在数值计算中的重要性。

MARUN是专门用于地下水流和溶质运移模拟的数值程序，采用Galerkin二维有限元法求解对流-弥散方程，其可靠性已经得到各方面的验证［8～10］。Li等［11］利用该程序对均质海滩在海潮作用下的地下水盐动态进行了有效的定量模拟，Li等［12］、Guo 等［13］、Xia等［14］成功定量拟合了阿拉斯加不同海湾地下水位、盐分和示踪剂浓度的野外实测值，同时对海滩地下水动态水化学特征进行了定量刻画。

本研究以有限元法数值计算为例，采用MARUN模拟程序对变异Henry问题进行模拟求解，得到了不同网格剖分及水动力弥散系数特选节点处的浓度穿透曲线，发现了浓度振荡现象的存在。同时，分析找到了浓度振荡的原因及合理的消除方法。模拟结果及分析表明网格Peclet数能够有效地判定给定的网格剖分是否会引起浓度振荡，对有限元法数值计算的网格剖分具有指导意义。

1 模型介绍及网格剖分

1.1 模型介绍

变异Henry问题的概念模型具体为:承压含水层均质各向同性，长2m高1m;上下顶板、底板为隔水边界;内陆边界为定水头边界，淡水水头值为2.0632m;海岸边界为定水头定浓度的海水边界，海水水位为2m。原Henry问题的左边界流量为6.6×10-5m/s，本文提出的变异Henry问题所对应的左边界流量为1.2×10-4m/s，流速大约是原Henry问题的两倍。模拟区域及相关模型参数取值见图1。

图1 变异Henry问题的模拟区域及参数取值Fig．1 Model domain and parameters used in the modified Henry problem

Henry问题使用的水动力弥散系数为D1=1.8857×10-5m2/s，本研究考虑不同水动力弥散系数对模拟结果的影响，故加入水动力弥散系数D2=4.0×10-5m2/s。

1.2 网格剖分

本研究选用MARUN模拟程序采用Galerkin有限单元法计算模拟区域内各节点的浓度。剖分区域为长2m高1m的二维规则矩形，故采用直角三角形单元剖分，长度方向上的单元数是高度方向上的2倍，保证单元在整个区域上均匀分布;且所有的三角形单元为“等腰直角三角形”，不存在钝角三角形单元，是理想的三角形网格。研究共选取两种网格剖分:［Ⅰ］20行、40列;［Ⅱ］40行、80列。网格剖分Ⅰ共861个节点，1 600个单元，网格间距0.05m;网格剖分Ⅱ共3 321个节点，6 400个单元，网格间距0.025m。

2 结果与讨论

2.1 模拟求解

使用MARUN程序对以下三种情况的变异Henry问题进行模拟求解，得到各情况模拟区域浓度的时空变化情况。各节点的初始浓度为0mg/L，初始水头由左右边界水头线性插值得到。设置初始时间为0s，时间步长以1.1倍的因子增大或缩小，初始时间步长为10s，最大为 1 200s。

Case1:采用网格剖分Ⅰ，水动力弥散系数为D1=1.8857×10-5m2/s;

Case2:采用网格剖分Ⅱ，水动力弥散系数为D1=1.8857×10-5m2/s;

Case3:采用网格剖分Ⅰ，水动力弥散系数为D2=4.0 ×10-5m2/s。

Case1与Case2仅网格剖分不同，模型参数、初始条件及边界条件均相同。Case1与Case3仅水动力弥散系数不同，其他模型参数、网格剖分、初始条件及边界条件均相同。

2.2 结果

淡水从内陆边界流入，海水从海岸边界驱入含水层，与排泄的淡水混合。内陆和海岸边界条件始终不变，所以入侵海水与淡水流场最终可以达到动态平衡状态，即得到承压含水层中盐水楔形体与淡水流场间的稳态解。本研究模拟区域内节点初始浓度均为0mg/L，故各节点的浓度应先由0mg/L逐渐增加，达到稳定状态后不再随时间变化［7］。

有限元法在处理以对流为主的地下水流和溶质运移问题时，网格剖分往往引起数值计算的数值振荡，当出现很陡的浓度锋面时，也就是以对流为主时，振荡尤为明显［1］。为更清楚地研究和说明浓度振荡的普遍存在，选取流速较小的节点A(x=1.5m，z=0.05m)处的浓度穿透曲线进行研究。

2.2.1 不同网格剖分

如图2所示，Case1、Case2为不同网格剖分节点A浓度穿透曲线的对比，两种情况的节点A浓度由初始时刻的0mg/L近直线上升，至2h左右不再呈现明显上升趋势(图2a);Case1节点浓度呈现随时间在某值附近振荡，并未稳定于某一浓度值;Case2节点浓度最终稳定于0.4mg/L(图2b)。

如前所述，入侵海水与淡水流场达到动态平衡后，浓度不再随时间变化。而Case1节点浓度并未稳定，说明使用网格剖分Ⅰ有限元法数值计算出现了数值振荡，表现为节点处的浓度振荡;Case2节点浓度达到稳定，说明使用网格剖分Ⅱ有限元法数值计算未出现数值振荡，表现为节点处浓度最终稳定于某一浓度值，而不随时间振荡。Case1节点A浓度数值解明显偏离Case2浓度数值解，表明使用网格剖分Ⅰ有限元数值计算存在数值误差。即使用网格剖分Ⅰ有限元数值计算时既存在数值误差，又出现了浓度振荡。

图2 三种情况的节点A处浓度穿透曲线Fig．2 Breakthrough curve of point A at(a)initial time，and(b)whole time of three cases

2.2.2 不同水动力弥散系数

Case1、Case3为不同水动力弥散系数节点A浓度穿透曲线的对比，两种情况的节点A浓度由初始时刻的0mg/L近直线上升，至2h左右不再呈现明显上升趋势(图2a);Case1浓度呈现随时间在某值附近振荡，并未稳定于某一浓度值;Case3浓度最终稳定于0.38mg/L(图2b)。说明水动力弥散系数为 D1=1.8857×10-5m2/s时，有限元法数值计算出现了数值振荡，表现为节点处的浓度振荡;水动力弥散系数增大为D2=4.0×10-5m2/s时，有限元法数值计算得到了盐水楔形体与淡水流场的稳态解，数值解未出现数值振荡，表现为节点处浓度最终稳定于某一浓度值。另外，水动力弥散系数为D2时，若对网格加密后(用网格Ⅱ)进行模拟计算，则粗细网格剖分两种情况下的节点处浓度穿透曲线几乎完全重合，网格加密不再影响计算的精度。说明此时的网格剖分Ⅰ已满足模拟计算精度要求，得到了精确的稳态数值解。

2.3 讨论

2.3.1 浓度振荡消除方法

Case1、Case2仅网格剖分不同，数值计算得出的特选节点处浓度穿透曲线特征明显不同。Case1节点A浓度振荡，并未稳定于某一浓度值;Case2节点A浓度最终达到稳定。即使用网格剖分Ⅰ时数值计算出现浓度振荡，网格剖分Ⅱ在网格剖分Ⅰ的基础上将网格单元整体加密一倍，消除了浓度振荡。

Case1、Case3仅水动力弥散系数不同，特选节点处浓度穿透曲线特征明显不同。Case1节点浓度振荡，Case3节点浓度最终达到稳定。即水动力弥散系数为D1时数值计算出现浓度振荡，而水动力弥散系数D2在D1的基础上增加约一倍，消除了浓度振荡。

故有限元法求解地下水流和溶质运移问题时，若出现浓度数值解在某值附近振荡的现象，可以通过加密网格将其消除，同时须考虑水动力弥散系数设置是否合理，通过增加水动力弥散系数也可消除浓度振荡。另外，Case1与Case2采用不同网格处理同一问题，说明同一问题对网格精度的要求不同;Case1与Case3均采用网格剖分Ⅰ，仅模型参数不同，说明同一网格对不同模型参数的有效性不同。已有研究结果表明［3～7］，Henry问题在采用网格剖分Ⅰ时能得到精确的数值解，并不出现浓度振荡，说明即使是研究区域相同，不同的边界条件对网格精度的要求不同。网格剖分要具体问题具体分析，不能单纯认为边界条件相同，网格就适用于不同模型参数取值。

2.3.2 理论依据

网格Peclet数(Pe)是衡量对流项在溶质运移问题中主导程度的无量纲参数:

式中:v——达西速度;

D——水动力弥散系数;

Δx——网格间距。

对于纯对流问题，D→0，Pe→∞。随着物理弥散变得明显，Peclet数则变小。同时，Peclet数还依赖于网格间距Δx;Δx减小时，Peclet数也减小。即网格剖分越细，Peclet数也越小，这意味着采用较小的网格间距可以减小数值振荡。若Δx对应的Pe＜2，振荡就消除了［1］。

三种情况下对应的最大网格Peclet数如下:

Case1:采用网格剖分Ⅰ，Δx=0.05m;D=1.8857×10-5m2/s。速度v最大值为0.0015m/s，故模拟区域内最大网格Peclet数为Pe1=3.97。

Case2:采用网格剖分Ⅱ，Δx=0.025m;D=1.8857×10-5m2/s。速度v最大值为0.0015m/s，故模拟区域内最大网格Peclet数为Pe2=1.98。

Case3:采用网格剖分Ⅰ，Δx=0.05m;D=4.0×10-5m2/s。速度v最大值为0.0015m/s，故模拟区域内最大网格Peclet数为Pe3=1.88。

变异Henry问题的内陆和海岸边界条件不变，模拟区域内各点的浓度最终应达到动态平衡状态。而Case1浓度数值解出现振荡现象，与真实的物理意义不符，此时Pe1＞2;Case2浓度数值解最终稳定于某一浓度值，达到了动态平衡状态，符合Henry问题真实的物理意义，此时Pe2＜2;Case3浓度数值解最终稳定于某一浓度值，符合真实的物理意义，此时Pe3＜2。即不满足Pe＜2时，使用给定网格剖分进行数值计算时，数值解会出现浓度振荡(Case1);满足Pe＜2时，使用给定的网格剖分进行数值计算时，数值解不出现浓度振荡(Case2、Case3)，说明Peclet数能够有效地判定给定的网格剖分是否会引起浓度振荡。使用有限元法进行数值计算时，必须考虑Peclet数，使其尽量小于2。

2.3.3 变异Henry问题与Henry问题的对比

Henry问题与变异Henry问题的区别是，Henry问题的内陆边界为定流量边界，流量为6.6×10-5m/s，变异Henry问题为定水头边界，水头值为2.0632m;Henry问题的海水水位为1m，变异Henry问题的海水水位为2m;变异Henry问题左边界平均流速为1.2×10-4m/s，大约是原Henry问题左边界流速的两倍，故变异Henry问题的对流作用更强，为满足Peclet数小于2，要求网格剖分更密，即对网格剖分精度要求更高，更适于检验网格剖分是否会引起浓度振荡。海水入侵范围较Henry问题小，50%等氯线与底板交于1.55m处(Henry问题50%等氯线与底板交于1.40m处［3～7］)。变异Henry问题的求解为消除地下水流和溶质运移问题数值求解中的数值振荡提供了实例。

地下水数值模拟软件常用的求解方法是有限元法和有限差分法，二者在数值计算时均可能出现数值振荡。若出现数值解在某值附近振荡的现象，可以通过加密网格将其消除，同时需要考虑模型参数设置是否合理，增加水动力弥散系数也可消除。进行网格剖分时，即使研究区域相同，不同的边界条件、不同的水动力弥散系数对网格精度的要求不同，同一网格对不同模型参数的有效性也不同，要具体问题具体分析，不可一概而论。

3 结论

(1)浓度振荡可以通过加密网格、增加水动力弥散系数消除。

(2)即使研究区域相同，不同的边界条件、不同的水动力弥散系数对网格精度的要求不同，同一网格对不同模型参数的有效性也不同。

(3)网格Peclet数能够有效地判定给定的网格剖分是否会引起浓度振荡，对网格剖分具有指导意义。

(4)在地下水流和溶质运移有限元数值计算中，要综合考虑模型参数的合理性，分析地下水流场，对流作用较强的区域要加密剖分，满足Peclet数小于2，这样既可以消除浓度振荡，又保证计算速度。

［1］ 郑春苗，Gordon D Bennett.地下水污染物迁移模拟［M］.2版.北京:高等教育出版社，2009:109-118.［ZHENG C M，Gordon D Bennett.Applied Contaminant TransportModeling［M］.2nd ed.Beijing:Higher Education Press，2009:109-118.(in Chinese)］

［2］ 冯绍元.解地下水中溶质运移问题的特征有限元法［J］.武汉水利电力学院学报，1991，24(4):410-417.［Feng S Y.Method of Characteristic Finite Element for Solving Solute Transport in an Aquifer［J］.Journal of Wuhan institute of water conservancy and electric power，1991，24(4):410-417.(in Chinese)］

［3］ Simpson M J，T P Clement.Improving the worthiness of the Henry problem as a benchmark for densitydependent groundwaterflow models［J］. Water Resources Research，2004，40:W01504.

［4］ Henry H R.Effects of dispersion on salt encroachment in coastal aquifers［J］.US Geological Survey，1964，Water-Supply Paper 1613-C.

［5］ Frind E O.Simulation of long-term transient densitydependent transport in groundwater［J］.Advances in Water Resources，1982，5:73-88.

［6］ Huyakorn P，Andersen P F，Mercer J，et al.Saltwater intrusion in aquifers:development and testing of a three-dimensional finite element model［J］.Water Resources Research，1987，23(2):293-312.

［7］ Abarca E，Carrera J，Vila X S，et al.Anisotropic dispersvie Henry problem［J］.Advances in Water Resources，2006，30:913-926.

［8］ Boufadel M C，M T Suidan，A D Venosa.A numerical model for density-and-viscosity-dependent flows in two-dimensionalvariably-saturated porous media［J］.Journal of Contaminant Hydrology，1999a，36:1-20.

［9］ Boufadel M C，M T Suidan，A D Venosa.Numerical modeling of water flow below dry salt lakes:effect of capillarity and viscosity［J］.Journal of Hydrology，1999b，221:55-74.

［10］ Boufadel M C.A mechanistic study of nonlinear solute transport in a groundwater-surface water system under steady state and transient hydraulic conditions［J］.Water Resources Research，2000，36:2549-2565.

［11］ Li Hailong，Boufadel M C，Weaver J W.Tideinduced seawater-groundwater circulation in shallow beach aquifers［J］.Journal of Hydrology，2008，352:211-224.

［12］ Li Hailong，Boufadel M C.Long-term persistence of oil from the Exxon Valdez spill in two-layer beaches［J］.Nature Geoscience，2010，3(2):96-99.

［13］ Guo Qiana， Li Hailong， Boufadel M C， et al.Hydrodynamics in a gravel beach and its impact on the Exxon Valdez oil［J］. JournalofGeophysical Research，2010，115:C12077.

［14］ Xia Yiqiang，Li Hailong，Boufadel M C，et al.Hydrodynamic factors affecting the persistence of the Exxon Valdez oil in a shallow bedrock beach［J］.Water Resources Research，2010，46:W10528.