舟山岱山岛燕窝山码头工程潮流泥沙数值模拟

李玉婷，宋志尧，李瑞杰，彭国强，陈 鹏

(1.南京师范大学 虚拟地理环境教育部重点实验室，南京 210023； 2.江苏省地理信息资源开发与利用协同创新中心，南京 210023；3.大规模复杂系统数值模拟江苏省重点实验室，南京 210023；4.河海大学 海岸灾害及防护教育部重点实验室，南京 210098； 5.河海大学 环境海洋实验室，南京 210098)

图1 燕窝山码头工程所在位置示意图Fig.1 Topographic map of simulated domain in Yanwoshan

近年来，为了开发利用深水岸线资源，修建码头和防护工程较为普遍，尤其在受风浪作用明显的海岸带，防护工程需要突出于海岸线一定距离才能满足防浪挡沙的要求，保证船舶靠泊条件[1]。在海区建设码头工程与防波堤，需要首先进行泥沙冲淤问题的研究。我国在这方面的研究成果很多。例如，李瑞杰等[2]通过建立潮流泥沙运动二维数学模型模拟河口的悬沙运动冲淤变化，为河口整治工程提供科学依据；宋志尧等[3]通过建立海岸河口二维水沙数学模型研究了长江口悬沙的分布特征和各主要入海通道的分沙比变化，为海岸河口相关工程的建设服务；柯杰等[4]在前人基础上通过建立多重嵌套的二维潮流泥沙数学模型，对枸杞岛拟建防波堤工程所引起的潮流场、含沙量场以及海床冲淤变化进行了分析和预测；曹成林等[5]以烟台八角海域为例，分别采用数值模拟计算与实测值对比两种方式分析评价了防波堤工程； 解鸣晓等[6]建立波流共同作用下的泥沙数学模型，研究了连云港口门防波堤建设对年平均含沙量场及进港航道回淤的影响。

岱山岛位于舟山群岛中部，在北纬30°13′～30°21′，东经121°3′～122°13′之间，北接长江口，西临杭州湾，东濒东海，南邻舟山本岛，是舟山市第二大岛，面积119.3 km2。岱山岛位于长江口外国际航线范围内，地理位置优越, 具有通江达海的区位优势。为满足舟山北向的水路客运需求，拟在岱山岛北侧建设燕窝山码头工程(图1)，包括码头和防波堤。根据燕窝山海港工程的规划，拟建码头为高桩式码头，栈桥为透空过水式，码头与陆域通过透空栈桥相连，具体方案详见本文第四部分。本文采用数值模拟计算分析码头工程实施后的流场以及泥沙淤积情况, 对现有防波堤平面布置方案对水流和泥沙输运的影响进行分析，选出合理解决方案。

1 工程海域概况

舟山市岱山县岱山岛属于北亚热带南缘季风海洋型气候。岛上存在沙滩、海涂、低山等多种地貌类型，受冷暖空气交替影响，灾害性天气频繁。

1.1 气象条件

此区域多年平均降水量为873.6 mm，6月份月降雨量最多，达121.9 mm，最大年降水量为128.6 mm(1996年)，全年大于25 mm的降水天数为6.6 d。多年平均风速为6.6 m/s，最大风速为38.6m/s。台风是影响工程区域的主要灾害性天气系统。

1.2 潮汐潮流条件

工程海域潮汐以M2分潮为主。工程海域的潮汐性质为不规则半日潮，潮汐不对称现象明显，三个临时潮位站均有浅水分潮存在。工程海域潮流运动形式为往复流。实测最大涨潮流流速为2.28 m/s，流向为216°；最大落潮流流速为2.52 m/s，流向为77°；垂向平均的最大涨潮流流速为2.17 m/s，流向为237°；垂向平均的最大落潮流流速为2.34 m/s，流向为78°。

1.3 泥沙条件

工程海域底质以细颗粒的粉砂和粘土为主，其中粉砂约占69%、粘土约占25%，此外还有少量的砂约占6%。底质中值粒径范围在3.72～6.51 μm之间，平均为5.17 μm。底质中值粒径的空间分布有东面粗、西面细的特征。实测最大含沙量为1.874 kg/m3，最小含沙量为0.171 kg/m3。垂向平均含沙量最大值为1.495 kg/m3，最小值为0.342 kg/m3(图2)。

图2 岱山岛北部海洋工程水文泥沙工程站位Fig.2 Engineering station of northern Daishan island

2 水沙动力学方程

采用建立大、小两套模型进行计算，大模型为岱山岛海域，小模型为码头工程局部海域。大模型为小模型提供水动力模型所需要的边界条件，以保证码头工程局部流场计算符合岱山岛海域潮流场的整体物理特征。

模型采用垂向平均二维浅水方程，采用 ADI法(Alternating Direction Implicit)和DS法(Double Sweep)格式离散，采用非结构三角形网格和有限体积方法进行数值离散和求解。高桩码头支撑桩对工程海域流场阻流作用的模拟采用等效阻力的方法。

动力学方程包括：

连续性方程

(1)

式中：d=h+η为总水深，η为水位，h为静水深；t为时间；u、v为垂线平均流速分别在x、y方向上的分量。

水深平均x、y方向的动量方程分别为

(2)

(3)

式中：g为重力加速度；f=2ωesinφ为科氏参数，ωe为地球自转角速度，φ为当地纬度；εxx、εxy、εyx和εyy为不同方向上涡粘系数；CZ为谢才系数。

二维悬沙输移扩散方程为

(4)

式中：C为垂向平均含沙量；εx、εy为扩散系数；Fc为水沙界面通量，由下式给出

(5)

式中：α为沉积系数；ω为泥沙沉速；M为冲刷系数；τb为底部切应力；τe为临界冲刷切应力；τd为临界淤积切应力。

地形变化方程为

(6)

式中：γ0为泥沙干密度，Zb为床面高程。

不考虑海面风作用,陆边界

Vn(x,y,z,t)=0(n为陆边法向)

(7)

开边界处给定水位,由东中国海潮波模型[7]计算得到。

悬沙闭边界条件采用法向零通量边界条件，悬沙开边界条件如下式

(8)

式中：Γ为水域开边界；C*(x,y,t)为已知悬沙浓度。

3 模型建立及水沙模拟

3.1 模型的构建

图3 大、小范围数学模型计算区域Fig.3 The large and small range of calculated area with mathematical model

模型计算区域如图3所示，计算区域北面以长江口以南、南汇嘴以东的东海海域为北开边界，南面以六横岛以南约10 km的东海海域为南开边界，东面以嵊山岛以东约40 km的东海海域为东开边界，西面以黄湾镇南杭州湾口为西开边界，其他边界均为固边界。

大、小模型均采用非结构三角形网格,其中大模型最小网格尺寸为50 m，节点数为15 468，网格数为28 987。小模型最小网格尺寸为5 m , 网格单元34 808个，网格节点17 289个。结果输出时间步长为1 800 s，计算时段为2017年3月21日00:00时～2017年4月7日00:00；经率定，舟山海域糙率取为0.02；本文利用考虑亚尺度网格效应的Smagororinsky(1963)公式计算水平涡粘系数，取为0.28；由东中国海潮波数学模型确定模型开边界的外海开边界条件。大、小模型的计算域及网格剖分见图4、图5所示。

3.2 模型的验证

图4 大模型计算区域网格剖分图 图5 小模型计算区域网格剖分图Fig.4 Mesh map of large model areaFig.5 Mesh map of small model area

模型计算的水位、流速、振幅以及过程线，均与原型验证良好，符合《海岸与河口潮流泥沙模拟技术规程》[8]的要求，限于篇幅，不再一一列出。关于含沙量，实测值依据舟山北向客运枢纽(燕窝山)工程水文泥沙测验分析报告(2017年)[9]，图6是6个测站的含沙量验证图。含沙量的验证较好地反映工程海域泥沙的输运特征。

6-a SW1测站6-b SW2测站6-c SW3测站

6-d SW4测站6-e SW5测站6-f SW6测站图6 各测站大潮含沙量验证图Fig.6 Sediment concentration verification

3.3 水动力模拟

工程区潮流数学模型计算流场如图7所示，以大潮为例，大潮涨潮时，潮流由东南方向进入计算区域，流经衢山岛和长涂山时发生分流。水流绕过燕窝岛后进入工程海域，近岸海域流速较小。大潮落潮时，来自杭州湾的水流绕过大鱼山向西北方向流去，部分水流进入工程海域，进一步绕过燕窝岛和鲞蓬山，穿过岱山、长涂山与衢山岛中间的海域与南侧落潮流汇合进入东海。

7-a 涨急7-b 落急 图7 小模型工程实施前大潮涨落急流场图Fig.7 Rapid flow field of tidal fluctuation before the implementation of small model engineering

4 工程应用计算

4.1 工程方案情况

根据舟山市交通规划设计院提供舟山北向客运枢纽工程工程可行性研究报告(2017年)[10]中的燕窝山陆岛交通码头工程总体布置图，拟建的燕窝山码头位于舟山市岱山县岱山岛北部，地理位置约北纬30°20′48″，东经121°10′23″。基本布置方案为两种如下：

方案一：防波堤为折线型布置，堤长1 200 m，防波堤西侧折线段长度100 m(方位角为N46°～N226°)，东侧折线段长度1 100 m(方位角为N76°～N256°)，拟建码头与拟建防波堤的东侧折线段以N76°～N256°的方位角平行布置。码头陆域位于防波堤南侧陆地，码头与陆域通过透空栈桥相连，陆域场地位于拟建码头东南侧(图8)。

方案二：防波堤轴线布置同方案一。码头陆域位于燕窝山西侧，码头与陆域通过透空栈桥相连，陆域场地位于拟建码头东侧(图9)。

图8 方案一平面布置示意图Fig.8 Layout of plan I

图9 方案二平面布置示意图Fig.9 Layout of plan II

具体建设方案如表1所示。

表1 各方案建设内容Tab.1 Construction content of each plan

4.2 方案实施后流场变化

4.2.1 工程实施后潮流运动特征

图10为方案一建设后工程海域潮流的大潮涨落急流场图，图11为方案二建设后工程海域潮流的大潮涨落急流场图。由图可知，大潮涨潮时，水流由燕窝山东侧流入工程海域，流经燕窝山后水流发生分离，在燕窝山西侧形成回流区。由于防波堤对水流的阻挡，水流由防波堤东侧口门处进入港池及近岸海域，水流流量与防波堤建设前相比大幅减少，同时在防波堤后侧港池水域形成回流区，流速均小于0.5 m/s。水流流至东垦山时再次发生分离，主流部分继续向西运动，小部分水流向南侧近岸海域运动形成大面积回流。大潮落潮时水流由工程海域西北侧流入防波堤附近海域，流经东垦山时主流继续向燕窝山运动，小部分水流流入近岸浅水海域及港池水域。由于陆地边界处水深较浅，落急时刻部分边界区域露出水面。防波堤的建设导致水流流经防波堤时发生分离，并在紧靠防波堤北侧边界区域形成回流区。由于工程海域地形较为复杂，导致防波堤南北两侧流速差异较大。防波堤北侧海域流速较大，西侧堤头以北2 km处涨急流速达1.75 m/s，落急流速达2.01 m/s。防波堤以南浅水海域流速基本小于0.5 m/s。中潮与小潮涨落潮流情况与大潮相似，但潮流动力较弱，并呈现落潮流流速略大于涨潮流流速的趋势。由于岱山北侧附近地形复杂，水深梯度大、岛屿礁石众多，无论涨潮流还是落潮流均会形成许多水流漩涡，工程海域附近潮流仍以往复流流态为主。

10-a 涨急10-b 落急 图10 方案一工程实施后大潮涨落急局部流场图Fig.10 The tide flood maximum map after the project I implementation

11-a 涨急11-b 落急 图11 方案二工程实施后大潮涨落急局部流场图Fig.11 The tide flood maximum map after the project II implementation

4.2.2 工程建设前后流速变化比较分析

工程前大潮涨潮时，防波堤南、北侧特征站位平均流速和最大流速均呈现由东向西递增的趋势。大潮落潮时，防波堤南、北侧平均流速和最大流速均呈现由西向东递减的趋势。图12为方案一建设后工程海域潮流的大潮涨落急流速等值线图，图13为方案二建设后工程海域潮流的大潮涨落急流速等值线图。根据图12～图13和方案一、二实施前后小范围潮流数学模型的数值模拟结果，以大潮为例，除了部分站位流速有所增大以外，其余流速均有所减小，最大流速变化趋势和平均流速变化趋势基本一致。

12-a 涨急12-b 落急13-a 涨急13-b 落急图12 方案一工程实施后大潮涨落急局部流速等值线图Fig.12 The velocity contour map after the plan I implementation图13 方案二工程实施后大潮涨落急局部流速等值线图Fig.13 The velocity contour map after the plan II implementation

4.3 方案实施后淤积变化

各主要方案年淤积强度如图14所示，其中“正”表示淤积增加或冲刷减弱，“负”表示淤积减弱或冲刷加强。不同方案工程实施后航道码头区域年最大淤积强度和年回淤量见表2。由图表可知，方案一工程实施后，防波堤附近海域泥沙冲淤强度有所变化，由于防波堤纵向轴线与等深线存在一定夹角，水流流至防波堤附近时流速减小，泥沙形成落淤。落潮时由于防波堤的建设，回流区流速减小，泥沙淤积强度较大；回流区两侧底床有轻微冲刷现象。防波堤东侧口门区域由于防波堤的束水作用，导致防波堤堤头区域流速增大，泥沙冲刷强度增大。防波堤工程实施后航道和港池水域由于开挖疏浚将导致泥沙形成落淤，淤积强度最大处位于航道弯曲段。方案二防波堤的布置方案与方案一一致，仅陆域场地设置在东侧燕窝山附近，对工程区域泥沙冲淤强度影响较小，其冲淤强度变化趋势与方案一类似。

由上述分析可以看出，方案一实施后航道和港池水域年淤积总量为8.65万m3，比选方案中，方案二防波堤布置形式引起的航道和港池水域回淤量为9.53万m3。

图14 不同方案工程实施前后淤积强度变化分布图Fig.14 Distribution of siltation intensity before and after implementation of different plans

表2 各方案工程实施后航道和港池水域最大淤积强度及年淤积总量Tab.2 The maximum siltation intensity and the total annual sedimentation volume of waterways and harbor basins after each plan implementation

4.4 方案实施后航道水流条件分析

对工程前和各方案工程实施后航道口门段沿航道各特征站位的涨潮和落潮期间横流流速进行统计分析，结果见表3和表4。防波堤工程实施前，落潮期间的横流起控制作用。方案一航道口门段横流流速沿航道向港池方向逐渐减小。在涨潮期间，方案二航道口门段横流流速沿航道向港池方向逐渐减小。工程实施后，各方案仍是落潮期间的最大横流大于涨潮期间的最大横流，落潮期间的横流起控制作用。

方案一实施后，在涨潮期间，航道口门段横流流速沿航道向港池方向先增大后减小，在C2处取最大值1.11 m/s，C5处取最小值0.24 m/s；在落潮期间，航道口门段横流流速沿航道向港池方向逐渐减小，大小为0.56～1.33 m/s；连续两个潮周期内，横流大于0.50 m/s的时长沿航道向港池方向逐渐减小，最大为17.50 h，最小为2.50 h。与方案实施前相比，方案实施后涨潮期间航道口门前半段(C1、C2、C3)横流流速增大，最大增量为0.21 m/s(C2)，后半段减小，最大减小量为0.57 m/s(C5)；涨潮期间各点均减小，最大减小量为0.48 m/s(C5)。

方案二实施后，航道口门段在涨落潮期间横流流速沿航道向港池方向变化趋势与方案一相似，涨潮期间在C3处取最大值1.11 m/s，C5处取最小值0.34 m/s；落潮期间在C1处取最大值1.40 m/s，C5处取最小值0.65 m/s；连续两个潮周期内，横流大于0.50 m/s的时长沿航道向港池方向逐渐减小，最大为17.50 h，最小为4.00 h。与方案一相比，方案二航道口门段横流流速除涨潮期间点C1和C2外均大于方案一，最大差值为0.34 m/s，且横流流速大于0.50 m/s的时长也均大于方案一。

表3 工程前航道特征站位横流流速Tab.3 Cross current velocity of channel characteristic stations before project implementation

表4 各方案实施后航道特征站位横流流速Tab.4 Cross current velocity of channel characteristic stations after project implementation

5 结论

综合泥沙数学模型计算结果的主要结论如下：

(1)工程海域附近潮流以往复流为主，燕窝山和东垦山的北海域流速较大，超过2.0 m/s，近岸海域水深变化梯度较大，5 m等深线以内水流流速较小，在0.6 m/s以内。

(2)各方案防波堤工程实施后，对防波堤附近海域潮位的影响很小。方案一实施后，流速变化最大的特征站位与工程前相比增大32.88%，各特征站位最大流速均在1.0 m/s以内。码头前沿特征站位最大横流流速为0.15 m/s，港池水域回流强度为0.37 m/s，航道特征站位最大横流流速为1.33 m/s；方案二实施后，码头前沿特征站位最大横流流速为0.12 m/s，港池水域回流强度为0.49 m/s，航道特征站位最大横流流速为1.40 m/s。

(3)计算结果表明，方案一实施后，航道和港池水域最大年淤积强度为0.71 m/a，年淤积总量为8.65万m3，方案二实施后, 航道和港池水域最大年淤积强度为0.69 m/a，年淤积总量为9.53万m3。

综上所述，从水流、泥沙计算结果看，燕窝山码头工程实施后并未对岱山北部海域的整体流态产生较大影响。但根据航道横流流速变化特征调整航道设计方案，应增大航道设计宽度，减小航道轴线与航道口门附近水流流向的夹角，保证船舶航行安全；从水深条件角度考虑，方案二码头位置水深较方案一好,挖泥浚深量较小。从淤积角度考虑，方案一年淤积总量略小于方案二。从文中模拟出的相关流速数据角度考虑，方案二仅码头前沿特征站位最大横流略小于方案一，而方案二港池水域回流强度和航道特征站位最大横流流速均大于方案一。从建设角度考虑，方案一利用燕窝山遮挡了东向风浪，故减少了防波堤长度，从而减少工程造价；陆域东侧为消浪区，码头泊稳条件好。而方案二，港区全由防波堤掩护，防波堤较长，工程造价较高；且码头前沿受反射波浪影响，不利于泊稳。故综合评价为方案一优于方案二。