基于粒子追踪的孔隙地下水污染物三维运移模拟方法

摘要

一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法，包括：设定粒子总数、每个粒子的初始坐标、释放时间，设定计算时间步长；在任意时刻、任意位置释放粒子；对每个被释放的粒子，计算其在一个时间步长后的位移；根据计算得到的粒子的新坐标，判断该粒子是否位于计算网格区域内，若已跳出，则通过边界控制使其回到计算区域边界上；重复上述步骤，以此迭代，计算得出每个时间节点上每个粒子的坐标；通过统计特定区域内粒子数量，即可得到该区域的地下水污染物浓度分布。与传统方法相比，本发明的方法可大大降低计算成本，且利用该方法得到的模拟结果与对流弥散方程的解析解有很好的匹配。

说明书

技术领域

本发明涉及地下水数值模拟计算领域，特别是涉及一种基于粒子追踪 的孔隙地下水污染物三维运移模拟方法。

背景技术

孔隙地下水作为北方重要的饮用水源，其水质直接关系到饮水安全。 地下水埋藏在地下，其污染具有很强的复杂性和隐蔽性，对地下水污染的 掌握需要借助数值模拟方法，通过模型的高效运算，可以节约大量经济和 时间成本，使管理决策人员迅速掌握地下水中污染物浓度变化情况。

地下水数值模拟主要包括地下水水流模拟和溶质运移模拟，前者数值 求解地下水水流方程，后者数值求解对流-弥散方程。描述孔隙地下水溶质 运移的对流-弥散方程可表述为：

$\frac{\partial C}{\partial t}=-▿\left(vC\right)+▿(D▿C);$

其中C表示地下水中溶质的浓度，v表示地下水流动孔隙速度，D表 示水力弥散系数，t表示时间，为微分算子。

孔隙地下水溶质运移模拟求解方法常用欧拉法、拉格朗日法，以及二 者的结合。欧拉法以空间中固定坐标系作为参照系，常见的有限差分法和 有限单元法属于欧拉法。应用有限元和有限差分方法进行地下水溶质迁移 模拟，有两个固有的缺陷，一是当网格Peclet数较大，即对流项强于弥散 项时，容易受到数值弥散的影响；二是当模拟区域范围较大时，因计算网 格较多，计算成本比较昂贵。

拉格朗日法通过质点追踪，以运动坐标系作为参照系，常见的拉格朗 日法有粒子追踪法等。运动坐标系容易造成数值不稳定；另外在质点追踪 过程中对质点速度的连续性有较高要求，否则由于速度插值也容易造成局 部的质量不守恒。

欧拉-拉格朗日混合法用拉格朗日法解决溶质运移中的对流问题，用欧 拉法解决弥散问题，结合了二者的优点，却也同时存在缺点，且计算耗时。

发明内容

有鉴于此，本发明的目的在于提出一种孔隙地下水中污染物三维运移 模拟方法，该数值模拟方法是基于拉格朗日法中的粒子追踪方法，以解决 模拟过程中的数值弥散问题，模拟得到的数值解可以完全拟合对流-弥散方 程的解析解，且与有限元法相比可大大节约计算成本。

为实现本发明的上述目的，本发明提出了一种基于粒子追踪的孔隙地 下水污染物三维运移模拟方法，包括以下步骤：

步骤S1：设定粒子总数、每个粒子的初始坐标、释放时间，设定计算 时间步长；

步骤S2：在任意时刻、任意位置释放粒子；

步骤S3：对每个被释放的粒子，计算其在一个时间步长后的位移，粒 子在地下水流场中的运动由以下方程控制：

$x_{t+\Delta t}=x_t+(v_x(x_t,y_t,z_t,t)+\frac{\partial D_{xx}}{\partial x}+\frac{\partial D_{xy}}{\partial y}+\frac{\partial D_{xz}}{\partial z})\Delta t+\sqrt{2D_{xx}\Delta t}{Z}_1+\sqrt{2D_{xy}\Delta t}{Z}_2+\sqrt{2D_{xz}\Delta t}{Z}_3$

$y_{t+\Delta t}=y_t+(v_y(x_t,y_t,z_t,t)+\frac{\partial D_{yx}}{\partial x}+\frac{\partial D_{yy}}{\partial y}+\frac{\partial D_{yz}}{\partial z})\Delta t+\sqrt{2D_{yx}\Delta t}{Z}_1+\sqrt{2D_{yy}\Delta t}{Z}_2+\sqrt{2D_{yz}\Delta t}{Z}_3$

$z_{t+\Delta t}=z_t+(v_z(x_t,y_t,z_t,t)+\frac{\partial D_{zx}}{\partial x}+\frac{\partial D_{zy}}{\partial y}+\frac{\partial D_{zz}}{\partial z})\Delta t+\sqrt{2D_{zx}\Delta t}{Z}_1+\sqrt{2D_{zy}\Delta t}{Z}_2+\sqrt{2D_{zz}\Delta t}{Z}_3;$

其中x、y、z表示粒子的空间坐标，v表示粒子的对流运动速度，Δt 表示时间步长，D表示水力弥散系数，Z表示介于0到1之间的随机数， 这样就根据粒子在t时刻的坐标计算出其在一个时间步长后，也即t+Δt 时刻的新坐标；

步骤S4：根据计算得到的粒子的新坐标，判断该粒子是否位于计算网 格区域内，若已跳出，则通过边界控制使其回到计算区域边界上；

步骤S5：重复步骤S2-S4，以此迭代，计算得出每个时间节点上每个 粒子的坐标；

步骤S6：通过统计特定区域内粒子数量，即可得到该区域的地下水污 染物浓度分布。

基于上述技术方案可知，本发明的孔隙地下水中污染物三维运移模拟 方法，将地下水中的溶质抽象为大量粒子，以粒子的运动来模拟溶质在孔 隙介质中的运移，其中以粒子的有序运动来刻画溶质因对流引起的迁移过 程，以粒子的随机位移来刻画溶质的弥散过程，因此可以精细刻画溶质在 地下水中的运移过程。传统的有限元法在求解对流-弥散方程时，需要求解 的代数方程个数与计算网格节点数成正比，本发明的方法基于统计物理学 中的随机行走粒子追踪，不直接求解对流-弥散方程，而是求解粒子的位移， 计算量仅与粒子个数有关，因此在三维情景下与有限元法相比可大大降低 计算成本，且运算结果与对流-弥散方程的解析解完全拟合。本发明的方法， 粒子不会凭空出现或消失，粒子总数量在输入和输出及整个过程中是不变 的，因此该方法从本质上是质量守恒的，可消除计算过程中的数值弥散。

附图说明

图1是瞬时释放污染物在三维均质孔隙介质中的污染羽分布；

图2是粒子到达计算网格区域边界时的边界控制效果示意图；

图3是瞬时释放污染物在三维均质孔隙介质中的突破曲线。

具体实施方式

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

本发明公开了一种基于粒子追踪的孔隙地下水污染物三维运移模拟 方法，包括以下步骤：

步骤S1：设定粒子总数、每个粒子的初始坐标、释放时间，设定计算 时间步长；

步骤S2：在任意时刻、任意位置释放粒子；

步骤S3：对每个被释放的粒子，计算其在一个时间步长后的位移，粒 子在地下水流场中的运动由以下方程控制：

$x_{t+\Delta t}=x_t+(v_x(x_t,y_t,z_t,t)+\frac{\partial D_{xx}}{\partial x}+\frac{\partial D_{xy}}{\partial y}+\frac{\partial D_{xz}}{\partial z})\Delta t+\sqrt{2D_{xx}\Delta t}{Z}_1+\sqrt{2D_{xy}\Delta t}{Z}_2+\sqrt{2D_{xz}\Delta t}{Z}_3$

$y_{t+\Delta t}=y_t+(v_y(x_t,y_t,z_t,t)+\frac{\partial D_{yx}}{\partial x}+\frac{\partial D_{yy}}{\partial y}+\frac{\partial D_{yz}}{\partial z})\Delta t+\sqrt{2D_{yx}\Delta t}{Z}_1+\sqrt{2D_{yy}\Delta t}{Z}_2+\sqrt{2D_{yz}\Delta t}{Z}_3$

$z_{t+\Delta t}=z_t+(v_z(x_t,y_t,z_t,t)+\frac{\partial D_{zx}}{\partial x}+\frac{\partial D_{zy}}{\partial y}+\frac{\partial D_{zz}}{\partial z})\Delta t+\sqrt{2D_{zx}\Delta t}{Z}_1+\sqrt{2D_{zy}\Delta t}{Z}_2+\sqrt{2D_{zz}\Delta t}{Z}_3;$

其中x、y、z表示粒子的空间坐标，v表示粒子的对流运动速度，Δt 表示时间步长，D表示水力弥散系数，Z表示介于0到1之间的随机数， 这样就根据粒子在t时刻的坐标计算出其在一个时间步长后，也即t+Δt 时刻的新坐标；

步骤S4：根据计算得到的粒子的新坐标，判断粒子是否位于计算网格 区域内，若已跳出，则通过边界控制使其回到计算区域边界上；

步骤S5：重复步骤S2-S4，以此迭代，计算得出每个时间节点上每个 粒子的坐标；

步骤S6：通过统计特定区域内粒子数量，即可得到该区域的地下水污 染物浓度分布。

作为优选，在步骤S1中，设定粒子初始坐标时，允许多个粒子的初 始坐标完全相同；粒子释放时间表示该粒子初次被释放的时间，应为时间 步长的整数倍；时间步长需满足库朗数(Courantnumber)小于1。

作为优选，在步骤S2中，对任意一个粒子，若当前步的计算时间大 于或等于步骤S1中设定的粒子释放时间，则该粒子被释放；若在某一坐 标位置所有粒子的释放时间相同，则该释放过程是瞬时的，若粒子释放时 间不同且构成一个连续时段，则该释放过程是持续的。

作为优选，在步骤S3中，粒子的对流运动速度v由地下水流场计算 所得，并通过粒子所在单元节点插值获得；Z由计算程序中的随机函数产 生；对流弥散系数D由下式计算获得：

$D_{ij}={\alpha }_T\left|v\right|{\delta }_{ij}+({\alpha }_L-{\alpha }_T)\frac{v_i{v}_j}{\left|v\right|}+D_{ii}^d;$

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

作为优选，在步骤S4中，通过粒子位置坐标与计算网格坐标进行比 较，确定粒子是否位于计算网格之内；边界控制的实现是通过计算粒子在 该时间步长内的位移轨迹与计算网格边界的交点，将交点坐标作为粒子的 新坐标。

作为优选，在步骤S6中，特定区域内粒子数量与被释放粒子总量的 比值，即为该区域内地下水污染物的相对浓度。

作为优选，以一定数量的粒子代表地下水中污染物(溶质)的一定浓 度，以粒子的位置代表地下水中污染物(溶质)的空间分布。其中，该粒 子没有体积和质量；不同的粒子可以按组分类，属于不同组别的粒子可以 表示相同类型，也可以表示不同类型的地下水污染物。

下面通过本发明的一个优选实施例对本发明做进一步说明。

如图1所示，选取正方体结构的计算区域，该区域为均质孔隙介质， 地下水流向为对角线方向。在该区域上游即正方体一个顶点附近瞬时释放 粒子，模拟地下水中污染物在该区域的分布。具体步骤如下：

步骤S1-S2：设定粒子总数为500，每个粒子的初始坐标为(0.01，0.01， 0.01)，在初始时刻瞬时释放；设定时间步长为一个粒子通过一个计算网格 长度所需时间的1/10；

步骤S3：计算一个时间步长后每个粒子的位移，根据粒子在t时刻的 坐标计算出其在一个时间步长后，也即t+Δt时刻的新坐标；

步骤S4：根据计算得到的粒子的新坐标，判断粒子是否位于计算网格 区域内，若已跳出，则通过边界控制使其回到计算区域边界，如图2所示；

步骤S5：重复步骤S2-S4，以此迭代，计算得出每个时间节点上每个 粒子的坐标；

步骤S6：通过统计特定区域内粒子数量，该数量与500的比值即可认 为是任意位置的溶质浓度。

由该区域下游即正方体对角顶点处的污染物浓度变化，获取突破曲线， 模拟结果与解析解对比如图3所示，其中红色曲线为解析解，蓝色曲线为 模拟结果。可见该方法获得的模拟结果与解析解完全拟合。

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