刘公岛岸滩整治方案数值模拟

， ，

(中海石油环保服务有限公司，天津 300452)

刘公岛位于山东半岛东端、黄海之滨的威海湾内，面积3.15 km2，与威海市区相距3.9 km(轮渡距离)，岛西、岛南与威海市区隔海相望，是国家级海洋生态特别保护区和国家级海洋公园。然而随着日久天长的风化、海蚀、地表径流切割、主体滑坡、风暴潮等自然灾害以及近年来威海湾沿岸的开发，海岛岛体面积逐年减小，生态环境日益恶化。南部沿岸沙滩严重流失，沙滩较5年前整体向海方向后退了1.5 m左右，地表裸露，遍布卵石、小石块，已经不适宜游人亲海游玩。

文中利用MIKE21软件建立二维潮流模型，针对刘公岛南部东村沿岸的沙滩开展了海洋水动力环境和冲淤数模研究，给出沙滩保护工程设计方案，为刘公岛南部海域岸滩整治修复工作提供技术支持。

1 MIKE21模型简介

采用丹麦水力学研究所研制的平面二维数值模型MIKE21FM来进行预测与分析。该模型采用非结构三角网格剖分计算域，三角网格能较好地拟合陆边界，网格设计灵活且可随意控制网格疏密，该软件具有算法可靠、计算稳定、界面友好、前后处理功能强大等优点，已在全球70多个国家得到应用，有上百例成功算例，计算结果可靠，为国际所公认。MIKE21FM采用标准Galerkin有限元法进行水平空间离散，在时间上，采用显式迎风差分格式离散动量方程与输运方程[1]。

1.1 模型控制方程

1)质量守恒方程。

2)动量方程。

式中：ζ——水位；

h——静水深；

H——总水深，H=h+ζ；

u、v——x、y方向垂向平均流速；

g——重力加速度；

f——科氏力参数(f=2ωsinφ，φ为计算海域所处地理纬度)；

εx、εy——x、y方向水平涡动粘滞系数。

1.2 定解条件

边界条件：固定边界取法向流速为零，即V·n=0；在潮滩区采用动边界处理；水边界采用预报潮位控制。

式中：A0——平均海面;

F0i、(v0+u)i——天文要素;

Hi、gi——某分潮调和常数，即振幅与迟角[2-5]。

2 潮流数值模拟

2.1 模型计算范围

2.1.1 计算域设置

建立海域数学模型计算域范围见图1，即为图中A、B、C三点以及岸线围成的海域。模拟采用非结构三角网格，将刘公岛南部东村沿岸附近海域网格进行局部加密(D、E、F三点与岸线围成的海域)，数值模拟局部域见图2，整个模拟区域最小空间步长约为25 m。用动边界的方法对干、湿网格进行处理。

图1 大海域计算域及网格设置

图2 局部海域网格设置

2.1.2 模型水边界输入

开边界：本次模拟的开边界水位由黄渤海潮流模型提供，其开边界潮位由下式输入计算。

式中：fi、σi——第i个分潮(这里取4个分潮);

M2、S2、O1和K1——交点因子和角速度；

Hi和Gi——调和常数，分别为分潮的振幅和迟角；

Voi+Vi——分潮的幅角。

闭边界：以大海域和工程周边岸线作为闭边界。

2.2 潮流数值模型及验证

潮流场数值模型采用计算区域的潮位、潮流流速和流向的实测资料对模型进行验证，其中潮位验证采用2010年9月的实测资料，潮流验证分别采用2010～2011年两次的实测资料，共设置了3个验潮点。潮位、潮流验证点的位置图见图3。潮位验证曲线见图4，潮流验证曲线见图5～图6。

图3 潮流潮位验证点位置

图4 1号点大潮潮位验证

验证结果表明，数值模拟得到的潮位和潮流与实测潮位和潮流资料基本吻合，能够较好地反映关注海域潮流状况。

图5 2号点大潮流速、流向验证

图6 3号点大潮流速、流向验证

2.3 现状潮流计算结果分析

图7是研究海域大潮期间涨急时潮流场现状，研究海域内的涨潮流整体由E往W流，外海潮流在刘公岛北端分叉，一部分潮流在刘公岛西南处转向东南，另一部分潮水沿刘公岛东北部岸线由西北向东南流动，大潮涨潮流速在刘公岛北部的东西两侧偏大，在0.3～0.5 m/s之间;而在威海湾内部潮流偏小，流速普遍介于0.10～0.30 m/s之间。

图7 研究海域大潮期间涨急时流场图

图8是研究海域大潮期间落急时潮流场现状，研究海域内的落潮流整体由W往E流，刘公岛西南侧海域潮流由SE往NW流，流速普遍介于0.10～0.30 m/s之间。

图8 研究海域大潮期间落急时流场图

2.4 设计方案潮流场计算结果

根据数值模拟结果，可以看出涨落潮在研究海域流场主要是东西方向的顺岸流动，而且流速较强，这可能导致泥沙运移造成冲刷。根据这一数值计算结果，针对刘公岛南部研究海域沙滩冲刷情况，设计如下方案[6]。

考虑到刘公岛南侧岸线较长，长约1 000 m，2个丁坝形式的设计可能难以阻止泥沙冲刷，故如图9所示，设计3个丁坝工况，进行数值模拟测试阻沙效果。

图9 设计方案布置图

设计方案的大潮涨急流场与落急流场分别见图10、图11。可以看出，相比现状条件，由于研究海域3条丁坝的建设，使得沙滩海域流态有所变化，流向有所改变，近岸流速有所减小，尤其是在丁坝的顶端，流速加强，而在3条丁坝所围的内区流速甚小。可见设计方案对沙滩海域潮流运动影响显著，显著地减小了坝之间海域的流速。

图10 设计方案涨急流场图

图11 设计方案落急流场图

3 地形地貌冲淤数值模拟

3.1 泥沙运动控制方程

MIKE21FM采用标准Galerkin 有限元法进行水平空间离散，在时间上，采用显式迎风差分格式离散动量方程与输运方程[7]。

泥沙控制方程为

u，v——水深平均流速,m/s；

Dx，Dy——分散系数,m2/s；

h——水深,m；

S——沉积/侵蚀源汇项,g/(m3·s)；

QL——单位水平区域内点源排放量,m3/(s·m2)；

CL——点源排放浓度,g/m3。

3.2 地形地貌与冲淤数值模拟结果

3.2.1 现状地形地貌与冲淤数值模拟结果

现状情况下研究海域的地形地貌冲淤数值模拟结果见图12。威海湾和杨家湾整体呈微淤积状态，淤积厚度小于4 cm/年，刘公岛东南侧和西南侧附近海域由于地形改变导致海域面积变窄，流速较大，处于微冲刷状态，每年的冲刷厚度小于6 cm。研究海域小范围处于微冲刷，每年冲刷程度小于4 cm，其余附近海域处于微淤积状态，淤积厚度小于4 cm/年。

在研究海域，由于地理环境特征和海洋水文动力特征，决定了该海域目前存在微冲刷状态，年冲淤率为3～6 cm/年之间，沙滩外侧冲刷相对更为严重。

图12 现状年冲淤率数值模拟结果(m/年)

3.2.2 设计方案地形地貌与冲淤数值模拟结果

设计方案研究海域的地形地貌冲淤数值模拟结果见图13，除了局部区域计算结果相比现状情况有所差别之外，其他海域并没有显著变化。威海湾和杨家湾整体呈微淤积状态，淤积厚度小于4 cm/年，刘公岛东南侧和西南侧附近海域由于地形改变导致海域面积变窄，流速较大，处于微冲刷状态，冲刷厚度小于6 cm/年。研究海域小范围处于微冲刷，冲刷程度小于4 cm/年，其余附近海域处于微淤积状态，淤积厚度小于4 cm/年。

图13 设计方案年冲淤率数值模拟结果(m/年)

由于3条丁坝的建设，研究海域近岸沙滩冲刷的现状有很大改观，在挡水坝两侧出现了较为显著的淤积现象，冲刷现状几乎消失，仅在西侧两条丁坝之间出现了很小范围的微冲地带，且冲刷率约1 cm/年，大部分沙滩的淤积率在2 cm/年以上，建坝促淤效果显著。

4 结束语

根据数值模拟结果，设计方案减小了近岸潮流，有效减缓了丁坝之间的水动力能量，防冲刷效果也比较明显，仅有小段岸线略有冲刷，基本形成了促淤态势，能够有效防止沙滩冲刷、岸线侵蚀。希望进一步优化详细工程方案，加强海洋调查和监测，对该段沙滩岸线进行有效的整治修复工作。

[1] DHI.Mike 21&Mike 3 Flow Model FM Scientific Documentation[M].Denmark:2011.

[2] 安永宁，杨 鲲.MIKE21模型在海洋工程研究中的应用[J].海岸工程，2013(3)：1-10.

[3] 包伟斌.MIKE21水动力数学模型在桥梁工程中流场模拟中的应用[J].华东科技,2014(4):238-239.

[4] 郑志飞.MIKE21模型在码头工程流场模拟中的应用研究[J].海峡科学,2010(10):104-106.

[5] 许 婷.丹麦MIKE21模型概述及应用实例[J].水利科技与经济,2010(8):867-869.

[6] 倪会东.东辽河河道冲淤分析及整治[J].水利科技与经济,2014(3):113-115.

[7] 吴园园，娄安刚.海阳中心渔港建设对附近海域冲淤的数值模拟预测[J].海岸工程,2014(1):50-52.