港内地形对N波诱发的港湾共振影响研究

2021-05-08宋向荣周校军高俊亮王岗

哈尔滨工程大学学报 2021年3期

宋向荣，周校军，高俊亮,2，王岗

(1.江苏科技大学船舶与海洋工程学院，江苏镇江 212003; 2.江苏科技大学江苏省船舶先进设计制造技术重点实验室，江苏镇江 212003; 3.河海大学港口海岸与近海工程学院，江苏南京 210098)

港湾共振是指当传入港内的低频波浪的频率与港湾的本征频率相同时，港内水体会出现大幅波动的现象[1]。港湾共振可以诱发港内泊船的剧烈运动，显著影响系泊安全和货物的装卸[2]。甚至可能造成人员伤亡和重大的经济损失。夏季由于受太平洋上的热带风暴影响，我国东南沿海城市的很多港口也曾受到港湾共振的影响[3]。诱发港湾共振的原因有很多，其中主要能够诱发港湾共振的有气象波、次重力波、地震波、海啸、海洋内波和剪切流等[4]。

为了对海啸波诱发的港湾共振进行预报和研究，Sobey[5]提出了正交模态分解法，港湾的本征频率、共振模态形状及各模态的响应幅值可以由这一方法求得。但该方法在全反射边界条件的处理上不够精确。Gao等[6]提出镜像法对其进行了改进，显著提高了预测结果的精度。随后Gao等[7]使用改进的正交模态分解法并结合Boussinesq数值模型系统研究了在不同入射波幅和港底剖面形状条件下孤立波对港内波幅演化、共振波幅及港内总波能和相对波能分布的影响。Gao等[8-9]进一步使用改进的正交模态分解法和Boussinesq数值模型研究了连续孤立波和N波的波形参数对港内最大爬高及港内相对波能分布的影响，上述研究中港内均设置为平底地形。随后，郑子波等[10]引入折线型地形，系统研究了该地形条件对孤立波诱发的港内最大爬高及波能的影响。Gao等[11]利用Boussinesq数值模型研究了高模态下双色短波群在近海岸礁地形附近引起的狭长港内低频振荡。

虽然绝大多数学者采用孤立波来研究海啸波[12-13]，但人们从许多实际观察中发现大多数传至海岸的海啸波波形像“N”，可以看作由一个大波峰和一个大波谷组成，这类海啸波被称为N波[14-15]。然而，目前对由N波诱发的瞬变港湾共振现象的研究还极为有限。因此，本文进一步系统研究不同港内折线型地形对N波入射诱发港湾共振条件下的港内最大爬高及波能分布的影响。另外，不同N波入射波幅对于港内波况的影响也被研究。

1 数值模型和分析方法

1.1 FUNWAVE-TVD数值模型

FUNWAVE-TVD 是Shi等[16]开发的完全非线性Boussinesq波浪模型，该模型使用混合有限元有限差分格式求解了Chen等[17]提出的完全非线性Boussinesq方程，并且结合了Kennedy等[18]提出的动参考面理论。该模型在频散中保留到O(μ2)，在非线性a/h保留了所有阶(这里μ=h/l，h为局部静水深，l为波长，a为波幅)。其控制方程表达式为：

ηt+·M=0

(1)

V2-V3-R}+gηh

(2)

式中：g、η和H=h+η分别表示重力加速度、局部表面位移和总局部水深；=((∂/∂x),(∂/∂y))为水平梯度运算符；M为水平体积通量：

(3)

(4)

式中：A=·(huα)，B=·uα。V1和V2为色散项。

(6)

V3为竖直涡量场的O(μ2)分布：

(7)

其中：

ω0=(×uα)·iz=vα,x-uα,y

(8)

zα,x(Ay+zαBy)-zα,y(Ax+zαBx)

(9)

式中：iz为竖直方向的单位向量；R为包含底部摩擦和次网格横向紊流混合的耗散和扩散项。

FUNWAVE-TVD通过一个高阶激波捕捉总变差下降(total variation diminishing，TVD)格式求解，这使得该模型不依赖经验公式就可以模拟波浪破碎并结合波能耗散[15]。此外，该模型采用保持强稳定性的三阶Runge-Kutta格式来自动调整时间步长，通过带有非阻塞通信的消息传递接口(MPI)实现并行运算。通过这些改进，该模型在预测近岸波浪过程时变得更加稳定，这些波浪过程包括波浪浅水、折射、衍射、破碎及在平台和自然沙滩上的爬高和下冲[15]。由于Gao等[9]对FUNWAVE-TVD模型在模拟N波诱发的瞬变港湾振荡的适用性进行了验证，故本文没有再对该模型进行验证。

1.2 分析方法

本文各实验均使用正交模态分解法来对港内波况进行分析处理。该方法由Sobey[5]提出，用来预测港湾的本征频率和与之对应的本征模态以及确定海啸诱发的港湾内各共振模态的响应幅值。该方法包括2个计算步骤。第1步为计算港口的本征频率和模态形状。第2步通过一个多维优化问题，对由海啸波所诱发的港湾共振过程中的各模态响应幅值进行预测。本征频率和模态形状将作为已知量参与到第2步的预测中。Gao等[6]随后对该方法进行了改进以提高本征频率和模态形状的计算精度。由于正交模态分解法是基于线性理论得到的，只有在港内波浪非线性参数较小的情况下才能精确地预测港内各模态的响应幅值[13]。故本文只针对入射N波波幅较小的实验进行了分离，然后再研究各实验的港内波能分布。

2 数值实验

图1呈现了所有数值实验所使用的狭长型矩形港口俯视图、港内各地形的剖面图以及坐标系的布置。所有港口的长度为1 500 m，宽度2b=30 m；坐标系的x、y轴布置在静水面上，x轴在沿港口长度方向的中心线上，y轴在港口口门处，z轴竖直向上。本文所研究的港内地形为折线型，分别为凸型底、凹型底、平台型底和斜直线型底。对于凸型底与凹型底地形而言，有且只有一个折点，且折点位于沿港口长度的垂直平分线上，各组中折点处的水深有所区别。对于平台型地形来说，可以将其看作是具有2个折点的折线型地形，2个折点处的水深一致并且关于港口长度垂直平分线对称。图中h1、h0、hp和l分别代表港口口门处水深、港口底墙处水深、平台处水深和平台的长度。在本文中，港口口门处的水深和外海水深均为12 m，底墙处的水深为3 m，平台处的水深为7.5 m；各组实验的底部剖面形状各不相同。港内的平均水深为：

(10)

图1 不同底部剖面形状的狭长型矩形港口概图Fig.1 Elongated harbors with different bottom profile shapes

图2 数值波浪水槽示意Fig.2 Schematic diagram of numerical wave flume

3 结果和分析

3.1 波面时间序列

图3显示了A组、C组、E组实验中入射N波波幅A0=0.04 m时测点G01、G76和G151处的波面时间序列，为了更直观地观察波面时间序列所以通过A0对其进行了无因次化处理。Ani和Anr分别表示的是入射和反射波波幅，其中n代表的是测点的位置。由于N波入射波幅均相同，港口口门外各水体条件也均一致，所以当波浪传播到测点G01处时各组的波幅基本均于t1=44.28 s达到峰值。随后入射N波逐渐向港内传播，因为港内地形不同，各组的波面时间序列也会有不同的变化。A～E组的平均水深逐渐减小，港内波浪的平均传播速度也逐渐减小。故从t3开始，A～E组的波幅达到峰值的时间逐渐增大，且下冲的极值受平均水深的影响比爬高要更大一些。由于底墙不可穿透，因此t5和t6时刻对应的是港内最大的爬高和下冲。例如，由于浅水效应，t5从A组的209.48 s逐渐增加到E组的249.16 s，增幅为39.68 s。对于波峰而言，入射N波的放大因子Am/A0逐渐从A组的4.171增加到E组的4.274；对于波谷而言，放大因子Am/A0的变化规律不是很明显。在波峰到达底墙之前，不管是入射还是反射，越靠近底墙处波峰与波谷的放大因子Am/A0之差越大。通过观察G01和G76测点处的波面时间序列，不难发现经过底墙反射后波浪的波峰和波谷的幅值均小于入射时波浪的波峰和波谷的幅值。这是因为随着波浪向港内传播，港内水深均逐渐减小，波浪越往港内传播,其发生的反射越多，波能在不断的被消耗，故反射回来的波幅较小。

图3 A、C、E组中A0=0.04 m时测点G01、G76和G151处的波面时间序列Fig.3 Time histories of the wave elevations at gauges G01, G76 and G151 for A0=0.04 m in Group A, C, E

3.2 港内最大波面

图4呈现了各组底墙处的最大爬高随入射波波幅的变化曲线。图中Am表示入射波在港内的最大爬高。为了更好地区分各组实验的最大爬高，通过A0对其进行了无因次化处理。由图4可以看出，入射N波波幅、港内地形的平均水深以及折点的数量都对最大爬高有一定的影响。对于各组而言，随着入射波波幅的增加最大爬高，均呈现先增大后缓慢减小的趋势。从纵向来看，可以得出港内平均水深越大，则最大爬高越大，C组与F组的平均水深相等，但C组的放大因子要略小于F组，所以折点的存在将会略微减小最大爬高。

3.3 港内波能

由于正交模态分解法只有在港内波浪非线性较弱的时候才能精确分离各共振模态的响应幅值[12]。故本文所分析的入射波幅为0.01～0.07 m，间隔为0.01 m。图5为波幅A0=0.02 m时B组和D组通过FUNWAVE-TVD模型模拟和使用正交模态分解法拟合得到的自由波面的对比。本文只研究了最低的40个共振模态。由图5不难发现两者吻合良好。B组实验中模拟和拟合的自由波面达到最大的时刻分别为235.6 s与235.1 s，最大爬高分别为0.083 3 m和0.083 5 m，最大下冲均为-0.093 9 m。D组实验中模拟和拟合的自由波面达到最大波面的时刻分别为255.6 s和254.7 s，最大爬高分别为0.085 1 m与0.084 8 m，最大下冲均为-0.093 7 m。可以看出，模拟和拟合的自由波面存在一定的误差。考虑到这些最大爬高和下冲的特殊性，将正交模态分解法的数值拟合误差NFE重新定义为模拟和拟合的最大爬高和下冲的最大相对误差，即：

NFE=

(11)

式中：Au指模拟的最大爬高;(Au)fitted指拟合的最大爬高；Ad指模拟的最大下冲;(Ad)fitted指拟合的最大下冲。NFE反映了不同共振模态响应幅值的精确性。B、D各组的NEF分别为0.24%、0.35%。图6呈现了所有实验中正交模态分解法的数值拟合误差。如图6所示，各组实验的NEF均随着入射波波幅的增大而增大，这是由于波幅的增加使港内共振的非线性变强，从而使数值拟合误差增大。这与Gao等[12]的研究结果一致。为保证不同共振模态响应波幅模拟结果的精确性，将各组实验的NFE控制在5%以内。

图4 各实验组次的港内最大爬高Fig.4 Maximum runups inside the harbor for all groups

图5 A0=0.02 m时模拟与拟合的自由波面的对比Fig.5 Comparison of simulated and fitted free wave surfaces for A0=0.02 m

本文将相对幅值定义为各共振模态响应幅值与相应的入射波波幅之比。图7显示了在入射波波幅分别为0.1、0.3、0.5 m条件下各组实验最低的40个共振模态的相对幅值分布。通过对图7的分析可以得出，当入射波幅相同时，港内地形的平均水深对港内波能的分布有一定的影响。总的来说，港内波能分布会随着平均水深的逐渐增加而向高模态偏移，港内波能分布更加均匀。以A0=0.05 m为例，A～E组港内水深逐渐增加。A组的港内的波能主要分布在前8个共振模态上，占据最高波能的为第3模态。E组的港内的波能主要分布在前10个共振模态上，占据最高波能的为第4模态。F组为直线型地形，其平均水深与C组相同，且与其他组相比不存在折点，不难发现F组中各模态的相对幅值分布较其他组要略小。

港内波能的分布不仅与港内平均水深有关，与入射N波波幅的大小也有一定的关系。图8呈现的是在A组、C组和E组中不同入射波幅对于港内响应幅值分布的影响，如图所示，当入射波波幅较小时，港内的波能主要分布在最低的几个共振模态上，只有很小的一部分波能分布在更高的共振模态上，且在较高的共振模态上，其波能的分布更加均匀。随着入射波幅的不断增加，波能的分布开始向较高的模态上偏移。即入射波幅越大，港内相对波能分布越均匀。以C组为例，当入射波幅较小时，港内的波能主要分布在前6个共振模态上，占据最高波能的共振模态为第2模态。

随着入射波幅的增加，港内的波能主要分布在前21个共振模态上，占据较高波能的共振模态为第4模态和第12模态。

图6 各组数值实验的拟合误差(NEF)Fig.6 Numerical fitting errors (NEF) for all cases

图7 A～F组中港内共振模态的相对幅值分布的比较Fig.7 Comparison of relative amplitude distributions of resonance modes in group A～F

图8 在A组、C组、E组中不同入射波幅对于港内响应幅值分布的影响Fig.8 The influence of different incident wave heights on the distribution of response amplitude in pot in group A, group C, group E

由于将港内自由波面视为各共振模态线性叠加，故港内总波能为：

(12)

式中Ai(i=1,2,…,40)表示的是第i个共振模态的响应幅值。图9显示了入射波波幅A0为0.01、0.03、0.05 m时各组实验港内总波能随着平均水深的变化情况及A～E组港内总波能的线性拟合曲线。为了更直观观察图形，本文对总波能进行了无因次化处理(除以入射波幅的平方)。从图中可以看出对于入射波幅为0.01、0.03、0.05 m时A～E组的折线型地形而言，港内总波能随着港内平均水深增加而线性增加。C、F这2组实验的港内平均水深相等，当入射波波幅较小的时候，平台型地形的港内总波能要小于直线型，但当入射波幅较大时，前者反而大于后者。