射阳河口污染物输运数值模拟

魏 广

（南京扬子水利自动化技术开发总公司，江苏 南京 210012）

0 引言

射阳河是江苏省苏北地区排水入海的最大天然干河，随着经济发展和人口的不断增长，以及废水排放量逐年增加，射阳河流域水环境压力越来越大，近海水域水环境污染日益受到了人们的重视。射阳河口处设有闸门，开闸排放污染物的时间与污染物在河口近海地区的输运扩散规律和浓度分布特征有密切关系。为了有效地控制污染，合理地确定水环境容量及允许排放量，科学地制定污染物排放时间和过程，需要深入进行水动力状况及污染物在水环境中的输运迁移规律的研究，而数值模拟是解决上述问题的重要途径。

1 二维水流水质耦合数学模型建立及求解

1.1 控制方程

守恒型二维浅水方程[1]与对流扩散方程耦合的矢量表达式为：

式中：q = [ h, h u, h v, h Ci]T为守恒物理量；F (q) = [ h u,h u2+ g h2/ 2, h u v, h u Ci]T为 x 向通量；G (q) = [ h v, h u v,h v2+ g h2/ 2, h v Ci]T 为 y 向通量；b (q) = [0, g h (s0x- sfx),g h (s0y- sfx), Δ(DkΔ) h Ci)) - KCig h Ci+ Si]T为源汇项；x，y，t 分别为空间及时间坐标系；h 为水深；u 和 v 分别为 x 和 y 向沿水深积分平均流速分量；g 是重力加速度；s0x和 sfx分别是 x向的水底和摩阻坡度；s0y和 sfy分别是 y 向的水底和摩阻坡度；Ci为污染物沿水深积分的垂线平均浓度；本模型包括 8 个组分（COD，NH3-N，BOD，DO，温度，盐度，TP 及 TN）方程，可同时也可单独模拟。其中，KCi是污染物综合降解系数；Si是污染物源汇项，其内容随组分而变；Dk为污染物的扩散系数；Δ为梯度算子。

1.2 方程离散

在任意形状的单元 Ω，对式（1）进行积分离散求得[2]

式中：A 为单元 Ω 的面积；m 为单元边总数； Lj为第 j 单元边的长度。

b (q) = [0, A g h (s0x- sfx)，A g h (s0y- sfy)，∑Dk(Δ hCi)nLj+ Si- AKCihCi]T，q = T(φ) · q，T (φ) 和T (φ)-1为坐标旋转变换及逆变换矩阵[1-2]，其中 φ 为坐标轴旋转角度。

式 (2) 采用显格式时间离散，其求解归结为如何确定法向通量，可通过解算局部一维黎曼初值问题的外法向数值通量（记为ƒLR）而得到[3]。

1.3 水利设施处的通量计算

使用计算域内的水利设施作为模型的内边界，水利设施处的法向通量的计算不能应用黎曼近似解，而要采用与水利设施相应的出流公式[3]。其通量公式概况如下：

式中：Qm是水工建筑物出流；h 为水深；un，vτ分别为法向及切向流速，满足 un= u cos φ + v sin φ 和 vτ=v cos φ - u sin φ。模型中包含的水工建筑物为闸门。本模型中采用经验公式：Qg= C1g· C2g· Hg· Ng· Lg·式中：C1g为闸门系数；C2g是淹没系数；Hg为闸门启闭高度；Ng是总的闸门个数；Lg是每个闸门的长度；当尾水位高于闸门溢洪道堰顶高程时，Δ Z 为闸上水位和尾水位之差，反之为闸上水位和溢洪道堰顶高程之差。

1.4 动边界的处理

在河口地区，随着潮水涨落，浅滩的漫、露交替出现，高潮位时淹没，低潮位时露出。在进行该区域数值模拟时，模型必须具备捕捉这种动态过程的能力。动边界是水平计算域中有水与无水区域的界线。在该二维模型中，根据有限体积法特点，采用干湿动边界法，以水深为判别标准设定单元界面的类型，并运用相应的方法计算跨界面的法向通量，以保证水量平衡。在每一时间步判断单元界面2 侧单元的水深 hL，hR，当 hL和 hR均为 0 时，则通量为 0；当 hL和 hR有 1 个为 0 时，则结合 2 侧单元的底高程，交界面类型可能为固壁边界、跌水、漫流等形态的边界，根据不同类型可选择瞬时溃坝解析解、堰流公式等方式估算通量。当 hL和 hR均不为 0 时，则按正常的方法估算通量。在实际建模时，设定 1 个水深限值，以此来判断单元的干湿情况。若相邻单元的水深均小于水深限值，则认为该单元为干单元，不进行计算。

2 模型验证

2.1 射阳河口原型观测

为了考察河口二维水流水质耦合数学模型的可靠性、单元概化处理的合理性和模型参数选择的适当性，采用射阳河口 2006年8月24日 6:00 ～ 8 月25日 10:00（大潮）和 2007年1月3日 17:00 ～1月4日 19:00（大潮）的潮位、流向和平均流速等实测资料对模型进行验证。利用平均绝对误差、均方差和相对均方差，分析潮位、流速和流速的计算值与实测值的拟合程度。

本次计算范围全长约 74 km，射阳河口以北25 km，以南 49 km，河口至外海处 48 km，计算区域面积约 4070 km2。前处理工作实现了网格的剖分，对计算水域采用无结构任意四边形网格布置，以贴合黄海天然岸线边界。对河口附近区域采用网格局部加密技术。网格总数为 5943，节点总数为 6076，网格尺度变幅范围为 350 m×500 m～900 m×1000 m。

根据网格单元的大小，设定时间步长为 Δ t = 10 s。

2.2 水动力计算结果验证及分析

采用上述网格布置及边界条件等，得出水动力模拟计算结果。大潮潮位、流向和垂线流速计算值与实测值误差统计表如表1 所示。

表1 潮位计算值与实测值误差统计表

上述验证结果分析表明：射阳河口二维水流水质耦合数学模型对区域的概化处理合理，模型选用的参数基本重现了射阳河口区域的潮位及潮流变化过程，计算结果可信。能够较好地模拟射阳河口的水流运动，可在此基础上模拟污染物的输运扩散。

3 射阳河口污染物输运数值模拟

3.1 计算条件选取

3.2 浓度场模拟

计算步长：水质模拟时取 10 s。

边界条件：输入污染物 CODMn排放量 823.31 g/s。

初始值的选取：单元内初始污染物浓度取为0 mg/L。

针对射阳河口实际情况，高潮位时不开闸，低潮位时开闸，在排放量保持一定的情况下，计算模拟 2006年8月 31日落潮不同时段污染物浓度场分布，确定何时开闸排放污染物对环境影响最小。设计工况条件表如表2 所示。

表2 工况条件

通过数值模拟结果分析可以得出在整个潮流周期中：

1）工况 1 污染物变化幅度为1.0～7.0 mg/L，4.0 mg/L 的污染带面积在开闸 3 h 时最大，为 3.0 km2。由流场分布可知，退潮时，流场是由西北至东南向，污染物在潮流的作用下向东南方移动。涨憩过后刚开始退潮时，水流流速较小，污染物扩散速度偏慢，随着时间的推移，水流流速逐渐变大，污染物扩散速度也变快。由于前 3 h 内污染物是连续不断排放，此后关闸，因此开闸 3 h 时污染物浓度和影响范围大。污染物不断向东南方移动，但污染物浓度和影响范围逐渐变小。

2）工况 2 污染物变化幅度为1.0～6.0 mg/L，4.0 mg/L 的污染带面积在开闸 3 h 时最大，为 2.1 km2。开闸的 3 h 和关闸后的前 2 h 还属退潮时期，污染物在潮流的作用下向东南方移动。关闸的第 3 小时，污染物不再继续向东南方移动。由潮流特征，可知落潮 3 h 后，水流流速较大，污染物扩散速度也随之变大。所以相同位置对应时刻污染物的浓度较工况 1 小，且相应污染物影响范围比工况 1 小。当闭闸不排污时，由于潮流和水体自净功能的共同作用，之前 3 h 所排污染物浓度和影响范围逐渐变小。

3）工况 3 污染物变化幅度为 5.0～12.0 mg/L。开闸 3 h 时浓度和影响范围最大，9.0 mg/L 的污染带面积为 3.8 km2。此时污染物排放的时间属于落急至落憩点之间这段时间。污染物在横、纵向的扩散都加大，且浓度值相对较大，但关闸后紧接着涨潮，流向相反，污染物在潮流作用下向西北方移动。先前排放扩散的污染物回荡、叠加而使污染物浓度变大，且污染物向河口聚积，不易扩散，污染物浓度和影响范围缩小的速度慢，因而对污染排放不利。

4）比较各个工况的最大污染带面积，工况 2 中最大污染带面积小，4.0 mg/L 的最大污染带面积为 2.1 km2；工况 1 中 4.0 mg/L 的最大污染带面积为3.0 km2，比工况 2 情况差，但还能为环境所接受，没超标；工况 3 中 9.0 mg/L 的污染带面积为 3.8 km2，超标。且由 3 种工况比较得出，在工况 2 中，闭闸后污染物浓度和影响范围缩小最快，因此工况 2 对环境影响最小，在实际中，可按此种情况开闸排污。这为有效控制射阳河口近岸海域污染、科学制定污染物总量控制规划及射阳闸调度运行具有重要的应用价值。

4 结语

本文采用基于有限体积法及黎曼近似解的二维水流水质模型，模拟射阳河口的水动力场和污染物浓度场分布。将潮位、流速和流向的计算值与实测值进行比较，验证结果表明所建立射阳河口二维水量水质耦合模型较好地重现了射阳河口区域的潮位及潮流变化过程，说明模型中的参数取值合理，计算结果可信，为进一步研究该海域提供一定技术支持。

[1]Zhao D H, Shen H W, Tabios Ⅲ GO, et al. Finite-Volume Two-Dimensional Unsteady-Flow Model for River Basins[J].Hydr. Engrg, ASCE, 1994, 120 (7): 863-883.

[2]赵棣华，徐葆华. 平面二维水流—水质有限体积法及黎曼近似解模型[J]. 水科学进展，2000, 11(4): 368-373.

[3]Zhao D H, Shen H W, Lai J S, et al. Approximate Riemann Solvers in FVM for 2D Hydraulic Shock Waves Modeling[J].Hydr. Engrg, ASCE, 1996,122 (12): 693-702.