波流联合作用下广东阳西电厂海域二维泥沙数值模拟研究

谭 超，黄本胜，邱 静，黄广灵，刘 达

(1.广东省水利水电科学研究院，广州 510630；2.河口水利技术国家地方联合工程实验室，广州 510630；3. 广东省水动力学应用研究重点实验室，广州 510630)

0 引 言

在河口海岸地区, 引起泥沙运动的主要动力因素是潮流和波浪, “波浪掀沙、潮流输沙”被公认为是泥沙运动机制[1,2]。近岸水动力环境对海岸带及港工建筑物的影响十分重要，而潮流、波浪和泥沙等相互作用又使得海岸带附近的水动力环境非常复杂。因此，在复杂动力条件下，海岸工程建设对沿岸泥沙输运的影响过程研究备受国内外学者关注，也取得了系列研究成果[3-9]。在对其诸多研究方法中，以数学模型最为方便经济，并成为研究大范围水动力条件变化的主要手段之一[10]。 电厂的建设改变了岸线形态，导致水动力条件和泥沙运动的变化。同时，华南海岸经常遭受不同程度台风的侵袭，台风期间产生的大浪和风暴潮极易掀起海床泥沙，使含沙量剧增，泥沙被水流搬运至近岸处发生骤淤，对近岸港口、航道及取、排水口的安全等造成较大影响[11,12]。

广东阳西火力发电厂位于广东省阳江市海陵湾的外湾-青海湾的西端海岸地段。青海湾呈半圆形，为岬角控制的弧形海湾，其东侧为海陵岛闸坡岬角，西端为白虎岭岬角。2个岬角在地貌上均属于侵蚀剥蚀丘陵，由花岗岩基岸及风化壳组成。岬角附近海域分布有礁石和小岛屿，较大的有双山岛。湾口大体向南，东西宽约14 km，湾顶为沙质海滩，其东北端为与海陵湾内湾相连潮汐水道。清海湾湾口水深约10 m，海湾东部大部水深为2～5 m，西部则一般为5～10 m。工程附近海域双山岛至马村-白虎岭沿岸水深多变，一般为2～5 m，最深处为白虎岭南侧岬角附近，其10 m等深线距岸仅400 m左右， 研究区域位置见图1。

图1 研究区域位置

本文拟建立波流联合作用下的二维潮流泥沙数学模型，对阳西电厂附近海域的潮流场和泥沙场进行模拟研究，计算取水口及防波堤区域冲淤量，分析大风浪作用下泥沙运动规律，其成果不仅可加深对海岸泥沙输运物质输移的关键过程的认识，同时可为阳西电厂及类似工程设计及运行维护提供技术支撑。

1 数学模型与工程方案

1.1 基本方程

1.1.1 潮流数学模型

连续方程:

(1)

动量方程:

(2)

(3)

(4)

应力项Ax、Ay为包括水平黏滞应力、表面风应力、底部切应力和波浪辐射应力。模型的空间离散是使用单元中心有限体积法，时间差分格式采用显式迎风格式。

1.1.2 波浪数学模型

波浪数学模型模拟的物理过程包括：地形和海流空间变化导致的波浪折射作用；地形和海流空间变化导致的浅水变形作用；逆向流造成的障碍和反射作用；障碍物的阻挡或部分传播作用。其基本方程为：

(5)

N=N(σ,θ,x,y,t)=E(σ,θ,x,y,t)/σ

式中：N为波作用密度谱；σ为相对波频；θ为波向；Cx和Cy为波浪在x和y方向的传播速度分量；Cσ和Cθ为σ、θ空间的波浪传播速度；S是以谱密度表示的源项，包括风能输入、波与波之间的非线性相互作用和能量耗散。

该模型采用了全隐式有限差分格式对波作用动量守恒能量平衡方程进行离散，然后在4个象限中用迭代的方法进行求解，其计算是无条件稳定的，因而允许较大的时间步长。

1.1.3 泥沙数学模型

(1)基本方程。根据海区泥沙场分析结果可知，本海区泥沙粒径较小，表层沉积物以粉沙为主，泥沙输运以悬移为主，用二维悬沙输运模型来计算泥沙输运及泥沙骤淤。

黏性泥沙输运模型涉及泥沙在水体中的运动以及泥沙与底床的相互作用。悬移泥沙的输运一般建立在水动力模型中的对流项中，可用以下方程来描述：

(6)

(2)波流共同作用下的床面剪切应力。海岸地区的泥沙输移大多是在波浪和潮流共同作用下完成的，波浪主要起掀沙作用，但本身搬运泥沙能力较弱，而泥沙一旦处于悬浮状态，一般潮流就能搬运，也即前面提到的“波浪掀沙、潮流输沙”。所以计算时必须考虑波流共同作用下的泥沙运动。在模型中，床面剪切应力可采用以下公式进行计算。

波浪作用下的底床平均剪切力计算公式为：

(7)

式中：fw为波浪摩擦系数，0.047 (Swart 1974)；Ub为波浪水质点底部轨迹运动速度：

(8)

式中：Hs为有效波高；T为周期。

纯流作用下的底床平均剪切应力计算公式为：

(9)

(10)

式中：V为平均流速；fc为流阻系数：h为水深；k为底床糙度。

波流相互作用的底床剪切力计算公式(soulsby eta al.1993)为：

(11)

(12)

式中：τc为纯流作用底床剪切力；τw为纯波作用底床剪切力；τmax为波流相互作用最大底床剪切力；τmean为波流相互作用平均底床剪切力；a、b、m、n、p、q为常数,由波流理论来确定。

1.2 数学模型构建

1.2.1 计算范围

为了更好的提供水动力模型所需要的边界条件以保证局部流场计算符合潮流场的整体物理特征，潮流模拟采用嵌套计算。大模型范围包括整个海陵湾，小模型范围和网格均与波浪模型的范围和网格一致，2层模型的计算区域见图2。

图2 模型网格布置

根据实测海流的流向及大小，大模型东边界定在海陵岛下塘见村附近，在电厂以东约34 km，西边界定在茂名市电白县博贺镇北山村附近，在电厂以西约38 km，南边界定在厂址以南约47 km之处，模型包括整个海陵湾水域。模拟原体水域面积约2 560 km2。模型采用能够很好地贴合岸线和方便对工程区域加密的非结构三角形网格单元，模型区域共剖分了44 191 个网格，共有22 709 个节点，其中最小的网格单元边长约为50 m，厂址附近海域使用最小网格，可以较好地模拟岸线变化及拟建电厂平面布置。

1.2.2 边界选取

(1)波浪边界条件。本次计算选择本项目工在工程海域测波站2004 年1月-2005年1月的实测波浪统计资料，测波站的水深为18 m，与波浪模型的南边界水深基本一致，因此本次波浪计算可以直接采用测波站的波浪统计资料作为计算边界。根据波浪统计结果，本海区的常浪向为ESE、SE、SSE和S向，其年内出现的频率分别为25.3%、40.73%、13.27%和13.48%，因此进行波流耦合计算时，需分4组工况进行。各组波浪边界情况见表1。

表1 波浪模型计算边界

(2)潮流边界条件。模型计算区域为近岸水域，因此边界由潮位来给。鉴于本海域的潮汐比较复杂，本次大模型试验的边界采用调和常数换算成水位来试算，然后通过与实测潮位的不断验证、调整和耦合，最终得出一个能反映实测情况的潮位边界，最终应用此边界作为大模型的计算边界。大模型计算结果为小模型提供水位或者流量边界。

(3)泥沙边界条件。本次建立泥沙数学模型的目的是计算工程后该海域的年内冲淤演变情况，因此，计算水动力的选择必须按年来考虑。在泥沙模型计算过程中，潮流边界由调和常数计算出的水位提供，计算边界包含整个全潮过程。波浪计算边界则以实测统计资料为基础，以各个主浪向的波要素为边界条件，分别计算出各个主浪向波浪作用1 a的冲淤量，然后根据它们各自的年出现频率进行加权统计，最后得到年内的冲淤变化量。

1.2.3 模型参数取值

(1)水动力计算参数。模型海域的初始潮位取各边界潮位的平均值，初始流速取0。经反复调试，大、小模型中的涡黏系数取0.28，曼宁糙率取0.025，模型计算时间步长60 s，最小时间步长取0.01 s，CFL数取0.8。

(2)泥沙计算参数。

①海床床面泥沙级配。通过对沉积物样品的粒度分析表明，工程海域沉积物中值粒径分布总体上表现为东细西粗、南细北粗的特征，近岸及双山岛附近沉积物中值粒径大部分小于2 Ф(即d50>0.25 mm)，泥沙颗粒较粗，在5 m等深线以外海域中值粒径一般大于4 Ф(即d50<0.063 mm)。 沉积物类型以沙、黏土质粉沙为主，其他类型为砾石质沙、粉沙质沙、粉沙质黏土。各类型沉积物空间分布见图3。

图3 沉积物类型空间分布图

②泥沙参数。悬沙中值粒径d50=0.004 mm，由于泥沙颗粒较小，需考虑泥沙的絮凝作用，在考虑泥沙絮凝作用下的悬沙中值粒径d50=0.024 mm；絮团泥沙沉速w=0.000 45 m/s；泥沙沉降几率根据经验取值a=0.45；悬沙干容重取值r0=720 kg/m3。

1.3 数学模型验证

潮流数学模型采用2013年夏季的水文资料进行验证，2013年夏季小潮期实测时间为7月16日12∶00至7月17日13∶00；中潮期实测时间为7月19日12∶00至7月20日13∶00。大潮期实测时间为7月22日14∶00至7月23日15∶00。2013年夏夏季小、中、大潮潮位验证相对误差分别为0.044、0.053、0.058 m(见表2)；流速验证平均绝对误差控制在0.12 m/s以内(见表3)；从含沙量验证情况看(见图4)，计算含沙量与实测含沙量在总体趋势和峰值的波动变化与实测含沙量过程较为同步且吻合较好。验证结果表明该模型可用于本工程海域水动力与泥沙模拟。

表2 潮位误差统计 m

表3 流速误差统计 m/s

图4 夏季中潮悬沙验证曲线

1.4 工程方案

本次拟建工程包括电厂二期取排水工程、码头扩建工程。二期工程新建5～8号机组取、排水口，每台机组取、排水量均为78 m3/s。二期工程同时对码头进行扩建，码头的扩建方案为在一期卸煤码头的基础上向外延伸310 m，增加一个7 万t的泊位。

2 计算结果分析

2.1 工程扩建前、后流场变化

2.1.1 流态变化

工程附近海区的影响主要体现在3个方面：防波堤延长改变了原有的海岸线，使流场发生变化；防波堤西南侧新加的7、8号取水口改变附近流向，增大流速；电厂北侧排水口亦改变附近流向、增大流速，见图5。

图5 工程扩建前后流场对比

延长防波堤对海岸线改变比较大，因而对流态的影响也最为显著，防波堤两侧流速减小，堤头外由于挑流作用流速增大。

2.1.2 流速变化

图6为夏季大潮涨急、落急时刻的工程扩建前后流速变化等值线。由图6可见，工程建设以后，防波堤两侧流速明显减小，而防波堤端口东侧由于挑流作用，流速增大；取排水口出流速明显增大。工程实施后，港池内平均流速有小幅度降低，范围为3.1%～9.9%；受到防波堤对涨落潮流的遮挡作用，防波堤北面与南面流速大幅降低，范围为25.8%～60%；排水口处，受电厂排水影响，流速增大可达167%；排水口附近，但离排水口有一定距离，流速变化不大，为0.8%～10.9%，可见电厂排水对海域流速大小的影响范围较小；点9位于延长防波堤东侧，由于挑流作用，流速增大幅度可达28.7%；取水口处，受电厂取水影响，流速增大，增幅最大可达58%。

图6 工程扩建前后流速变化等值线

2.2 年回淤分析

年回淤量计算采用的波浪边界为实测统计值(见表1)，以各个主浪向的波要素为边界条件，分别计算出各个主浪向波浪作用1 a的冲淤量，然后根据它们各自的年出现频率进行加权统计，最后得到年内的冲淤变化量。计算中采用的潮流动力为典型的全潮过程。

泥沙数学模型计算的波浪、潮流共同作用下工程附近海域地形冲淤变化见图7[其中图7(a)为工程扩建前的冲淤结果，图7(b)为工程扩建后的冲淤结果]，图7中正值表示淤积，负值表示冲刷。

图7 波流共同作用下地形冲淤结果

由图7及表4可见，工程扩建前后工程海域的总体冲淤趋势基本一致，防波堤南、北两侧均表现为淤积，防波堤东侧为冲刷。工程后由于防波堤加长，东侧冲刷区相应的东移，冲刷的最大深度有所增加(增加约0.02 m/a)，而防波堤南、北侧的淤积区由于防波堤的加长而淤积强度稍有增强(约0.01 m/a)。

表4 工程特征部位的冲淤强度

工程扩建后，取水口东侧、南侧均为泥沙淤积区，取水口附近年冲淤厚度范围为0.1～0.2 m/a；加长码头防波堤后，堤头东侧水域则会有所冲刷，最大冲刷为0.2 m/a左右；而码头港池、泊位、回旋水域则稍有淤积，年淤积厚度为0.04～0.10 m/a。

2.3 骤淤计算

骤淤是恶劣天气条件(台风、风暴潮、热带气旋)下短时间内产生的泥沙淤积。根据工程海区的地质资料，厂址附近浅海区泥沙较粗，为沙粒，水深在10 m以上的外海区泥沙粒径较小为粉沙及黏土等，这种泥沙在风暴潮时很容易起动，水体含沙量将显著增大，有可能造成取水口及港池的大量淤积，因此必要进行骤淤分析。

本次骤淤计算本试验选取200 a一遇高潮位4.16 m、遭遇9713号台风(特大台风，S向，H4%=6.6 m(H有效=5.5 m)，且连续作用2 d作为台风骤淤的试验条件。

从图8及表5可见，工程扩建后，取水口东南侧均为泥沙淤积区，最大淤积厚度可达0.5 m，距取水口位置约1.6 km，取水口附近淤积厚度为0.20～0.25 m；加长码头防波堤后，堤头东侧会有所冲刷，最大冲刷可达0.5 m左右；而码头港池、泊位、回旋水域则稍有淤积，淤积厚度为0.2～0.5 m。

图8 9713号台风过后地形冲淤状况

3 结 论

本文建立工程海区波流共同作用下的平面二维泥沙数学模型，在利用工程海区实测水沙资料对模型进行验证的基础上，计算分析了波流联合作用下广东阳西电厂海域潮流泥沙过程及暴风骤淤特征，主要得出以下结论。

(1)阳西电厂二期工程方案实施后，流态变化比较大的区域主要是取排水口附近、防波堤两侧以及防波堤东侧水域。排水口处流速明显增大，增幅可达167%， 取水口处与防波堤东侧水域流速也增大，增幅最大达到58%，防波堤两侧流速则明显减小，降幅为25.8%～60.0%。

(2)阳西电厂二期工程方案实施后，取水口东侧、南侧均为泥沙淤积区，取水口附近年冲淤厚度范围为0.1～0.2 m/a；加长码头防波堤后，堤头东侧水域则会有所冲刷，最大冲刷为0.2 m/a左右；而码头港池、泊位、回旋水域则稍有淤积，年淤积厚度为0.04～0.10 m/a。

(3)台风作用下，海域含沙量明显增大，进入港池和取水口的含沙量也随之增大，因此可造成港池和取水口的较强淤积。骤淤条件下电厂取水口含沙量为0.08～0.10 kg/m3。取水口东南侧均为泥沙淤积区，最大淤积厚度可达0.5 m距取水口位置约1.6 km，取水口附近淤积厚度为0.20～0.25 m；加长码头防波堤后，堤头东侧会有所冲刷，最大冲刷可达0.5 m左右；而码头港池、泊位、回旋水域则稍有淤积，淤积厚度为0.2～0.5 m。