平面二维泥沙输移模型及其应用

第 11 卷 2011 年
第8期 8月
中 国 水 运 Chi na W er Tr a ns por t at
V . 11 ol Augus t
N 8 o. 2011
平面二维泥沙输移模型及其应用
夏雪瑾，高程程
（上海市水务规划设计研究院，上海 200232） 摘 要：为研究海岸工程对周边水沙环境和河床冲淤的影响，运用 Mik e2 1 建立了平面二维水沙数值模型， 并将模型应用
工程实例。

通过实测资料验证流场和含沙量场，结果表明模型能较好地模拟含沙量场。

在此基础上，从工程实施后引起的 水流变化以及年冲淤变化强度这些指标来判断工程实施对海域的影响。

研究成果可为今后类似的工程提供一定的参考。

关键词：二维悬沙模型；数值模拟；含沙量；泥沙回淤 中图分类号：TV 14 文献标识码：A 文章编号：1006- 7973（2011）08- 0082- 03 为零；开边界条件（水边界条件） s (x , y, t ) = s (x, y, t ) ： 二、应用实例 应用该 模型，建立舟山群岛海域二 维泥沙数学模型，计 算分析工 程给周边海域泥沙带来变化 ，为工程的决策提供技 术依据。

1．模型范围 考虑计 算水边界取在基本上不受本 工程影响的海域，同 时兼顾到 水文条件等相关资料获取的 方便，选取模型（右图 1 ）开边界的北边界为 3 0 °19 ′00 ″N，东边界 1 22 °4 0 ′0 0″E ，南边界 2 9°2 6′0 0″N，西边界 1 21 °37 ′0 0 ″E ； 模型闭边界为自然岸线， 包括象山港全域， 最西在 12 1 °2 5′2 4 ″E。

模型范围约 1 20 km ×10 0 km ，采用三角形 网格， 共布有 46 43 3 个三角形单元， 域内最大水深达 1 00 m 以上，最小空间步长约为 35 m 。

在河口海岸地 区兴建工程时，特别是兴 建重要的航运或 围垦工程时， 往往要求评估工程实施对周 边水域水沙环境的 影响。

河口海 岸泥沙运动规律研究比较复 杂，而在规划设计 阶段，通常要 求迅速给出比较概括性的解 答。

因此，数学模 型研究泥沙运 动被越来越广泛地应用于解 决近岸工程实际问 题。

水流泥沙 数学模型可对岸滩的演变、 海床冲淤等提供长 期预报，为工程的规划设计提供科学依据。

目前在悬沙、 底沙输移以及河床演变研 究中，二维泥沙 数值模型应用 最为广泛。

国内外相继出现 了一批功能强大， 通用性好、 成熟的商业化综合数学模型， TRIM 2D 模型[1 ]、 如 美国的 Miss ips ipp i 大学水科学计算中心 CC HE 2D 模型 、 丹麦水利所的 M iKE2 1 模型等， 大多数情况下具有足够的精 度能满足工程要求，被广泛应用。

本文主要利用 Mike 2 1 建 立二维平面水 沙数值模型研究工程对大范 围水域泥沙运动的 影响，为今后类似工程的规划设计提供一些借鉴。

一、潮流泥沙数学模型 1．基本方程 基本方程包括二维浅水方程（此处略）和二维悬沙方程。

悬沙基本方程：
S t + u S x + v S y = 1 h x h Dx S x + 1 h y hD y S y + Fs h
[2]
图 1 模型范围示意图 2．边界条件 本计 算域潮位开 边界由东中 国海潮波数学 模型提供[3 ]。

悬沙模块 由于只有同步的实时实测含 沙量资料，整个边界取 所有实测含沙量的平均值 0.2 5k g/ m 3 作为悬沙的边界条件。

3．计算参数的选取 （1 ）糙率 n ：由于模型计算域较大，根据模型研究区域 实际床面情况通过试算作适当调整，n 取值在 0 .02 5~ 0.0 4 之间。

（2 ）沉速：工程区域悬移质平均粒径为 0 .02 9m m ，易 受 海水 影响 发生 絮凝而 加速 沉降 ，因 此沉 速计 算公式 为： 。

ω = ω FD （F 为絮凝因子；D 为衰减系数） 0 （3 ）泥沙临界淤积切应力： 经过反复调试，选取临界淤 积切应力为 0 .1 3N/ m 2。

（4 ）泥沙临界冲刷切 应力τ 经过反复调试，选取泥 e ： 沙临界冲刷切应力在 0.1 5- 0.2 5N/ m 2 之间。

其中：S 为悬沙浓度，D x、D y 分别为 x、y 方向上的泥 沙扩散系数；F s 为泥沙冲淤函数。

2．泥沙冲淤函数 模型中泥沙冲 淤函数采用的是切应力法 ，由床面临界淤 积切应力和临界冲刷切应力确定源汇项： ω cb (1 ττ ) s d
Fs = 0 E (1 ττ)n e
τ< τ d τ < τ< τ d e τ≥ e τ
式中： 为瞬时底床剪切应力， d 为临界淤积切应力， e τ τ τ 为临界冲刷切应力， 为床面泥沙冲刷系数， E 由率定计算确定。

3．定解条件
初始条件：整 个计算域内每一个节点的 水位和流速、流
向由流场模型计算结果提供， 悬沙浓度初始值在开始时取零。

边界条件：闭 边界（即陆地边界）取含 沙量的法向梯度 收稿日期：2 01 1- 0 6 - 14 作者简介：夏雪瑾，上海市水务规划设计研究院。

第8期 4．模型验证
夏雪瑾等：平面二维泥沙输移模型及其应用
83
域的泥沙 运动规律，因此本文中采用 河口挟沙力公式。

根据 实测资料相关性分析拟合出来的河口挟沙力公式为： 大潮：S = 64.845(
v2 gh )0.4388
2 v2 0 .80 62 ) ； 中潮：S = 87.191( v ) 0.85 97 ； 小潮： gh gh
模型采用了 2 0 07 年地形资料和工程区附近水域的大、 中、小潮水文 实测资料进行验证，潮位验 证取定海、岱山、 桃花岛、沈家门、镇海共 5 个潮位站的逐时实测潮位，流速 和含沙量验证取用 1 1 条水文测线的实测值。

（1 ）潮位和流速验证 根据潮流模型验证结果（由于篇幅所限，验证图略） ：潮 位、 流速流向的计算值与实测值吻合良好， 符合相关规程要求； 工程区测流断面的涨落潮量计算值与实测值较为一致，很好地 反映了工程区域进出潮量的实际情况； 模型能够较好地反映工 程区的水流特性，并为悬沙模型提供较为准确的水动力条件。

（2 ）含沙量验证 泥沙数学模型 的困难之处在于泥沙运动 本身的复杂性和 验证资料的匮 乏。

本文从两方面来对含沙 量进行验证：①对 有实测资料的 1 1 个测站的含沙量过程线进行局部验证；② 整体含沙量分布趋势要能反映工程区海域泥沙运动的趋势。

1 ）含沙量过程线验证
含沙量 （ kg/m 3） 0 4 . 0 .3 5 0 3 . 0 .2 5 0 2 . 0 .1 5 0 1 . 0 .0 5 时间 0 1 1 0 :0 0 1 2 0 :0 0 1 3 0 :0 0 1 4 0 :0 0 1 5 0 :0 0 1 6 0 :0 0 17 0 00 : 18 0 00 : 19 0 00 : 20 0 00 : 2 1 0 0 : 0
S = 2.34(
式中：s 为全潮平均挟沙力，v 为全潮平均流速，h 为全 潮平均水深。

比较挟 沙力公式和悬沙模型计算的 工程海域的全潮平均 含沙量场分布图（图 4 ，以中潮为例） ，可见挟沙力公式与悬 沙模型计 算的全潮平均含沙量场分布 相当相似，但公式计算 的结果量 级上偏大一些。

主要是因为 资料在测验期间有偏北 风，风浪 大于年平均情况，因此采用 此次实测数据来拟合挟 沙力公式 也隐含考虑了波浪影响因素 ，同时也表明波浪对当 地海域含 沙量实际分布有一定影响。

总的来说，群岛东南向 外部海域 的含沙量相当小，而靠近杭 州湾的海域含沙量则远 大于舟山 群岛内部水域及东南方向的 群岛外部水域，这种分 布情况与 实际情况是吻合的，各个海 域含沙量的量级也与实 测资料中含沙量的量级及分布相一致。

P2
计 算值 实 测值
图 2 p2 含沙 量验证曲线 结果（代表测点 p 2 验证见图 2）表明，模型泥沙参数选 取比较合理， 计算值基本可以反映含沙量 随潮变化过程，涨 落潮最大误差、平均误差基本控制在 50 %以内，符合有关规 范规程的要求，基本反映了该海域的悬沙输运特征。

2 ）泥沙运动趋势验证 由于本区缺乏 整个工程研究区域完整的 、长历时的泥沙 实测资料，而 且前人对该区研究成果也较 少，因此，本文中 从国内常用水 流挟沙能力公式（刘家驹公 式、窦国仁公式、 张瑞瑾公式和 河口挟沙力公式）中选用公 式来计算得出大中 小潮平均含沙 量场，然后与用模型模拟的 大中小潮平均含沙 量场作比较，以做校核。

因为缺乏波浪 资料，因此没有选用刘家 驹公式和窦国仁 公式，主要在 张瑞瑾公式和河口挟沙力公 式之间考虑。

根据 实测大中小潮 数据并结合这两种公式对流 速与含沙量的相关 性进行拟合分析，相关分析拟合线（以小潮为例）见图 3 。

0
（a）拟合公式
（b）悬沙模型
图 4 中潮全潮平均含沙量 综上分 析，模型计算的含沙量场不 仅与实测资料吻合较 好，而且 能较好反映工程区泥沙运动 整体趋势，该模型可用 来进一步研究工程区海域泥沙淤积问题。

4．工程对附近海域的影响 在工程区海域平面二维水沙模型验证较好的基础上，应用 模型对工程（右图 5 ）后水流、泥沙进行计算，并采用涨落流 速的变化、 海床冲淤变化等指标分析工程实施后对海域的影响。

0
- 1
-1
- 2
-2 实测
- 3
实测 相关线 -3
相关 性
图5
对附近海域水流流场的影响
- 4
- 5 -6 -5 - 4 -3 -2 - 1 0
-4 -1 0 -8 -6 -4 -2 0
（1 ）对附近海域水流流场的影响 根据计 算，由于工程后截断了登步 岛与桃花岛间的潮通 道，使得 工程区近旁水域流速减小， 主要的流速减弱区为工 程区西北 和东南水域，此两水域为涨 、落潮流的流影区。

而 工程两侧 进潮通道流速均有所增大， 其中朱家尖与登步岛间 的流速增大较多 ，最大增大 0 .1m / s ，这种状况将更有利于 航道水深 的维护。

但是，螺头水道的 流速略有减少，主要是 由于总进潮量的减少导致，需加以注意。

图 3 小潮流速与含沙量相关性分析图（张瑞瑾公式左，河口 挟沙力公式右） 从相关系数 R 可以看出，张瑞瑾公式拟合相关系数太低 在 0 .3 以下， 而河口挟沙力公式拟合相关系数都在 0 .6 以 上， 表明全潮平均含沙量 S 与全潮平均流速 v 以及全潮平均水深
h 之间存 在相关的函数关系，公式能较好的拟合出工程区海


84 表1
分析点 虾峙门锚地 福利门水道 虾峙门水道 条帚门水道 佛渡水道 崎头洋 螺头水道 工程前水深 39 .1 81 .7 91 .5 52 .3 17 .9 31 .8 102.7
中 国 水 运 工程后各水道淤积强度比较
变幅 0.41% 不淤 不淤 不淤 不淤 0.13% 0.02% 分析点 象山港外航道 马峙锚地 沈家门渔港 0.16 不淤 不淤 不淤 不淤 0.04 0.02
第 11 卷 （单位：m a） /
工程前水深 13.6 13.9 3.9 39.4 8.6 14.1 49.6 淤积强度 不淤 不淤 不淤 0.13 0.20 0.16 0.27 变幅 不淤 不淤 不淤 0.33% 2.28% 1.14% 0.54%
淤积强度
工程区附近
（2 ）对附近海域冲淤的影响 在二维水沙数 模的基础上，采用半经验 半理论的回淤强 度公式来预测 工程给周边海域海床带来的 变化。

泥沙回淤强 度计算公式[4]采用： P =
αω * n TS S [1 ( *2 )] γ S*1 d
是通过含 沙量过程线和整体含沙量分 布趋势两方面的验证， 表明此模 型能基本正确反映计算区域 的含沙量场分布特征， 可利用此 模型进行工程对海域影响的 分析计算研究，从而为 研究类似工程建设中的海岸动力环境变化提供了有效手段。

需要说 明的是，由于模型以及实测 资料所限，模型未考 虑波浪的 作用，此外鉴于本文研究的 海域是无径流的浓度较 低的海域 ，在长江口这些泥沙运动机 理更加复杂的重要河口 水域的水沙模拟，该模型是否适用，还需要今后进一步研究。

参考文献 [1] Casulli V. Semi- implicit finite difference methods for two dimensional shallow water equations. J of Computational
Physic，1990，86：56～74
n 是一 年中潮数；α是沉降机率；P 是年回淤强度；S *1 和 是年平均全潮平均含沙量。

计算结果（表 1） 表明，工程的实施主要对工程附近海 域的回淤强度 影响最大，总的趋势是主要 在工程所在的清慈
门水道两端附 近产生较大的回淤强度，离 工程区越远，回淤
式中： ω为泥沙沉速；γ 是泥沙干容重；T 为潮周期； d
S *2 为工程前后对应于不同流速和水深的全潮平均含沙量；S *
强度越小，工 程对周边航道、锚地和码头 造成的冲淤影响较 小，年平均回淤强度基本上都在 0 .05 m / a 以下。

三、结语 本文运用 Mike 2 1 建立了适用于河口海岸区域的平面二 维水沙数学模 型，并将模型应用于工程实 例，分析工程建设 对水动力场的 影响以及泥沙回淤情况。

实 例应用表明模型计 算方法可靠， 验证基本反映工程区海域水 沙运动规律，特别 （上接第 8 1 页）
[2] Wang S S Y. C C ED2D Manual. C enter for C omputational Hydro science and Engineering .USA： The University of M ississippi， 1998 [3] 林珲，闾国年， 宋志尧.东 中国海潮波系统 与海岸演变模 拟研究[M].北京：科学出版社，2000 [4] 王义刚，林祥等，河口边滩围垦后淤积计算方法研究[J]， 海洋工程，2000 （3）
通过检验，河网多边形生成的正确率 1 00 %。

河网多边形的生成使得平原区产流分配的计算方法得到了 改进，优于传统的通过求算陆域宽平均分配入河道的方法。

但 因为河网多边形的产流分配只考虑河网多边形中网格与河道间 的距离和河道过水能力两个因素外，并没有兼顾地区的高程， 因此通过河网多边形进行产流的分配依然会存在一些不合理情 况，比如会将低洼处的产水分配入高程较高处的河道。

随着近 年来地理信息系统技术的日益成熟，将数字高程信息带入河网 多边形，从而来纠正产流分配中出现的上述不合理的状况是非 常有必要的，这也将是本课题下一步继续研究的工作方向。

参考文献 [1] 程文辉，王船海，朱琰. 太湖流域模型[M]. 南京：河海大 学出版社，2006：133~136.
图 5 太湖流域河网多边形的生成 五、小结 在平原区产汇流的背景下，为模拟平原区的坡面汇流，从 而提出河网多边形的概念。

河网多边形的构建，解决了平原区 产流河网分配以及下垫面信息的数值化的问题。

河网多边形的 生成是平原区坡面汇流计算中一项不可缺少的环节。

利用本文 所提出基于空 间方位角的最短路径搜索 的河网多边形生成算 法，经过河网数据的预处理、两类河网多边形生成及最后的后 处理， 最终以实例—太湖流域在 Mic r oso ft Vis u a l C+ + 2 00 5
的平台下实现了河网多边形的自动生成。

算法执行效率良好，
[2] 江斌，黄波，陆锋. GIS 环境下的空间分析和地学视觉化 [M]. 北京：高等教育出版社，2002. [3] 黄杏元，马劲松，汤勤. 地理信息系统概论[M]. 北京：高 等教育出版社，2001. [4] 周晓峰， 王志坚. “数字流域”剖析[J]. 计算机工程与应用 [J]. 2003， （3） ，104- 106. [5] 张欧阳，张红武. 数字流域及其在流域综合管理中的应用 [J]. 地理科学进展，2002， （1） ：66- 72. [6] Al G ore，T he Digital Earth ：U nderstanding our planet in the 21st Century.















。