CN104316658B - 一种模拟地下水一维溶质运移过程的方法 - Google Patents

Abstract

本发明公开了模拟地下水一维溶质运移过程的方法，属于地下水动力学领域。该方法利用马尔科夫链的基本性质，通过利用粒子的跳动实现对溶质运移的刻画，根据大量的随机游走的粒子的概率分布对溶质的运移过程进行模拟，完全不同于传统的数值解法，利用本发明所得到的结果和传统的对流弥散方程的解析解有很好的匹配，而本发明所带来的计算资源的消耗，相比于传统的一些方法，则大大降低。

Description

一种模拟地下水一维溶质运移过程的方法

技术领域

本发明涉及地下水动力学领域，涉及溶质运移的一维对流弥散方程的模拟方法。

背景技术

中国华北和西北的总体供水量有50％以上来自地下水。由于地下水环境所表现出来的隐蔽性和系统复杂性，致使地下水的污染修复极其困难和昂贵。地下水污染所引起的生态环境破坏和人体健康的危害是社会经济可持续发展的巨大挑战。初步调查显示，我国有90％的城市地下水遭到不同程度的有机和无机有毒有害物质的污染，并且已呈现出由点状向面状扩展的趋势。因此，研究地下水中溶质运移的规律，揭示地下水中污染物的迁移机理，是具有非常重大的实际意义的。

地下水溶质运移最基本的控制方程为对流弥散方程，其一维形式为：

$\mathrm{\θ}\frac{\∂C}{\∂t}=-v\frac{\∂C}{\∂x}+D\frac{{\∂}^{2}C}{\∂x^2}---\left(1\right)$

如果将式子两边同时除以θ，可以得到：

$\frac{\∂C}{\∂t}=-v^{*}\frac{\∂C}{\∂x}+D^{*}\frac{{\∂}^{2}C}{\∂x^2}---\left(2\right)$

其中 $v^{*}=\frac{v}{\mathrm{\θ}},D^{*}\frac{D}{\mathrm{\θ}}.$

在表达式(1)中，C表示溶质浓度，D表示弥散系数，v表示对流速度。对流弥散方程广泛地被用于溶质运移问题。由于解析解一般仅限于比较简单和理想的情况下，因此数值解成为解决实际问题的最基本的和最常用的方法和手段。各种解对流弥散方程的方法可以归结为三类：欧拉型，拉格朗日型以及二者的混合。欧拉型是以在空间中的固定坐标系作为参考系，主要包括有限差分法、有限元法、边界单元法和有限分析法；拉格朗日型则是以跟随着流体运动的运动坐标系作为参考系，在运动的坐标系中观察对流—弥散的现象，通过控制坐标系的运动速度，可以使对流项不出现或者使其变得很小。对流—弥散方程中同时包含了双曲项以及抛物项，这个方程的数值解已经被很多科研工作者研究过。基于欧拉坐标系，这个方程的数值解最常用的方法是有限元和有限差分，而这两种方法在实际应用中碰到的最常见的问题莫过于数值弥散，同时当所给定的计算条件不满足收敛性条件时，甚至会出现数值震荡，从而偏离我们所需要的正确解[薛禹群，谢春红水文地质学的数值方法北京：煤炭工业出版社，1980]。此外，当弥散项占主导时，各种数值解法均能得到满意的结果；反之，若对流项占主导时，各种常见的数值解法解的精度与效果不同，有些会产生数值弥散过量及解不收敛的问题。拉格朗日方法能够有效地消除数值弥散，却不能保持质量守恒，且难以处理条件比较复杂的、实际的问题[S.P.Neumann,Adaptive Eulerian Lagrangian finite-element method foradvection dispersion,Int.J.Numer.Meth.Eng.20321(1984)]。针对数值解中出现的种种问题，人们做了很多研究，发展出隐格式优先差分算法、特征方程法等算法，推动了算法领域的革新；同时也有利用随机游走理论、格子玻尔兹曼方法对该方程在微观尺度求解的理论。一系列的算法和理论发展，已经基本上克服了收敛性的问题。但是所有这些方法和理论，都没有提供一个直观的数学模型，其物理意义在水文地质学领域也没有得到充分体现。另外，即使解决了收敛性的问题，但是数值弥散的存在也是这些理论不得不面对的另一个问题。隐式差分格式虽然可以有效解决数值解收敛性的问题，但是在解决收敛性问题的同时却带来了显著的数值弥散[孙讷正地下水污染——数学模型和数值方法北京：地质出版社1989；N.Sun.andW.Yeh,A proposed upstream-weight numerical method for simulating pollutant transport ingroundwater,Water Resour.Res.191489(1983)]。特征方程的方法则是分别利用欧拉法和拉格朗日法来处理弥散项和对流项，然而这种方法最大的问题在于当初条件和边条件不满足特定要求时，将会出现质量不守恒的情况，而这在污染物溶质运移过程中是不可能发生的。

发明内容

本发明的目的在于提出一种模拟地下水一维溶质运移过程的方法，该模拟方法是基于马尔科夫过程的随机游走方法，得到的数值解可以和经典的对流弥散方程的解析解完全匹配，和经典的有限差分方法相比，在精度相同的情况下，可以大大节省计算资源，提升计算效率。

为达到本发明的上述目的，本发明提出了的模拟地下水中一维溶质运移过程的方法，具体步骤：

(1)在初始时刻释放溶质，此释放过程可以是瞬时的，也可以是持续的，既可以是点源，又可以是线源；

(2)设定时间步长Δt及空间步长Δx，且时间步长和空间步长仅需满足稳定性条件即可，其中D^*是对流弥散方程中的弥散系数；

(3)在离散化过程中Δx为空间维度上的单位1，且取Δt为时间维度上的单位1，设在一维尺度范围内，粒子向左跳动的概率是μ，向右跳动的概率是λ，根据公式计算，v^*和D^*已知的情况下，λ和μ则可以相应求得，从而可以得到在任意时刻的溶质运移的对流弥散方程的表达式：

$\frac{\∂P(x,t)}{\∂t}=P(x-1,t)F+P(x,t)R+P(x+1,t)B---\left(1\right)$

其中，F＝λ,B＝μ,R＝-λ-μ。这样，只需要知道前一时刻一维尺度上任一点的浓度值，即可得到后一时刻一维尺度上任意点的浓度值；而在时刻0时，以为尺度上任一点的浓度值是已知的。

(4)将0时刻的数据带入式子(1)，进行迭代，即可得到后面任意时刻的任意点的浓度值，从而实现对溶质运移的过程进行模拟。

本发明提出了一种全新的对一维溶质运移过程进行模拟的方法，在求解过程中集中体现了马尔科夫链的基本性质，即若当前状态已知，未来的状态独立于过去的状态。这种方法不同于以往的有限差分法和格子玻尔兹曼算法等，是对随机游走方法的一种扩展和延伸。

本发明通过利用粒子的跳动实现对溶质运移的刻画，使用粒子的概率分布密度对溶质浓度的分布进行预测，同时这种方法所得到的结果和传统的对流弥散方程的解析解有很好的匹配，而采用这种方法所带来的计算资源的消耗，相比于传统的一些方法，则大大降低。

传统的数值解法，在解决对流弥散方程的过程中难免会出现严重的数值弥散，而本发明则完全克服了这样的弊端，因为本发明所采取的办法是使用粒子的随机游走，根据大量的随机游走的粒子的概率分布对溶质的运移过程进行模拟，完全不同于传统的数值解法的求解过程，故而完美地解决了这个问题。

附图说明

图1是粒子运动过程中从x_i运动到相邻点的运动示意图；

图2是本发明利用马尔科夫随机过程得到的穿透曲线和采用有限差分法及解析解得到的穿透曲线的对比效果图；

图3是实施例公开的在某些时刻的一维浓度分布状况图。

具体实施方式

下面通过实例对本发明做进一步说明。需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。

本发明根据一维对流弥散方程

$\mathrm{\θ}\frac{\∂C}{\∂t}=-v\frac{\∂C}{\∂x}+D\frac{{\∂}^{2}C}{\∂x^2}---\left(2\right)$

如果将式子两边同时除以θ，可以得到：

$\frac{\∂C}{\∂t}=-v^{*}\frac{\∂C}{\∂x}+D^{*}\frac{{\∂}^{2}C}{\∂x^2}---\left(3\right)$

其中 $v^{*}=\frac{v}{\mathrm{\θ}},D^{*}\frac{D}{\mathrm{\θ}}.$

又由科尔莫格罗夫方程可得：

$\frac{\∂P(x,t)}{\∂t}=\mathrm{\λP}(x-1,t)+\mathrm{\μP}(x+1,t)-(\mathrm{\λ}+\mathrm{\μ})P(x,t)---\left(4\right)$

也即：

$\frac{\∂P(x,t)}{\∂t}=-(\mathrm{\λ}-\mathrm{\μ})\frac{P(x+1,t)-P(x-1,t)}{2}+\frac{\mathrm{\λ}+\mathrm{\μ}}{2}(P(x+1,t)-2P(x,t)+P(x-1,t))---\left(5\right)$

又因为：

$\begin{array}{c}\frac{\∂P(x,t)}{\∂x}\≈\frac{\frac{1}{2}(P(x+\mathrm{\Δx},t)+P(x,t))-\frac{1}{2}(P(x-\mathrm{\Δx},t)+P(x,t))}{\mathrm{\Δx}}\\ =\frac{P(x+1,t)-P(x-1,t)}{2}\\ \frac{{\∂}^{2}P(x,t)}{\∂x^2}\≈\frac{P(x+\mathrm{\Δx},t)-2P(x,t)+P(x-\mathrm{\Δx},t)}{{\left(\mathrm{\Δx}\right)}^2}\\ =P(x+1,t)-2P(x,t)+P(x-1,t)\end{array}$

则

$\frac{D^{*}}{\mathrm{\Δ}{x}^2}=\frac{\mathrm{\λ}+\mathrm{\μ}}{2},\frac{v^{*}}{\mathrm{\Δx}}=\mathrm{\λ}-\mathrm{\μ}$

从而推导出公式(1)。

本发明中粒子的跳动如图1所示，无关于溶质运移过程中的弥散系数和对流速度，仅仅和粒子向各个方向的跳动概率有关，且此概率必须为小于1的正值。

以下结合具体实施例进一步详细描述本发明：

(1)在初始时刻释放溶质质量为1的点源溶质；

(2)设定时间步长Δt＝0.25秒及空间步长Δx＝1米，但是需要说明的是这里的时间步长和空间步长可以不给出具体的单位，而仅仅采用单位长度即可(例如，可以设定0.25秒为一单位的时间长度，设定1米为以单位的空间长度)，设定弥散系数D^*＝0.15米²/秒(注释：弥散系数在实际的应用中并非是固定的数值，基本上要根据实际情况来确定，但是在模拟过程中，一般指定了弥散系数的话，就会用这个指定的弥散系数，所以，在模拟过程中，这是一个确定的已知的数值)、对流速度v^*＝0.2米/秒(流速也是根据实际情况确定的，一般情况下的合理范围是0.001米/秒到1米/秒，本例采用的数值是给定的一个数值而已)，则在这样的设定下，满足稳定性条件。

(3)在运移过程中，据本发明所提出的计算关系(本关系的推导见于具体实施方式)，则可以设粒子向左跳动的概率是μ＝0.05，向右跳动的概率是λ＝0.25(μ和λ是根据D^*和v^*计算得到的，但是在一个具体的模拟过程中，每次给定了D^*和v^*，则μ和λ就可以相应确定，在该时间点和空间点的话μ和λ也就可以根据D^*和v^*来唯一确定)，从而可以得到在任意时刻的溶质运移的对流弥散方程的表达式：

$\frac{\∂P(x,t)}{\∂t}=P(x-1,t)F+P(x,t)R+P(x+1,t)B---\left(1\right)$

其中，F＝λ,B＝μ,R＝-λ-μ。这样，只需要知道前一时刻一维尺度上任一点的浓度值，即可得到后一时刻一维尺度上任意点的浓度值；而在时刻0时，以为尺度上任一点的浓度值是已知的。

(4)将0时刻的数据带入式子(1)，进行迭代，即可得到后面任意时刻的任意点的浓度值，从而对对流弥散方程进行求解，进而对溶质运移的过程进行模拟。

当然，本模拟过程不仅仅可以针对瞬时点源，也可以针对持续性的点源，即在源处维持一个恒定的溶质的浓度，并且在后续的溶质运移过程中点源的的浓度不随时间发生改变。此时，若再进行迭代的话，那么在x＝0处将是一个定浓度的边界条件和初始条件，并不随着迭代过程发生任何变化。

首先，以一维瞬时点源的马尔科夫过程模型得到的结果，和解析解以及有限差分得到的数值解进行对比。设在t＝0时，在位于x＝0的点处释放溶质总量为单位1的溶质，同时取距离源点10个单位长度的一个点作为观测点，观测这个点300个单位时间内的溶质浓度的变化(即画出“穿透曲线”)，采用三种方法分别计算这个点的穿透曲线，进行对照，也就是图2所展示的内容。

图3所显示的内容是在某一个时间点时，在所研究的区域范围内，任何一点的溶质浓度的对比。这个模型的初始条件和边界条件都与图2中所示的模型一模一样，图3显示了本发明利用马尔科夫随机游走方法所得的结果和解析解以及优先差分数值解所得的结果相一致。

虽然本发明已以较佳实施例披露如上，然而并非用以限定本发明。任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案作出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (3)

1.一种模拟地下水一维溶质运移过程的方法，具体包括以下步骤：

(1)在初始时刻释放溶质；

(2)设定时间步长Δt及空间步长Δx，且时间步长和空间步长满足稳定性公式其中D^*是对流弥散方程中的弥散系数；

(3)在离散化过程中Δx为空间维度上的单位1，且取Δt为时间维度上的单位1，在任意时刻的溶质运移的对流弥散方程的表达式为：

$\frac{\∂P(x,t)}{\∂t}=P(x-1,t)F+P(x,t)R+P(x+1,t)B---\left(1\right)$

其中，F＝λ，B＝μ,R＝-λ-μ；即在一维尺度范围内，粒子向左跳动的概率是μ，向右跳动的概率是λ，根据公式计算得到λ和μ，其中，v^*为对流速度；

(4)将0时刻的数据带入方程式(1)，进行迭代，即可得到后面任意时刻的任意点的浓度值，从而实现对溶质运移的过程模拟。

2.如权利要求1所述的模拟地下水一维溶质运移过程的方法，其特征在于，所述步骤(1)中释放过程是瞬时的，或是持续的。

3.如权利要求1所述的模拟地下水一维溶质运移过程的方法，其特征在于，所述v^*的取值范围为0.001米/秒到1米/秒。