CN106886682A - 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 - Google Patents

Abstract

一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法，包括：建立裂隙介质概念模型，将裂隙结构离散成二维计算网格，所述计算网格属性包括隙宽；建立裂隙介质地下水水流数学模型，计算得到模拟地下水水流场分布，并与实验观测值进行对比验证，获得接近裂隙真实分布的模拟裂隙结构；针对所述模拟裂隙结构及其对应的模拟地下水水流场分布，用随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，获取溶质的分布。

Description

用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法

技术领域

本发明涉及裂隙地下水中渗流与溶质运移的数值模拟领域，特别是关于一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法。

背景技术

含水介质广义上可分为孔隙介质、裂隙介质与岩溶介质。裂隙介质本身是透水性非常微弱的基岩，很多卫生垃圾填埋、放射性核废料处置等地下工程选择裂隙介质作为天然防渗层。随着此类工程的发展，其对地下水的影响受到越来越多的重视，裂隙介质中地下水水流和溶质运移研究也因此成为水环境管理和保护方面一个重要领域。

裂隙介质中地下水水流和溶质运移的主要通道是裂隙，相较于孔隙介质，对裂隙介质中地下水水流和溶质运移的数值模拟具有更大的难度。一方面裂隙介质具有强烈的非均质性，从现场获取的裂隙数据很难完整且准确的表达裂隙的真实分布。另一方面，在数值模拟中裂隙表面粗糙度、裂隙的分布和连通性、裂隙与岩体的相互作用等因素都会对地下水水流和溶质运移过程造成影响。

通常对于裂隙介质的数值模拟有单裂隙渗流概念模型和裂隙岩体渗流概念模型，后者包括等效连续体模型和离散裂隙网格模型，以及混合模型等。单裂隙渗流概念模型是假设岩体被一条裂隙切割，可将岩体概化成两个表面光滑且平行的平板，平板间隔代表裂隙，间隔大小代表隙宽。

研究单裂隙的水力特征是研究岩体裂隙渗流的基础。天然岩体裂隙均为粗糙裂隙，很难满足单裂隙渗流概念模型中光滑平板的假定。考虑裂隙表面粗糙度对裂隙水流的影响，就需要对单裂隙渗流概念模型进行修正。有关裂隙的参数估计是裂隙介质数值模拟的关键，特别是渗透系数，对表征地下水渗流及污染物浓度分布和运移状况至关重要。

在裂隙介质溶质运移问题中，除了要考虑裂隙中的渗流作用，还需要考虑水力弥散和分子扩散过程。传统的有限元和有限差分方法需要较高的计算成本，在裂隙介质污染物运移数值模拟中难以发挥作用。

发明内容

鉴于现有方案存在的问题，为了克服上述现有技术方案的不足，本发明提出了一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法。

根据本发明的一个方面，提供了一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法，包括：建立裂隙介质概念模型，将裂隙结构离散成二维计算网格，所述计算网格属性包括隙宽；建立裂隙介质地下水水流数学模型，计算得到模拟地下水水流场分布，并与实验观测值进行对比验证，获得接近裂隙真实分布的模拟裂隙结构；针对所述模拟裂隙结构及其对应的模拟地下水水流场分布，用随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，获取溶质的分布。

从上述技术方案可以看出，本发明具有以下有益效果：

考虑单裂隙介质的非均质性，采用地质统计学方法获取模拟隙宽分布，并通过重复用实验观测值验证模拟得到的水流场分布，获取最佳裂隙结构，相较于单裂隙渗流概念模型更接近岩体和裂隙的真实分布；

基于随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，采用粒子模拟污染物，通过计算粒子在水流场中的位移获得粒子在不同时间点的位置，采用统计的方式获得污染物分布，相较于有限元和有限差分法可以更精细的刻画污染物的运移。

附图说明

图1为本发明实施例提供用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法的流程图；

图2为图1中步骤S2中的具体流程图；

图3为图1中步骤S3中的具体流程图。

具体实施方式

本发明某些实施例于后方将参照所附附图做更全面性地描述，其中一些但并非全部的实施例将被示出。实际上，本发明的各种实施例可以许多不同形式实现，而不应被解释为限于此数所阐述的实施例；相对地，提供这些实施例使得本发明满足适用的法律要求。

在本说明书中，下述用于描述本发明原理的各种实施例只是说明，不应该以任何方式解释为限制发明的范围。参照附图的下述描述用于帮助全面理解由权利要求及其等同物限定的本发明的示例性实施例。下述描述包括多种具体细节来帮助理解，但这些细节应认为仅仅是示例性的。因此，本领域普通技术人员应认识到，在不悖离本发明的范围和精神的情况下，可以对本文中描述的实施例进行多种改变和修改。此外，为了清楚和简洁起见，省略了公知功能和结构的描述。此外，贯穿附图，相同附图标记用于相似功能和操作。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明实施例提供了一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法，采用地质统计方法获取裂隙结构，更接近岩体和裂隙的真实分布，具体地，由地质统计学的方法来获得密度分布函数，以此来表征单裂隙介质计算网格的隙宽值，根据隙宽计算渗透率，从而求解单裂隙介质地下水渗流过程，基于随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，具体地，将地下水中的溶质抽象为大量粒子，以粒子的运动来模拟溶质在单裂隙介质中的运移，其中以粒子的有序运动来刻画溶质因对流引起的迁移过程，以粒子的随机位移来刻画溶质的分子扩散过程，可以更精细的刻画污染物的分布。

图1为本发明实施例提供的一种单裂隙中溶质运移数值的模拟方法的流程图，如图1所示，该模拟方法具体包括：

步骤S1：建立单裂隙介质概念模型。

具体地，建立裂隙结构离散成二维计算网格的裂隙介质模型，选取四边形平面结构的计算区域代表单裂隙介质，将裂隙结构离散成二维计算网格，该网格可以是规则网格，也可以是不规则网格，将隙宽作为每个网格的属性，而不需要作为网格的维度。

步骤S2：采用数值模拟方法获取地下水水流场分布，并进行验证，获得接近岩体和裂隙的真实分布的模拟裂隙结构。

步骤S2具体包括，如图2所示：

步骤S201：采用一个密度分布函数来表征计算网格的隙宽b，该密度分布函数由地质统计学的方法来获得。

具体地，用统计学方法生成一组符合高斯分布的数据，数据个数等于步骤S1中生成的二维计算网格个数，以生成数据的值代表隙宽，用地质统计学方法生成隙宽分布函数时，应结合实验方法，以裂隙实际测量值作为参考。

步骤S202：计算每个计算网格的渗透率，求解地下水水流方程，获得模拟二维地下水流场分布；

具体地，每个计算网格的渗透率k由下式计算获得：

其中k代表渗透率，b代表计算网格的隙宽，

该地下水水流方程为：

其中h为地下水水头，x，y为空间坐标，S_S为单位储水量，K为渗透系数，t为时间，

渗透系数K与渗透率k的关系：

ρ为水的密度，g为重力加速度，μ为水动力黏度。

步骤S203：以实验观测值对步骤S202生成的模拟二维地下水流场分布进行验证；

具体地，由于裂隙中的实际流场很难获取，可以用流入边界和流出边界的地下水水头来做验证，所述验证采用选取流入与流出边界的水头，将数值模拟结果与实验观测值进行比较，若误差值超过可接受的限度，则去除该组数据，若误差值在可接受的限度内，则保存该组数据及相应的误差值。

步骤S204：重复步骤S201-S203，获得最优组数据，获得最优模拟裂隙结构。

具体地，多次重复步骤S201-S205，获得多组数据及相应的误差值，选取误差值最小的一组数据，可以认为该组数据对应的隙宽分布能够有效代表单裂隙结构的隙宽分布，进而获得模拟裂隙结构，并获得该隙宽分布下的模拟地下水流场。

步骤S3，针对所述模拟裂隙结构及其隙宽分布下的地下水流场，用随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，更精细的刻画污染物的分布。

步骤S3具体包括，如图3所示：

步骤S301：以粒子代表污染物，设定粒子总数、每个粒子的初始坐标、每个粒子的释放时间等参数；

具体地，利用步骤S2中获得的最优模拟裂隙结构，以随机行走粒子追踪方法进行单裂隙介质中溶质运移模拟，以一定量的粒子代表污染物，设定粒子总数、每个粒子的初始坐标、释放时间等参数，粒子在计算初始时刻为静止状态，每个粒子是相对独立的，设定粒子初始坐标时，允许多个粒子的初始坐标完全相同；

比如，在四边形平面结构的流入边界设置一定数量的均匀分布的粒子以代表污染物，例如10000个，设定所有粒子的释放时间为t＝0s。

步骤S302：计算出每一个时间节点上每个粒子的坐标。

针对每一个粒子，当计算时间大于或等于粒子释放时间时，该粒子被释放，可随地下水迁移，并基于分子扩散进行随机行走；粒子的释放过程可以是瞬时的，也可以是持续的。

对于每个被释放的粒子，计算其在一个时间步长后的位移，该粒子在二维地下水流场中的运动由以下公式控制：

其中x，y代表粒子的空间坐标，v代表粒子的对流运动速度，t代表任意时刻，Δt代表时间步长，D代表对流弥散系数，Z₁、Z₂是一个介于0到1之间的随机数，

其中粒子的对流运动速度v等于粒子所在位置的地下水流速，对流弥散系数D由下式计算获得：

其中δ_ij为克罗内克符号，α_L为纵向弥散系数，α_T为横向弥散系数，为分子扩散系数，v_i为i方向平均渗流速度。

根据上述公式计算出每一个时间节点上每个粒子的坐标。

步骤S303：通过统计特定区域内粒子数量，即可求出任意位置的溶质浓度。

具体地，统计某一时间节点特定区域内粒子数量，特定区域内粒子数量与被释放粒子总量的比值，即为该区域内地下水污染物的相对浓度。比如，流出边界上的粒子数量与之前设定的粒子的数量的比值即为流出边界上的溶质浓度。

本实施例中提出的方法以统计学方法获取裂隙隙宽，解决了实际工程应用中裂隙隙宽难以测量的问题；同时以密度分布函数表征隙宽值，不同于单裂隙渗流概念模型中的光滑平板假设，更加接近实际；随机行走粒子追踪方法本身具有统计学特性，加之粒子的随机行走特性，可以模拟水力弥散和分子扩散过程，相较于有限元和有限差分法更加适合裂隙介质中污染物运移模拟。

还需要说明的是，本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值。

需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法，其特征在于，包括：

步骤S1：建立裂隙介质概念模型，将裂隙结构离散成二维计算网格，所述计算网格属性包括隙宽；

步骤S2：建立裂隙介质地下水水流数学模型，计算得到模拟地下水水流场分布，并与实验观测值进行对比验证，获得接近裂隙真实分布的模拟裂隙结构；

步骤S3：针对所述模拟裂隙结构及其对应的模拟地下水水流场分布，用随机行走粒子追踪方法对单裂隙介质中溶质运移进行模拟，获取溶质的分布。

2.根据权利要求1所述的随机行走粒子追踪方法，其中，步骤S2包括：

步骤S201：采用一个密度分布函数来表征所述计算网格的隙宽；

步骤S202：计算每个计算网格的渗透率，求解地下水水流方程，获得模拟二维地下水流场分布；

步骤S203：以实验观测值对所述模拟二维地下水流场分布进行验证；以及

步骤S204：重复步骤S201-S203，选取与实验观察值误差最小的模拟二维地下水流场分布，采用其对应的模拟裂隙结构作为最接近岩体和裂隙的真实分布的模拟裂隙结构。

3.根据权利要求2所述的随机行走粒子追踪方法，其中，所述密度分布函数由地质统计学方法获得。

4.根据权利要求2所述的随机行走粒子追踪方法，其中，所述渗透率采由下式计算获得：

$k=\frac{b^2}{12}$

其中k代表渗透率，b代表计算网格的隙宽。

5.根据权利要求2所述的随机行走粒子追踪方法，其中，所述地下

水水流方程为：

$S_S\frac{\∂h}{\∂t}=K(\frac{{\∂}^{2}h}{\∂x^2}+\frac{{\∂}^{2}h}{\∂y^2})$

其中h为地下水水头，x，y为空间坐标，S_S为单位储水量，K为渗透系数，t为时间，渗透系数K与渗透率k的关系为：

$K=\frac{k\mathrm{\ρ}g}{\mathrm{\μ}}$

其中，ρ为水的密度，g为重力加速度，μ为水动力黏度。

6.根据权利要求2所述的随机行走粒子追踪方法，其中所述验证采用选取流入与流出边界的水头，将数值模拟结果与实验观测值进行比较。

7.根据权利要求1所述的随机行走粒子追踪方法，其中，步骤S3包括：

步骤S301：以粒子代表污染物，设定粒子总数、每个粒子的初始坐标、每个粒子的释放时间；

步骤S302：每个粒子在计算初始时刻为静止状态，当计算时间大于或等于该粒子释放时间时，该粒子被释放，计算出每一个时间节点上每个被释放粒子的坐标；

步骤S303：通过统计特定区域内粒子数量，求出任意位置的污染物浓度。

8.根据权利要求7所述的随机行走粒子追踪方法，其中，每个被释放粒子的坐标根据其原始坐标及其在一个时间步长内的位移计算得到，每个粒子在二维地下水流场中的运动由以下公式计算：

$x_{t+\mathrm{\Δ}t}=x_t+\left(v_x\right(x_t,y_t,t)+\frac{\∂D_{xx}}{\∂x}+\frac{\∂D_{xy}}{\∂y})\mathrm{\Δ}t+\sqrt{2D_{xx}\mathrm{\Δ}t}{Z}_1+\sqrt{2D_{xy}\mathrm{\Δ}t}{Z}_2$

$y_{t+\mathrm{\Δ}t}=y_t+\left(v_y\right(x_t,y_t,t)+\frac{\∂D_{yx}}{\∂x}+\frac{\∂D_{yy}}{\∂y})\mathrm{\Δ}t+\sqrt{2D_{yx}\mathrm{\Δ}t}{Z}_1+\sqrt{2D_{yy}\mathrm{\Δ}t}{Z}_2$

其中x，y代表粒子的空间坐标，v代表粒子的对流运动速度，t代表任意时刻，Δt代表时间步长，D代表对流弥散系数，Z₁、Z₂是介于0到1之间的随机数。

9.根据权利要求8所述的随机行走粒子追踪方法，其中，对流弥散系数D由下式计算获得：

$D_{ij}={\mathrm{\α}}_T\left|v\right|{\mathrm{\δ}}_{ij}+({\mathrm{\α}}_L-{\mathrm{\α}}_T)\frac{v_i{v}_j}{\left|v\right|}+D_{ij}^d$

其中δ_ij为克罗内克符号，α_L为纵向弥散系数，α_T为横向弥散系数，为分子扩散系数，v_i为i方向平均渗流速度，对流运动速度v等于粒子所在位置的地下水流速。

10.根据权利要求1所述的随机行走粒子追踪方法，其中，所述计算网格为规则网格或不规则网格。