第11卷第8期中国水运Vol.11No.82011年8月ChinaWaterTransportAugust2011

收稿日期：6

作者简介：夏雪瑾，上海市水务规划设计研究院。平面二维泥沙输移模型及其应用

夏雪瑾，高程程

（上海市水务规划设计研究院，上海200232）

摘要：为研究海岸工程对周边水沙环境和河床冲淤的影响，运用Mike21建立了平面二维水沙数值模型，并将模型应用

工程实例。通过实测资料验证流场和含沙量场，结果表明模型能较好地模拟含沙量场。在此基础上，从工程实施后引起的

水流变化以及年冲淤变化强度这些指标来判断工程实施对海域的影响。研究成果可为今后类似的工程提供一定的参考。

关键词：二维悬沙模型；数值模拟；含沙量；泥沙回淤

中图分类号：TV14文献标识码：A文章编号：1006-7973（2011）08-0082-03

在河口海岸地区兴建工程时，特别是兴建重要的航运或

围垦工程时，往往要求评估工程实施对周边水域水沙环境的

影响。河口海岸泥沙运动规律研究比较复杂，而在规划设计

阶段，通常要求迅速给出比较概括性的解答。因此，数学模

型研究泥沙运动被越来越广泛地应用于解决近岸工程实际问

题。水流泥沙数学模型可对岸滩的演变、海床冲淤等提供长

期预报，为工程的规划设计提供科学依据。

目前在悬沙、底沙输移以及河床演变研究中，二维泥沙

数值模型应用最为广泛。国内外相继出现了一批功能强大，

通用性好、成熟的商业化综合数学模型，如TRIM2D模型[1]、

美国的Missipsippi大学水科学计算中心CCHE2D模型[2]、

丹麦水利所的MiKE21模型等，大多数情况下具有足够的精

度能满足工程要求，被广泛应用。本文主要利用Mike21建

立二维平面水沙数值模型研究工程对大范围水域泥沙运动的

影响，为今后类似工程的规划设计提供一些借鉴。

一、潮流泥沙数学模型

1．基本方程

基本方程包括二维浅水方程（此处略）和二维悬沙方程。

悬沙基本方程：

hFyShDyhxShDxhySvxSutSsyx++=++11

其中：S为悬沙浓度，Dx、Dy分别为x、y方向上的泥

沙扩散系数；Fs为泥沙冲淤函数。

2．泥沙冲淤函数

模型中泥沙冲淤函数采用的是切应力法，由床面临界淤

积切应力和临界冲刷切应力确定源汇项：()

()≥<<<=

eneedddbs

sEcF

τττττττττττω

101

式中：τ为瞬时底床剪切应力，dτ为临界淤积切应力，eτ

为临界冲刷切应力，E为床面泥沙冲刷系数，由率定计算确定。

3．定解条件

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

边界条件：闭边界（即陆地边界）取含沙量的法向梯度为零；开边界条件（水边界条件）：()()tyxstyxs,,,,=二、应用实例

应用该模型，建立舟山群岛海域二维泥沙数学模型，计

算分析工程给周边海域泥沙带来变化，为工程的决策提供技

术依据。

1．模型范围

考虑计算水边界取在基本上不受本工程影响的海域，同

时兼顾到水文条件等相关资料获取的方便，选取模型（右图

1）开边界的北边界为30°19′00″N，东边界122°40

′00″E，南边界29°26′00″N，西边界121°37′00

″E；模型闭边界为自然岸线，包括象山港全域，最西在121

°25′24″E。模型范围约120km×100km，采用三角形

网格，共布有46433个三角形单元，域内最大水深达100m

以上，最小空间步长约为35m。

图1模型范围示意图

2．边界条件

本计算域潮位开边界由东中国海潮波数学模型提供[3]。

悬沙模块由于只有同步的实时实测含沙量资料，整个边界取

所有实测含沙量的平均值0.25kg/m3作为悬沙的边界条件。

3．计算参数的选取

（1）糙率n：由于模型计算域较大，根据模型研究区域

实际床面情况通过试算作适当调整，n取值在0.025~0.04

之间。

（2）沉速：工程区域悬移质平均粒径为0.029mm，易

受海水影响发生絮凝而加速沉降，因此沉速计算公式为：

FD0ωω=（F为絮凝因子；D为衰减系数）。

（3）泥沙临界淤积切应力：经过反复调试，选取临界淤

积切应力为0.13N/m2。

（4）泥沙临界冲刷切应力eτ：经过反复调试，选取泥

沙临界冲刷切应力在0.15-0.25N/m2之间。

2011-0-14第8期夏雪瑾等：平面二维泥沙输移模型及其应用83

4．模型验证

模型采用了2007年地形资料和工程区附近水域的大、

中、小潮水文实测资料进行验证，潮位验证取定海、岱山、

桃花岛、沈家门、镇海共5个潮位站的逐时实测潮位，流速

和含沙量验证取用11条水文测线的实测值。

（1）潮位和流速验证

根据潮流模型验证结果（由于篇幅所限，验证图略）：潮

位、流速流向的计算值与实测值吻合良好，符合相关规程要求；

工程区测流断面的涨落潮量计算值与实测值较为一致，很好地

反映了工程区域进出潮量的实际情况；模型能够较好地反映工

程区的水流特性，并为悬沙模型提供较为准确的水动力条件。

（2）含沙量验证

泥沙数学模型的困难之处在于泥沙运动本身的复杂性和

验证资料的匮乏。本文从两方面来对含沙量进行验证：①对

有实测资料的11个测站的含沙量过程线进行局部验证；②

整体含沙量分布趋势要能反映工程区海域泥沙运动的趋势。

1）含沙量过程线验证

P2

00.050.10.150.20.250.30.350.4

110:00120:00130:00140:00150:00160:00170:00180:00190:00200:00210:00时间含沙量（kg/m3）计算值实测值

图2p2含沙量验证曲线

结果（代表测点p2验证见图2）表明，模型泥沙参数选

取比较合理，计算值基本可以反映含沙量随潮变化过程，涨

落潮最大误差、平均误差基本控制在50%以内，符合有关规

范规程的要求，基本反映了该海域的悬沙输运特征。

2）泥沙运动趋势验证

由于本区缺乏整个工程研究区域完整的、长历时的泥沙

实测资料，而且前人对该区研究成果也较少，因此，本文中

从国内常用水流挟沙能力公式（刘家驹公式、窦国仁公式、

张瑞瑾公式和河口挟沙力公式）中选用公式来计算得出大中

小潮平均含沙量场，然后与用模型模拟的大中小潮平均含沙

量场作比较，以做校核。

因为缺乏波浪资料，因此没有选用刘家驹公式和窦国仁

公式，主要在张瑞瑾公式和河口挟沙力公式之间考虑。根据

实测大中小潮数据并结合这两种公式对流速与含沙量的相关

性进行拟合分析，相关分析拟合线（以小潮为例）见图3。

-5-4-3-2-10

-6-5-4-3-2-10实测相关

性

-4-3-2-10

-10-8-6-4-20实测相关

线

图3小潮流速与含沙量相关性分析图（张瑞瑾公式左，河口

挟沙力公式右）

从相关系数R可以看出，张瑞瑾公式拟合相关系数太低

在0.3以下，而河口挟沙力公式拟合相关系数都在0.6以上，

表明全潮平均含沙量S与全潮平均流速以及全潮平均水深

之间存在相关的函数关系，公式能较好的拟合出工程区海域的泥沙运动规律，因此本文中采用河口挟沙力公式。根据

实测资料相关性分析拟合出来的河口挟沙力公式为：

大潮：20.806264.845()vSgh=；中潮：20.859787.191()vSgh=；小潮：

20.43882.34()vSgh=

式中：s为全潮平均挟沙力，v为全潮平均流速，h为全

潮平均水深。

比较挟沙力公式和悬沙模型计算的工程海域的全潮平均

含沙量场分布图（图4，以中潮为例），可见挟沙力公式与悬

沙模型计算的全潮平均含沙量场分布相当相似，但公式计算

的结果量级上偏大一些。主要是因为资料在测验期间有偏北

风，风浪大于年平均情况，因此采用此次实测数据来拟合挟

沙力公式也隐含考虑了波浪影响因素，同时也表明波浪对当

地海域含沙量实际分布有一定影响。总的来说，群岛东南向

外部海域的含沙量相当小，而靠近杭州湾的海域含沙量则远

大于舟山群岛内部水域及东南方向的群岛外部水域，这种分

布情况与实际情况是吻合的，各个海域含沙量的量级也与实测资料中含沙量的量级及分布相一致。

（a）拟合公式（b）悬沙模型

图4中潮全潮平均含沙量

综上分析，模型计算的含沙量场不仅与实测资料吻合较

好，而且能较好反映工程区泥沙运动整体趋势，该模型可用

来进一步研究工程区海域泥沙淤积问题。

4．工程对附近海域的影响

在工程区海域平面二维水沙模型验证较好的基础上，应用

模型对工程（右图5）后水流、泥沙进行计算，并采用涨落流

速的变化、海床冲淤变化等指标分析工程实施后对海域的影响。

图5对附近海域水流流场的影响

（1）对附近海域水流流场的影响

根据计算，由于工程后截断了登步岛与桃花岛间的潮通

道，使得工程区近旁水域流速减小，主要的流速减弱区为工

程区西北和东南水域，此两水域为涨、落潮流的流影区。而

工程两侧进潮通道流速均有所增大，其中朱家尖与登步岛间

的流速增大较多，最大增大0.1m/s，这种状况将更有利于

航道水深的维护。但是，螺头水道的流速略有减少，主要是

由于总进潮量的减少导致，需加以注意。v

h84中国水运第11卷

表1工程后各水道淤积强度比较（单位：m/a）

分析点工程前水深淤积强度变幅分析点工程前水深淤积强度变幅虾峙门锚地39.10.160.41%象山港外航道13.6不淤不淤福利门水道81.7不淤不淤马峙锚地13.9不淤不淤虾峙门水道91.5不淤不淤沈家门渔港3.9不淤不淤条帚门水道52.3不淤不淤39.40.130.33%佛渡水道17.9不淤不淤8.60.202.28%崎头洋31.80.040.13%14.10.161.14%螺头水道102.70.020.02%工程区附近49.60.270.54%（2）对附近海域冲淤的影响

在二维水沙数模的基础上，采用半经验半理论的回淤强

度公式来预测工程给周边海域海床带来的变化。泥沙回淤强

度计算公式[4]采用：)](1[1*2**SSTSnPd=γωα

式中：ω为泥沙沉速；dγ是泥沙干容重；T为潮周期；

n是一年中潮数；α是沉降机率；P是年回淤强度；S*1和

S*2为工程前后对应于不同流速和水深的全潮平均含沙量；S*是年平均全潮平均含沙量。

计算结果（表1）表明，工程的实施主要对工程附近海

域的回淤强度影响最大，总的趋势是主要在工程所在的清慈

门水道两端附近产生较大的回淤强度，离工程区越远，回淤强度越小，工程对周边航道、锚地和码头造成的冲淤影响较

小，年平均回淤强度基本上都在0.05m/a以下。

三、结语

本文运用Mike21建立了适用于河口海岸区域的平面二

维水沙数学模型，并将模型应用于工程实例，分析工程建设

对水动力场的影响以及泥沙回淤情况。实例应用表明模型计

算方法可靠，验证基本反映工程区海域水沙运动规律，特别是通过含沙量过程线和整体含沙量分布趋势两方面的验证，

表明此模型能基本正确反映计算区域的含沙量场分布特征，

可利用此模型进行工程对海域影响的分析计算研究，从而为

研究类似工程建设中的海岸动力环境变化提供了有效手段。

需要说明的是，由于模型以及实测资料所限，模型未考

虑波浪的作用，此外鉴于本文研究的海域是无径流的浓度较

低的海域，在长江口这些泥沙运动机理更加复杂的重要河口

水域的水沙模拟，该模型是否适用，还需要今后进一步研究。

参考文献

[1]CasulliV.Semi-implicitfinitedifferencemethodsfortwo

dimensionalshallowwaterequations.JofComputational

Physic，1990，86：56～74

[2]WangSSY.CCED2DManual.CenterforComputational

HydroscienceandEngineering.USA：TheUniversityof

Mississippi，1998

[3]林珲，闾国年，宋志尧.东中国海潮波系统与海岸演变模

拟研究[M].北京：科学出版社，2000

[4]王义刚，林祥等，河口边滩围垦后淤积计算方法研究[J]，

海洋工程，2000（3）

（上接第81页）

图5太湖流域河网多边形的生成

五、小结

在平原区产汇流的背景下，为模拟平原区的坡面汇流，从

而提出河网多边形的概念。河网多边形的构建，解决了平原区

产流河网分配以及下垫面信息的数值化的问题。河网多边形的

生成是平原区坡面汇流计算中一项不可缺少的环节。利用本文

所提出基于空间方位角的最短路径搜索的河网多边形生成算

法，经过河网数据的预处理、两类河网多边形生成及最后的后

处理，最终以实例—太湖流域在MfV++5

的平台下实现了河网多边形的自动生成。算法执行效率良好，通过检验，河网多边形生成的正确率100%。

河网多边形的生成使得平原区产流分配的计算方法得到了

改进，优于传统的通过求算陆域宽平均分配入河道的方法。但

因为河网多边形的产流分配只考虑河网多边形中网格与河道间

的距离和河道过水能力两个因素外，并没有兼顾地区的高程，

因此通过河网多边形进行产流的分配依然会存在一些不合理情

况，比如会将低洼处的产水分配入高程较高处的河道。随着近

年来地理信息系统技术的日益成熟，将数字高程信息带入河网

多边形，从而来纠正产流分配中出现的上述不合理的状况是非

常有必要的，这也将是本课题下一步继续研究的工作方向。

参考文献

[1]程文辉，王船海，朱琰.太湖流域模型[M].南京：河海大

学出版社，2006：133~136.

[2]江斌，黄波，陆锋.GIS环境下的空间分析和地学视觉化

[M].北京：高等教育出版社，2002.

[3]黄杏元，马劲松，汤勤.地理信息系统概论[M].北京：高

等教育出版社，2001.

[4]周晓峰，王志坚.“数字流域”剖析[J].计算机工程与应用

[J].2003，（3），104-106.

[5]张欧阳，张红武.数字流域及其在流域综合管理中的应用

[J].地理科学进展，2002，（1）：66-72.

[6]G，TD：U

yicrosotisualC200AloreheigitalEarthnderstandingourplanetinthe

21stCentur.