非结构化网格下大范围波生流模拟和应用

王 平，张宁川

(大连理工大学海岸和近海工程国家重点实验室，辽宁大连 116024)

波浪在向近岸海域传播的过程中，由于地形的影响会发生变形、破碎等现象，伴随着质量、动量和能量的输运及传递，进而形成强度和影响范围都较大的近岸流。在近岸海域，波生沿岸流的现象非常普遍，尤其在浅水和河口区域更加显著，对近岸地貌变迁、泥沙输运以及近海污染物输运等都有重要的影响。故在近岸潮流、泥沙和污染物输运的模拟时考虑波浪的影响必不可少。

目前波生沿岸流数值模拟方法主要有两种，一种是基于Boussinesq方程的波生流模型［1］，该模型将波浪和流场的计算融合到一个控制方程中;第二种是建立在Longuet-Higgins和Stewart［2］提出的波浪辐射应力理论的基础上，该方法是将波浪场计算得到的辐射应力等波浪影响参量作为驱动因素添加到潮流场的控制方程中，很多学者对该方法进行了大量的数值模拟和研究，得到了一定的结果。已有的基于辐射应力理论模拟近岸波生流的方法，较多采用有限差分法对控制方程在矩形网格下进行数值离散，如包西林［3］、Sun［4］、Yoon［5］和 Tang［6］等，但实际工程中边界多为不规则形状，矩形网格对其拟合的精度较差。Wu［7］、李孟国［8］和唐军［9］等建立了基于三角形网格的波生流模型，但其波浪场都是基于缓坡方程，计算时网格步长受制于波长的限制，造成计算大范围波浪场时的效率较低，且不能和流场同步耦合。Choi［10］、朱首贤［11］、王彪［12］等将SWAN模型模拟得到的辐射应力添加到潮流模型中得到一种考虑波浪作用的流场模型，该方法解决了波生流计算时波浪场和潮流场不同尺度的问题，使得波浪和流场可以同步耦合，但受模式的限制在模拟中均未考虑波浪紊动的影响，同时在非结构化网格中SWAN模型并未考虑绕射的影响。

从含流的缓坡方程［13］出发，推出光程函数方程，联立光程函数方程和波作用守恒方程建立了基于非结构化网格下考虑绕射效应的波浪模型，将波浪场模拟得到的辐射应力和紊动系数等同步添加到流场模型(FVCOM［14］)中得到基于非结构化网格下的大范围波生流数值模式。其中波浪场模型的计算不受波长的限制，适用于近岸大范围波生流计算，同时使得波浪场和流场可以同步耦合计算。对所建立的模型进行了多个算例验证，数值模式的结果与实测吻合较好，验证了模型精度、可靠性以及对复杂边界的适应性。同时将模型应用于大连近岸区域，对常见浪在大连近岸区域产生的波生流进行了研究，验证了模式对大范围区域波生流计算的适用性。

1 数值模式

1.1 近岸波浪场控制方程

Kirby［12］推导出的含流非定常波浪方程，其表达式为:

其中方程左端第三项和第四项为波能密度在θ方向和σ方向的传播，右端为源汇项，可以包括波浪破碎，底摩阻等对能量的影响。利用波数矢K替代波数k的求解，可以计算绕射对波高的影响。改进后波速、波群速以及波能在θ和σ方向的传播速度分别为:

联立式(2)、式(4)建立包含绕射的波浪计算模式，该模式在空间网格上的划分不再受波长的限制，可以用于大范围波浪传播计算;同时利用有限体积法对控制方程进行离散，离散基于非结构化网格可以和目前比较通用的流场模型同步耦合，且对复杂的岸线能很好的拟合;而且光程函数方程可以有效改善模型在计算近岸复杂海域的波浪绕射效应［16］，提高了模型计算近岸波浪传播的精度。

1.2 近岸波生流控制方程

流场模型采用修改后的三维水动力模型FVCOM，其水动力方程采用在ζ坐标系下的三维Navier-Stokes方程组，同时引入三维波生时均剩余动量和波浪的紊动掺混效应对水体的影响，方程基本形式见式(5)、式(7):

式中:η为自由水面;x和y为水平坐标;ζ=(z-η)/D为转换后的垂向坐标;w为ζ坐标系下的垂向流速分量;u和v为水平流速分量;Kmc为波流共存的垂向紊动粘性系数;Amc为波流共存的水平紊动粘性系数;R为波浪辐射应力项。

对波流共同作用下的紊动系数采用文献［17］的方法，即分别单独求解水流与波浪引起的紊动系数A和K，并将其线性叠加，可表述为:

波浪引起的紊动系数Aw和Kw计算均采用Larson-Kraus公式计算，Aw=Kw=λumaxH=2λH2cosh(1+σ)kD)/Tsinh(kD)，式中umax为波浪底部质点最大流速;λ为无因次系数。

水流引起的水平紊动系数采用Smagorinsky公式(10)。式中C为常数;Ω为控制单元体面积。水流引起的垂向紊动系数Km采用Mellor-Yamada紊流闭合模型求解。

FVCOM模型采用外模和内模分别计算的方法，其中外模为沿深度平均的二维流场计算，内模为沿深度分层的三维流场计算。计算波浪对流场的影响时，在外模和内模均考虑辐射应力项。波浪辐射应力分项采用波能E等参数计算，二维辐射应力形式具体计算公式如下:

三维辐射应力采用郑金海［18］推导的公式:

上述公式将辐射应力的计算沿深度分为三个计算区域，R01、R02和R03的具体计算形式见文献［18］。

流场和波浪的耦合过程分为，流场和波浪单独计算以及参量相互传递两个过程。其中流场为波浪提供流速和水位参量，而波浪场则为流场提供辐射应力项和紊动系数。流场和波浪场在完成每一步计算后以及下一步计算开始之前进行参量传递，最终实现波流的相互耦合计算。

1.3 控制方程的离散

基于三角形网格，波浪场采用的有限体积法对波作用守恒方程进行离散;对光程函数方程采用网格中心格式的有限体积法离散。时间离散采用欧拉向前格式，空间采用格林公式，具体形式如下:

对式(2)采用网格顶点格式的有限体积法在控制体内积分，空间求导采用将面积分转化为线积分，并将线积分写成各控制边求和的形式，整理得到离散后的方程为:

式中:MT(i)为第i个控制体的边数，对三角形控制体取为3;Ai为控制体面积;标注m为控制体第m边上的参量;δ表示控制体边上参量在x或y方向上的梯度，采用上述格林公式计算;ζ表示参量在时间上的梯度，采用向后差分的一阶格式离散。

对式(3)分四步离散求解，空间离散采用单元中心顶点格式的有限体积法在控制体内积分;对波谱在频率方向上的离散采用FCT［19］离散方法;对波谱在传播方程上的离散采用二阶隐身Crank-Nicolson差分方法;对源项对波谱的影响采用二阶隐式中心差分离散求解。离散后见式(15a)～(15d)。

式中:NT(i)为第i个控制体的边数;Ai为控制体面积;标注m为控制体第m边上的参量;FCT中的φ1和A*计算公式见文献［19］;Cgx，i，m和 Cgy，i，m为控制体边 m 上 U+(Cg/k)K 在 x，y 方向上的分量。

流场控制方程的离散见文献［14］，对于增加的辐射应力项同样采用单元中心顶点格式的有限体积法在控制体内积分，同时用格林公式将面积分转换为各控制边的求和形式，具体形式见式(16a)～(16b)。

2 数值模型验证及应用

2.1 沿岸流计算及实验验证

波浪斜向入射近岸时，破碎引起的辐射应力会引起水体的顺岸移动而发生沿岸流。王淑平［20］进行了沿岸流的模拟实验，实验分别研究了1∶40和1∶100两种坡度下波浪斜向入射产生的沿岸流，实验工况见表1。

表1 沿岸流实验计算工况参数Tab.1 The calculation parameters of longshore current experiment

其中D为波浪入射处水深，H0为入射波高，T为入射波周期，γ为波浪破碎指标，Bfric为波流作用下的底摩擦系数。波浪沿正向入射，图1给出了两种工况计算稳定后的波生流平均流场图，图2和图3给出了2种计算工况下的波高、增减水和沿岸流平均流速的数值解和实测值的比较。从计算结果可以看出，所建立的模式可以很好地模拟波浪斜向入射近岸时引起的沿岸流现象。

图1 两种工况下波生沿岸流流场Fig.1 The wave-induced longshore current velocity field of two cases

图2 工况1数值解与实测值对比Fig.2 Comparison between the numerical and experimental results in case 1

图3 工况2数值解与实测值对比Fig.3 Comparison between the numerical and experimental results in case 2

在破碎带内外选取4个点得到其流速的垂线分布，见图4，Dx为点距岸边的垂线距离。从图中可以看出波生流的垂线分布区别于潮流的近对数分布形态，而垂向分布较为均匀。根据Visser［21］的研究，这种特征与波浪引起的附加垂向掺混有关，其中较大的附加粘性可使流速垂向剖面更加均匀。

图4 两种工况下四个点的流速垂向分布Fig.4 The vertical velocity distribution of four points in two cases

2.2 人工岛后波生环流的计算及验证

波浪传播过程中遇到人工岛等障碍物时，会在人工岛后发生绕射;由于波高的变化从而产生辐射应力，在辐射应力的驱动下会形成岛后环流现象。日本东京电力研究所于1994年利用大型水池实验研究了柏崎刈羽核电站人工岛内侧取水口周围波浪场以及波浪破碎而引起的波生流［22］。

试验中人工岛布置在坡度为1∶50的斜坡上，整个实验区域长为23.5 m，宽为19.2 m;人工岛大小为6.6 m×3.3 m。模型计算域选取19.2 m×15.0 m，入射波高为0.03 m，周期为1.0 s，波浪破碎指标选为0.76，流场底摩擦系数选为0.025。

图5给出实验得到的波生流实测流场分布，表2给出了A～J10个点的平均流场实测值和计算值对比，图6给出了波生流模型的计算网格和平均流场分布计算结果。根据对比可以看出本模型流场流速的数值解和实测值比较接近;模拟得到在人工岛内侧的对称轴附近产生了回流现象，同时在岛内侧区域形成一对方向相反的环流现象，该形态和实测得到人工岛后的波生流结果一致。

图5 人工岛后波生流实测流场Fig.5 The measured result of wave-induced current velocity field behind the artificial island

表2 人工岛后波生流场沿深度平均流速的计算值和实测值对比Tab.2 Comparison between the numerical and experimental results of the average velocity along the depth

图6 人工岛后波生流计算网格和流场计算结果Fig.6 Computational grids and wave-induced current velocity field

2.3 连续凹槽地形上的波生流计算及验证

Dingemans［23］和李孟国［8］等在验证近岸波生流模型时，均采用波浪在正向和斜向入射连续凹槽地形来验证模式的适定性。连续凹槽海区的水深可以用以下表达式表示:

模型的计算水深图和计算网格见图7，分别计算波浪正向入射(T=4.0 s，H0=0.92 m，θ0=90°)和斜向入射(T=4.0 s，H0=1.0 m，θ0=63°)两种工况。李孟国模拟得到的波生流场见图8。

图7 模型计算网格和水深Fig.7 Computational grids and the bottom topography

本模式模拟得到的平均波生流场见图8，得到波浪在正向入射连续凹槽地形时，会在凹槽的两边产生相向的沿岸流，沿岸流汇聚到凹槽顶点时会形成水体的雍高，同时雍高的水体则会通过破碎带回流海洋从而产生裂流;波浪在斜向入射连续凹槽时则会产生顺凹槽形状的沿岸流。对比图8和图9可知本模式模拟的结果和已有研究结果较吻合。

2.4 大连湾区域的波流耦合计算

由于本波浪模式的空间步长不再受制于波长的限制，波浪模型计算与流场计算可以共用一套网格，波浪和流场模型同步计算，每步计算完成后将流场实时得到的水深、流速提供给波浪场;同时将波浪场得到的辐射应力、紊动系数等反馈给流场，即实现了波浪场和流场之间参量的双向传递和同步耦合。

模型计算选择大连湾及附近海域作为研究对象，其中波浪采用不规则波，不规则波的能量采用JONSWAP谱计算。图10(a)和10(b)分别为计算区域地形和计算网格图。波浪采用S和SE向入射(为大连湾区域的常见浪向)，波高为2.0和4.0 m，周期为8 s;计算中波浪和流场的步长均选为5 s，最小空间步长为100 m。

模型分两步计算完成。首先计算波浪场，等波浪场计算稳定后，流场开始计算，并将波浪场计算所得的辐射应力、紊动系数等参量添加进流场;此时波浪和流场同步计算直至波生流场达到稳定形态。由于波生流多发生在近岸区域，模型选取图10(a)中方框中的波生流场分析。

图8 李孟国模拟得到的波生流场Fig.8 The wave-induced current velocity field simulated by LI Meng-guo

图9 本模式模拟得到的波生流场Fig.9 The wave-induced current velocity field computed by the present model

图10 大连湾海域分布及模型计算网格Fig.10 The area distribution and computational mesh of the Dalian Bay

图11(a)和11(b)为波浪S向入射下波高分别为2.0 m和4.0 m产生的波生流平均流场;图12(a)和12(b)为波浪SE向入射下波高分别为2.0 m和4.0 m产生的波生流场。从图可知波浪斜向入射破碎时会在波浪破碎处至岸边区域发生沿岸流，沿岸流的方向同波浪入射成小角度;当遇到凹槽地形时，若产生两股相对的沿岸流时(图11)，则会在相汇处形成逆流;波高越大，波浪破碎处离岸越远，沿岸流的影响范围就会越大，强度也会越大;同时沿岸流的强度和近岸潮流强度为同一数量级，在考虑近岸污染物、泥沙输运时，波生流为主要影响因素之一。

图11 S向波浪入射在近岸产生的沿岸流Fig.11 The wave-induced current velocity field by the south incident wave

图12 SE向波浪入射在近岸产生的沿岸流Fig.12 The wave-induced current velocity field by the south-east incident wave

图13为4.0 m有效波高的波浪S向和SE向入射后，形成的近岸增减水图。从图中可知波浪在近岸破碎时不仅会形成沿岸的波生流，同时在波浪破碎区形成增水，而在波浪未破区域形成减水。

从上述结果可以得到，波浪对流场的影响主要分为两个方面，一是会在近岸破碎带内产生沿岸流，二是会在近岸形成增减水;而流对波浪的影响主要表现为，波流同向时波高减小、周期变大，而波流逆向时波高增大、周期变小，当波浪与水流斜向相交时，波浪则会发生折射现象。

图13 S向和SE向波浪入射后在近岸形成的增减水Fig.13 The distribution of wave set-up and set-down by the south and south-east incident waves

3 结语

在非结构化网格下基于光程函数方程和波作用守恒方程，建立了考虑绕射效应的大范围波浪计算模型;利用波浪模型计算得到的辐射应力、波浪紊动系数等参量以及三维水动力模型(FVCOM)建立了非结构化网格下的大范围波生流数值模型。模型分别对波浪斜向入射产生的沿岸流、人工岛后的波生环流以及连续凸起地形下的逆流进行了模拟，数值结果和实测结果对比表明，本波生流数值模式可以高效精确地实现近岸波浪破碎产生的近岸流的模拟计算。

利用本模型对实际海域的计算结果，表明该模型可以用于大范围的波流同步耦合计算;同时非结构化网格可以很好地拟合复杂岸线的变化;且说明波生流为近岸水动力不可缺少的影响因素之一。由于缺少实际海域波生流的实测资料，本模式对实际海域波生流的模拟计算需要进一步的验证。

［1］ 卢 吉，余锡平.基于Boussinesq方程的近岸波流统一模型［J］.水动力学研究与进展:A辑，2008，23(3):314-320.

［2］ Longuet-Higgins M S，Stewart R W.Radiation stress and mass transport in gravity waves with application to surf beats［J］.Journal of Marine Research，1962，13:481-504.

［3］ 包西林，西村仁嗣.Hardy-Cross法波生流数值解析［M］.北京:海洋出版社，2006.

［4］ Sun T，Tao J H.Experimental and numerical study of wave-induced long-shore currents on a mild slope beach［J］.China Ocean Engineering，2005，19(3):469-484.

［5］ Yoon S B，Cho Y，Lee C.Effects of breaking-induced currents on refraction-diffraction of irregular waves over submerged shoal［J］.Ocean Engineering，2004，31(5-6):633-652.

［6］ Tang J，Shen Y M，Qiu D H.Numerical study of pollutant movement in waves and wave-induced long-shore currents in surf zone［J］.Acta Oceanologica Sinica，2008，27(1):122-131.

［7］ Wu C S，Liu P L F.Finite element modeling of nonlinear coastal currents［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1985，111(2):417-432.

［8］ 李孟国.波生近岸流的数学模型研究［J］.水道港口，2003，24(4):161-166.

［9］ 唐军，魏美芳.非结构化网格下近岸波生流数值模拟［J］.海洋学报，2010，32(6):41-46.

［10］ Choi J，Lim C，Lee J，et al.Evolution of waves and currents over a submerged laboratory shoal［J］.Coastal Engineering，2009，56(3):297-312.

［11］朱首贤.流、浪模式和物质长期输运分离研究［D］.上海:华东师范大学，2005.

［12］王 彪，沈永明，王 亮.长兴岛海区波流相互作用数值模拟研究［J］.海洋工程，2012，30(3):87-96.

［13］ Kirby J T.A note on linear surface wave-current interaction over slowly varying topography［J］.J.Geophys Res.，1984，89:745-747.

［14］Chen C，Liu H，Beardsley R C.An unstructured grid，finite-volume，three-dimensional，primitive equations ocean model:Application to coastal ocean and estuaries［J］.J.Atm.＆ Oceanic Tech.，2003，20:159-186.

［15］ Mei C C.The Applied Dynamics of Ocean Surface Waves［M］.New York:Word Scientific Co.Pte.Ltd.，1983:740.

［16］ Hong Guang-wen.Mathematical models for combined refraction-diffraction of waves on non-uniform current and depth［J］.China Ocean Eng.，1996，10(4):433-454.

［17］ Xia H Y，Xia Z W，Zhu L S.Vertical variation in radiation stress and wave-induced current［J］.Coast.Eng.，2004，51(4):309-321.

［18］郑金海，严以新，彭世根.波浪剩余动量流垂向分布研究［J］.河海大学学报，2000，28(1):8-13.

［19］ Boris J P，Book D L.Flux corrected transportⅠ，SHASTA，a fluid transport algorithm that works［J］.J.Comp.Phys. ，1973，11:38-39.

［20］王淑平.沿岸流的研究［D］.大连:大连理工大学，2001.

［21］ Visser P J.Laboratory measurements of uniform longshore currents［J］.Coastal Engineering，1991，15:563-593.

［22］李绍武，柴山知也.近岸区人工岛周围波生流系统数学模型的验证［J］.天津大学学报，1998，31(5):584-589.

［23］ Dingemans M W，Radder A C，De Vriend H J.Computation of the driving forces of wave-induced currents［J］.Coastal Engineering，1987，11:539-563.