沿程非均匀分布植物区对海啸波消减的数值研究

蒋昌波，徐 进，邓 斌，陈 杰，屈 科

(1.长沙理工大学 水利工程学院，湖南 长沙 410114; 2.水沙科学与水灾害防治湖南省重点实验室，湖南 长沙 410114; 3.洞庭湖水环境治理与生态修复湖南省重点实验室，湖南 长沙 410114)

海啸引起的极端波况容易造成巨额的财产损失和众多的人员伤亡[1-2]。2004年的印尼海啸直接导致了超过23万人的伤亡[3]，12万栋房屋被摧毁[4]，经济损失高达140亿美元[5]。2011年日本海啸的最大波爬高接近40 m，也造成1.8万人伤亡[6]，导致300座桥梁和16.2万栋房屋的损坏[7]。此外，随着南中国海的大规模开发，南中国海海啸的研究也成为热点[8-9]。为了防范海啸的破坏作用，人们在实践中采用了多种工程措施，如建设防波堤和海墙。这些海岸保护措施，一定程度上起到了防浪作用[10]，但也会直接或者间接破坏海岸输沙平衡，致使岸滩演变加剧[11]，影响海岸的美观。而海岸植物作为一种自然屏障，既有利于提高海湾的绿化，也可以有效的减少海啸波的强度，降低海啸波对海岸区域的破坏[12-13]。随着“蓝色海湾整治行动”的持续推进，工程与生态措施协调发展成为关注的热点。因此，学者们日益关注沿海植物对海岸区域的保护机理[14]。

过去几十年，学者们通过理论分析[15-16]、物理试验[17-19]和数值模拟[20-22]开展了大量植物消波作用的研究工作，结果表明植物的消波特性主要受到水动力条件和植物分布特征的影响。以往学者们多将植物概化成为均匀分布的圆柱群，研究不同的入射波要素条件下均匀植物的消波特性。事实上，海岸带植物分布具有非均匀性，这种植物密度分布的非均匀性对植物消波特性影响方面的研究较少，非均匀植物区中的消减机理研究还不够。有学者通过现场观测[23-26]研究植物分布的非均匀性对植物消波特性的影响，但其高成本和复杂的现场条件在很大程度上制约研究的深入开展[27]。植物区的不均匀性主要分为沿植物区宽度方向(水平方向)的不均匀性和植物区高度方向(竖直方向)的不均匀性。Iimura等[28]在物理模型试验中组合两种不同水平分布密度的植物(密度大的在前或者密度大的在后)，发现孤立波经过不同分布的植物区后波高和流速接近。陈杰等[18]在试验中将红树林概化成根、茎、叶三层，研究竖直方向上不同分布对植物消波的影响。但是受到水槽尺度和试验模型布置等因素的限制，试验研究中的植物分布密度变化范围有限，难以系统揭示植物密度分布非均匀性对植物消波特性的影响。

近些年虽然有学者认为孤立波不能准确地模拟海啸波，但是其在传播中能够维持稳定的波形以及与长时间传播的海啸波有相似的形态[29]，因此目前仍然有许多学者用孤立波模拟海啸波。本文在研究非均匀分布刚性植物对海啸波的消减作用时使用孤立波代替海啸波。采用非静压单相流模型NHWAVE计算不同水平密度分布植物与孤立波的相互作用，研究它们消波特性的差异性。对比分析不同植物密度分布情况下反射波、波高和最大流速沿程衰减，并与均匀分布植物进行对比。

1 理论模型

1.1 控制方程

NHWAVE控制方程为σ坐标系下的连续性方程和动量方程，如式(1)和(2)所示：

(1)

(2)

采用k-ε湍流模型，其控制方程为：

(3)

(4)

1.2 拖曳力系数

由于研究非均匀分布植物对海啸的消波特性，需要根据不同的植物分布密度选取相应的拖曳力系数。根据Huang等[31]的试验研究，雷诺数Re的范围1.72×103～5.74×103，与文中计算工况的雷诺数处于同一区间，可将Huang等[31]得到的拖曳力系数经验公式添加到模型中，其具体表达式为：

CD=1.245+4.587φi/(1-φi)

(5)

2 模型验证

为了验证数值模型计算植物消波的准确性以及拖曳力系数经验公式在模型中的适用性，计算了孤立波和非淹没刚性植物的相互作用，并且与Huang等[31]的试验数据进行了详细对比分析，由于文献中流速数据的缺失，因此本文只对自由液面进行了验证。Huang等[31]的试验在长32 m、宽0.55 m的试验水槽中进行，植物区布置在水槽的中间。刚性植物由有机玻璃管制成，单株植物外径为0.01 m，试验中设置了A、B、C三种排列方式，如图1(a)所示，这三种分布方式所对应的植物分布密度φ分别为0.175、0.087、0.044，与此对应的单位面积植物株数N分别为2 228、1 108和560。在试验中共设置了三种不同的植物区长度W，分别为0.545 m、1.090 m和1.635 m，具体的试验布置以及波高仪分布位置如图1(b)所示。

图1 植物区三种排列方式和波高仪安放的位置

文中共验证5组不同算例，所有工况的水深h=0.15 m，波高H分别为0.03 m、0.04 m和0.041 7 m，计算工况如表1所示。建立二维数值水槽，计算区域长度为15 m，水平网格大小0.01 m，所有的工况水深均为0.15 m。为了确定竖直网格的数量，对竖直网格的收敛性进行验证，如图2所示。图2分别展示了A1工况下，测点G1和G5处分别在竖直方向设置10、15和20层网格时的水位对比，可认为竖直网格在20层时已收敛。验证结果如图3所示。图3(a)中G1数据列为G1波高仪测得水位的时间序列，其第一个峰值为入射波波高，第二个峰值为植物区的反射波，而G5为波高仪G5测得孤立波经过植物区的透射波数据列。由于工况B2中植物区密度较小，所以反射波也较小。比对图3(a)和(b)中的反射波，图3(c)植物区长度分别为0.545 m、1.090 m和1.635 m时的相对波高沿程变化，可看出计算结果与试验测量结果整体吻合度较高，表明数值模型可以有效模拟植物区的波浪传播变形。通过上述5组算例的验证可以看出Huang等[31]提出的经验公式对于文中数值模型是适用的，使用该经验公式算出来的拖曳力系数可以准确地模拟植物区对孤立波的反射、透射以及孤立波波高在植物区中沿程波高的衰减趋势，因此，使用Huang等[31]提出的经验公式(5)计算拖曳力系数CD是可靠的，可以用于研究不同密度分布植物对孤立波的消波作用。

表1 验证的工况

图2 网格收敛性验证

图3 工况A1、B2自由液面验证和C1、C2、C3沿程相对波高验证

3 工况设置

为研究不同密度分布植物对孤立波的消波特性，文中建立了二维数值水槽，计算区域网格设置与验证算例相同，所有的工况水深均为0.15 m，左侧为造波边界，右侧设置有长度为6 m的消波区，波高H设置为0.02～0.05 m。为满足计算稳定性的要求，计算均采用自适应时间步长，最大CFL数设置为0.2。

刚性植物区布置在X0=6.0 m和Xn(Xn=X0+W，W为植物区的长度)之间，其中植物区长度W设置为0.545 m、1.090 m和1.635 m，测点布置如图1(b)所示。在植物区设置了四种不同的密度线性变化方式，分别是：1)L1，植物密度从X0处的N0增大到Xn处的Nn；2)L2，植物密度从X0处的Nn减小到Xn处的N0；3)L3，植物密度在X0和X0+W/2之间由N0增大到Nn，在X0+W/2和Xn之间由Nn减小到N0；4)L4，植物密度在X0和X0+W/2之间由Nn减小到N0，在X0+W/2和Xn之间由N0增大到Nn。其中N0=560 m-2，Nn=2 228 m-2。图4(a)以W=1.090 m为例分别展示了L1、L2、L3、L4植物密度N沿程变化，图4(b)为沿程拖曳力系数，横坐标x为沿水平正方向与G3测点的距离，x=0 m即为G3所在的位置，所有计算工况如表2所示。

图4 W=1.090 m时植物区密度和拖曳力系数沿程变化

表2 计算工况表

4 结果分析

4.1 反射波对植物分布密度的响应机制

图5展示了植物区长度W=1.635 m时植物区前的反射系数随波高的变化，图中H1为G1所测波高，即入射波高，h为水深，CR为反射系数。图5反映出不同分布密度，植物区波浪的反射效应并不相同，反射系数从大到小依次是L2、L4、均匀分布、L3和L1，说明植物越集中在前面部分反射波越大。其中L1与L2的差值最大，H1/h=0.13时L2的反射系数比L1大80.2%，H1/h=0.33时L2的反射系数比L1大78.3%。从图5中还看出所有分布的反射系数均与入射波高线性相关，随着入射波高的增大而增大。

图5 W=1.635 m时反射系数随波高变化

图6 W=1.635 m、H1/h=0.33时L1和L2反射波随时间的变化

4.2 波高衰减对植物分布密度的响应机制

4.2.1 入射波高的影响

图7展示了W=1.090 m，入射波高H1/h分别为0.13、0.20、0.27、0.33时孤立波在植物区沿程波高衰减。横坐标x为沿水平正方向与G3测点的距离，x=0 m即为G3所在的位置。对于不同的波高，孤立波在植物区中沿程消减规律相近，L1与L2的曲线组成类“0”字形，L3和L4的曲线组成类“8”字形，“8”字形的交点均在x/h=4.10～4.20之间，而均匀分布的曲线在“0”字和“8”字的中间。造成这种现象的原因是，L1和L2在植物区密度分布相反，L1的密度线性增大，在植物区最前端密度最小，所以植物区前端波高衰减的强度弱，而L2的密度线性减小，在植物区最前端密度最大，所以植物区前端波高衰减的强度大，但是随着密度沿程的变化L1波高衰减的强度逐渐大于L2，最终在植物区尾端波高之间的差小于1.5%。L3在植物区中点处密度最大，向两端线性减小，消波效率先增大后减小，而L4与此相反，在中点处密度最小，向两端线性增大，所以在植物区中点后会得到相同的波高即两条曲线的交点，最终L3、L4、L1、L2、均匀分布的透射波高十分接近，其之间的差距小于2.5%。

图8展示了W=0.545 m、1.090 m和1.635 m时透射系数随波高的变化，从图中可看出透射系数随着入射波高的增大而减小。而L1、L2、L3、L4和均匀分布的透射系数随波高变化的趋势是一致的，因此在植物区长度和波高相同的情况下5种分布的透射系数相同。

图7 W=1.090 m时波高沿程变化

图8 透射系数随波高的变化

4.2.2 植物区长度的影响

图9(a)、图7(d)、图9(b)分布展示了H1/h=0.33，植物区长度分别为0.545 m、1.090 m、1.635 m时植物区中孤立波沿程波高。从图中可看出随着植物带长度的增大，植物区中波高沿程消减的机制没有明显的变化，L1、L2、L3和L4分布的沿程曲线仍然呈现“0”字形和“8”字形，而均匀分布的曲线位于它们中间。可以发现W=0.545 m和W=1.090 m时L3和L4的交点分别在x/h=2.08和x/h=6.18处，皆在植物区中部偏后的位置。

图10展示了H1/h=0.33时，透射系数随植物区长度的变化，图中横坐标W为植物区长度。从图中可看出透射系数均随着植物区长度的增大而减小，但是4种非均匀分布和均匀分布植物的透射系数随植物区长度变化的规律相同，因此出现了当植物区长度相同时5种分布具有相近透射系数的情况。

图9 H1/h=0.33时波高沿程变化

图10 H1/h=0.33时透射系数随植物区长度的变化

4.3 植物区中流速对植物分布密度的响应机制

4.3.1 入射波高的影响

图11分别展示了W=0.545 m、1.090 m和1.635 m时流速透射系数UT随波高的变化。4种非均匀分布和均匀分布植物的流速透射系数UT均随着入射波高的增大而减小，例如W=1.090 m时，H1/h=0.33的流速透射系数只有H1/h=0.13的64.1%。然而5种分布的流速透射系数在随着波高变化的过程中保持一致性，即UT随波高变化呈现出一样的规律，在植物区长度和波高相同时具有相近的流速透射系数。

图11 流速透射系数随波高的变化

图12展示了植物区长度W=1.090 m，入射波高H1/h分别为0.13、0.20、0.27、0.33时，植物区沿程最大流速的变化，图中的Umax,1为入射波的最大流速，取波峰经过断面时该断面上的最大流速。从图中可看出沿程最大流速的曲线也呈现出类似的“0”字形和“8”字形，只是其在最前端出现开口，这与上述的沿程波高曲线不同。出现不同的原因是图7和图9中的波高在x=0处由于反射等原因会使得瞬时波高大于入射波高，而且x=0处的密度越大瞬时波高越大，因此会在波高衰减的过程中出现两两交叉的情况。但是最大流速则不同，在x=0处的植物密度越大会使得最大流速越小，所以在消减过程中并不会出现交叉。通过图12(a)～(d)的对比可发现入射波高的变化，对最大流速的沿程变化规律没有明显的影响。

图12 W=1.090 m时沿程最大流速变化

4.3.2 植物区长度的影响

图13(a)、图12(d)、图13(b)分别展示了H1/h=0.33，植物区长度分别为0.545 m、1.090 m、1.635 m时植物区沿程最大流速的变化。从图中可看出对于不同的植物区长度，最大流速沿程消减的机制没有明显的变化。W=0.545 m、1.090 m和1.635 m的情况下，L1和L2的曲线呈现出开口的“0”字形，而L3和L4的曲线组成开口的“8”字形，均匀分布位于它们中间，而且最终经过植物区后最大流速相近，最大差为5.4%。

图14展示了H1/h=0.33时，UT随植物区长度的变化。从图中可看出，UT均随着植物区长度的增大而减小，例如W=1.635 m时的流速透射系数只有W=0.545 m时的51.2%。然而4种非均匀分布和均匀分布的流速透射系数随植物区长度的变化规律是一致的，在植物区长度相同的情况下5种分布具有相近的流速透射系数。

图13 H1/h=0.33时最大流速沿程变化

图14 H1/h=0.33时，流速透射系数随植物区长度的变化

5 结 语

采用非静压数值模型开展了非均匀刚性非淹没植物对海啸波消减的数值模拟工作。通过与Huang等[31]的试验数据进行对比验证数值模型的计算可靠性。文中研究了5种不同密度分布植物与孤立波的相互作用，分析了不同密度分布植物对海啸的消波特性。通过对计算结果进行分析，主要结论如下：

1)在植物总株数相同的情况下，植物区的非均匀性对孤立波的消减作用几乎没有影响，只改变了沿程的消波强度分布。对于平均密度相同的5种分布，孤立波波高和流速在不同分布密度的植物中消减过程不同，L1与L2沿程波高和最大流速曲线组成类“0”字形，L3和L4沿程波高和最大流速曲线呈类“8”字形。不同分布密度的植物区沿程消波强度分布不同，但最终透射波的波高和流速几乎相同，不同分布的透射波高之差小于3%，最大流速的差小于5.4%。在所有的工况下，孤立波的透射波高和流速的差皆很小，说明透射波的能量相差较小，可以说在植物总株数相同的情况下，不同分布密度的植物对孤立波的消减作用是相同的。

2)孤立波的入射波高和植物区的长度变化，不改变波浪在非均匀分布植物区沿程消减机制，波高和流速沿程衰减曲线皆呈现出类“0”字形和类“8”字形。而且孤立波经过不同分布植物区的透射系数和流速透射系数随着入射波高和植物区长度变化的规律是一致的。

3)通过研究可见，在植物总株数相同的情况下，植物区前端密度较大时会引起海啸波反射增强，但却对最终的海啸波消减效率影响有限。因此在实际工程实践中，在植物总分布密度接近的情况下，可使植物带密度随与海岸的距离逐渐减小。这样不但可以有效的消减海啸波的破坏作用，并且可以进一步减少海啸波反射引起的次生灾害。