基于SPH-SWE方法对河流保守型污染物扩散的模拟

田丽蓉，余金玲，吴玉帅，顾声龙

(青海大学水利电力学院，青海 西宁 810016)

光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)是一种基于拉格朗日描述下的无网格方法，将连续介质离散成为许多粒子，通过追踪离散粒子的运动轨迹来描述连续介质的运动过程。SPH方法最初由Lucy等[1]和Gingold等[2]在研究天体物理问题中提出的一种方法。它已经成功应用到水动力学领域，而且是解决水力学领域诸多问题较成熟的数值模拟工具。SPH方法在处理大变形，自由表面流[3]和多相流[4]等问题上具有明显的数值优势。近年来，基于浅水波方程的SPH方法也得到了很大的发展[5]。Cleary等[6]将SPH方法应用于热传导问题，后续研究大多是基于此文献提出的SPH扩散模型来模拟，研究不同环境下污染物的扩散。Zhu等[7]研究了溶质在多孔介质中的传输。饶登宇等[8]建立描述孔隙水流动的SPH水动力模型和描述溶质分子扩散的扩散模型，得出弥散度与流速变异系数、迂曲度、迂曲路径差以及不均匀系数大致呈正相关，与孔隙率呈负相关。沈嘉渊[9]则使用SPH方法和一阶精度的欧拉向前方法模拟污染物扩散问题，通过不同时间积分方法的对比分析，得出精度最高的时间积分法。然而，在浅水环境下模拟污染物扩散的研究较少，采用SPH-SWE(Shallow Water Equation，SWE)方法研究污染物扩散，利用无网格性质求解浅水中污染物的输运过程比网格法求解SWE更具有优势[10]。由于污染物输运方程的二阶性质，传统的网格方法在模拟过程中存在网格畸形等问题。无网格方法中SPH法的发展迅速，在流体力学领域具有显著数值优势。它不依赖网格，在计算中以质点描绘整个模拟区域，因此在模拟污染物扩散时具有独特的优越性。本研究使用开源代码SWE-SPHysics(http：//www.sphysics.org)求解水流运动方程，并对其二次开发，近似离散求解污染物的输运方程，设置一维/二维均匀流在不同扩散系数下的污染物案例，并进行模拟，在验证结果合理的前提下，对实际地形污染物扩散的过程进行预测模拟。

1 模拟方法

1.1 控制方程

忽略地球自转引起的科氏力和流体粘度，拉格朗日描述下的SWE方程形式如下：

(1)

(2)

(3)

式中：ρ为密度，v为速度，g为重力加速度(g=9.81 m/s2)，b为河床粒子海拔，c为污染物浓度，D为污染物的浓度扩散系数,Sf=nv|v|/dw4/3，其中n为曼宁糙率系数。

1.2 水深求解

采用变密度法和变光滑长度法求解水深，利用牛顿迭代法求解流体粒子密度和水深。流体粒子密度ρl：

ρl=ρw/d

(4)

(5)

(6)

式中：ρw表示流体体积密度(ρw=1 000 kg/m3)，m为流体粒子质量，下标“i”和“j”分别为i粒子和j粒子，W为核函数，d为水深，ρ0和h分别为初始密度和光滑长度，“dm”为维度，下标“0”表示初始状态。

1.3 速度求解

浅水环境下，动量方程的SPH离散格式如下：

(7)

式中：ti=Ti/mi，bi为i河床粒子梯度，ki=(b)。

(8)

(9)

修正矩阵计算公式：

(10)

1.4 浓度求解

浅水环境下，污染物输运方程的SPH离散格式如下：

(11)

式中：D为污染物的浓度扩散系数，D也可是偏导数Dxx、Dxy、Dyx和Dyy，变扩散系数公式：

DL=6.09×u*×h

(12)

DT=0.6×u*×h

(13)

(14)

(15)

(16)

2 模拟工况设计

2.1 一维均匀流模拟

设置长10 000 m的矩形水槽，d=1 m，粒子间距△x=1 m，共10 000个流体粒子。模拟D为0、5、50 m2/s时污染物的扩散过程。出入流边界，vin/out=2 m/s，din/out=1 m，浓度为0 kg/m3，下标“in/out”表示入流/出流，初始流体充满计算域，v0=2 m/s，d0=1 m，c0设置如下：

(17)

一维分布源的解析解：

(18)

式中：c0为初始浓度，x1为浓度投放位置。

为验证SPH-SWE扩散模型在不同扩散系数的工况计算结果是否合理，设计模拟一维均匀流连续分布源。主要参数包括：其污染源类型都是分布源，在2 500～3 500 m处投放污染物浓度为1 kg/m3，扩散系数分别取0、5、50 m2/s的3种工况。

2.2 经典溃坝案例

设置经典干、湿溃坝案例，不考虑扩散影响。干、湿河床的计算区域长均为2 000 m，△x=10 m，上游d0和c0分别为10 m和0.7 kg/m3，下游d0和c0分别为5 m和0.5 kg/m3，下游无水，v0=0 m/s，模拟瞬时溃坝后污染物浓度扩散过程，与解析解[11]对比验证。

(19)

2.3 二维均匀流下的污染物扩散模拟

设置二维均匀流的验证工况，计算域1 200 m×400 m，△x=1 m，初始状态时，流体粒子充满计算域，坡度为0.001，初始vi和d分别为2.9 m/s和5 m，vin、din和cin分别为2.9 m/s、5 m和0 kg/m3，分别模拟D为0，10 000 m2/s时污染物浓度的分布情况。

2.4 实际地形模拟

选取位于青海省互助县南门峡河上南门峡水库作为研究区域，库区海拔高程在2 700 m以上，占地面积约为218 km2，地形图如图1所示。郑仙佩[12]对南门峡水库区域进行了溃坝模拟研究，取△x为10、15、20 m进行模拟，并对模拟结果进行了收敛性分析，得出的3种结果基本一致，验证了SPH-SWE模型的收敛性和结果的可靠性。

图1 南门峡水库Fig.1 Nanmenxia reservoir

模拟污染物扩散系数取0 m2/s和变扩散系数时，南门峡水库瞬间溃坝后污染物的扩散过程。计算区域为7 380 m×10 800 m，比降0.012，糙率0.02。河床△x=30 m，共89 167个河床粒子，流体△x=20 m，共2 664个流体粒子，污染物坐标位置(X：2 450 m，Y：9 255 m)，污染物投放半径为100 m，浓度为1 kg/m3，其它区域co为0 kg/m3，模拟时长600 s。

3 结果分析

3.1 一维均匀流模拟结果分析

图2为一维瞬时点源的解析解和模拟值对比图。图2a为当t=2 000 s时，计算域内流体粒子的水深-位移图和速度-位移图，水深和速度分别为常数1 m和2 m/s。图2b是Dx=0 m2/s，t为0、1 300、2 600 s时，解析解与模拟结果吻合。图2c和图2d分别是Dx取5 m2/s和50 m2/s，t为0、1 300、2 600 s，模拟值与解析解一致，上述3种工况均采用SPH-SWE扩散模型模拟污染物的扩散过程，其结果合理，验证了扩散系数大小对模拟结果无影响。

图2 一维瞬时点源的解析解和模拟值对比图Fig.2 Comparison of analytical and simulation results with 1D instantaneous point source

3.2 经典溃坝模拟结果分析

图3为t=50 s时干河床溃坝后水深和浓度分布图，图3a为干河床溃坝计算域内的水深变化趋势图，图3b为污染物浓度为常数0.7 kg/m3。图4为t=50 s时湿河床溃坝后水深和浓度分布图，采用MATLAB求得解析解，将解析解与模拟值进行对比。图4a中，水深的数值结果接近解析解，在不连续处未出现数值震荡，证明了水流计算结果的合理性。图4b为污染物浓度的模拟值和解析解一致，验证了SPH-SWE模型具有较准确的计算结果。

图3 干河床溃坝后水深和浓度分布图Fig.3 Distribution of water depth and concentration after the dam break of dry bed

图4 湿河床溃坝后水深和浓度分布图Fig.4 Distribution of water depth and concentration after the dam break of wet bed

3.3 二维均匀流模拟结果分析

图5为二维均匀流中D取0 m2/s时，仅考虑移流作用下污染物随时间变化的浓度分布图。由图5可知，当t为0、50、120、200、300 s时，污染物扩散到达最右端的距离是均匀水流速度与时间的乘积，分别对应200、345、548、786、1 076 m，验证了移流扩散过程中污染物浓度随时间扩散的模拟结果具有可靠性。

图5 移流作用下二维均匀流中污染物浓度分布图Fig.5 Distribution of pillutant concentration in two-dimensional uniform flow under the effects of advection

图6为D取10 000 m2/s时考虑移流和扩散作用下污染物随时间变化的浓度分布图。由图6可知，当t为0、50、120、200、300 s时的污染物浓度分布，与图5污染物扩散源的中心位置保持一致，由于扩散效应导致污染物的扩散范围比图5更大。根据模拟结果的合理性，验证了SPH-SWE扩散模型在模拟二维水流仍然具有很好的稳定性，而且不需考虑移流扩散项的求解，只需对其余部分离散求解，结果表明SPH-SWE方法在污染物扩散模拟方面具有很好的优势。

图6 移流扩散作用下二维均匀流中污染物浓度分布图Fig.6 Distribution of pillutant concentration in two-dimensional uniform flow under the effects of advection and diffusion

3.4 实际地形预测结果分析

图7分别是当t为600 s时，南门峡水库溃坝后在D取0 m2/s和变扩散系数条件下，污染物扩散的浓度分布图。郑仙佩[12]对南门峡水库溃坝的水流运动过程进行了细致的讨论，且验证粒子间距取20 m的方案是可行的，这里不做重复讨论，着重考虑溃坝时污染物扩散浓度分布情况。

图7 南门峡水库溃坝后浓度分布图Fig.7 Distribution of dam break concentration in Nanmenxia reservoir

图7a为当D=0 m2/s，时间取600 s时污染物在实际溃坝后的浓度分布图。由图7a可知，污染物的迁移过程主要与水流流速有着很大的关系，污染物主要集中在水流的中间位置。但仅考虑污染物的移流过程，浓度会出现局部不连续，实际水流不会出现此情况。由于此模拟工况的前提假设是不考虑扩散效应，污染物在输运过程中只与流体粒子的运动轨迹有关。图7b为在变扩散系数下时间取600 s时，南门峡水库溃坝后污染物浓度的分布图，由图7b可知，在沿河道水流间位置污染物浓度达到峰值0.075 6 kg/m3。其次，浓度值依次向上、下游递减，不会出现图7a中浓度不连续的情况，浓度未出现数值奇异点，符合实际河流中污染物扩散规律。图7b污染物的变扩散系数与地形起伏、水流速度和水深有关，接近实际宽浅河流中污染物扩散的过程，验证了SPH-SWE扩散模型在复杂水流环境中具有很好的稳定性，得出SPH-SWE方法在污染物输运过程的预测具有很好的研究价值，可较好预测可溶保守型污染物在河流中的浓度分布情况。

4 讨论与结论

污染物在水体中的迁移扩散规律是水环境分析的基础内容，对研究河流中污染物扩散规律具有重要意义。刘圣勇[13]采用有限差分法对一维河流污染物扩散进行模拟。许媛媛[14]采用有限体积法对二维浅水方程和物质输运方程进行求解，该模型采用高精度的迎风重构技术处理网格界面的污染物输运，虽然计算精度有所提高，但是计算效率较低和求解格式较复杂。本文采用SPH-SWE法对河流中可溶保守型污染物扩散进行模拟，由于SPH方法具有拉格朗日的特性，在模拟污染物对流扩散过程具有优势。通过设计一维和二维均匀流工况，模拟保守型溶质污染在不同扩散系数下的迁移扩散过程，数值模拟结果和解析解一致；同时选取经典溃坝案例进行验证，模拟结果与解析解高度吻合，表明了SPH-SWE法对河流中污染物的迁移扩散模拟研究具有很大潜力；通过预测南门峡水库溃坝后污染物浓度分布情况，得到在沿河道水流中间位置污染物浓度达到峰值0.075 6 kg/m3，可预测污染物在实际地形中的浓度变化趋势和规律。本研究主要有以下几点结论：

(1)SPH-SWE扩散模型模拟一维、二维均匀流和经典溃坝工况，避免了污染物浓度奇异值的出现，验证了该模型具有很好的稳定性。

(2)验证了出入流边界条件下，SPH-SWE扩散模型可以很好地模拟可溶保守型污染物的迁移扩散过程。

(3)实际地形预测污染物模拟结果合理，表明了SPH-SWE法可用于实际地形中河流污染物扩散规律的研究，可较好预测可溶保守型污染物在河流中的浓度分布情况。