港湾共振是指当传入港内的低频波浪的频率与港湾的本征频率相同时，港内水体会出现大幅波动的现象^[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(μ²)，在非线性a/h保留了所有阶(这里μ=h/l，h为局部静水深，l为波长，a为波幅)。其控制方程表达式为：

$ \eta_{t}+\nabla \cdot \boldsymbol{M}=0 $

(1)

$ \begin{array}{c} \boldsymbol{M}_{t}+\nabla \cdot\left[\frac{\boldsymbol{M} \boldsymbol{M}}{H}\right]+\nabla\left[\frac{1}{2} g\left({\eta}^{2}+2 h \eta\right)\right]= \\ H\left\{\overline{\boldsymbol{u}}_{2, t}+\boldsymbol{u}_{\alpha} \cdot \nabla \overline{\boldsymbol{u}}_{2}+\overline{\boldsymbol{u}}_{2} \cdot \nabla \boldsymbol{u}_{\alpha}-\boldsymbol{V}_{1}-\right. \\ \left.\boldsymbol{V}_{2}-\boldsymbol{V}_{3}-\boldsymbol{R}\right\}+g \eta \nabla h \end{array} $

(2)

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

$ \boldsymbol{M}=H\left\{\boldsymbol{u}_{\alpha}+\overline{\boldsymbol{u}}_{2}\right\} $

(3)

式中：u_α为参考面z_α=ζh+βη处ζ=-0.53、β=0.47时的速度；u₂为整个水平速度场平均水深的O(μ²)分布：

$ \begin{array}{c} \overline{\boldsymbol{u}}_{2}=\left[\frac{z_{\alpha}^{2}}{2}-\frac{1}{6}\left(h^{2}-h \eta+\eta^{2}\right)\right] \nabla B+ \\ {\left[z_{\alpha}+\frac{1}{2}(h-\eta)\right] \nabla A} \end{array} $

(4)

式中：A=▽·(hu_α)，B=▽·u_α。V₁和V₂为色散项。

$ \boldsymbol{V}_{1}=\left\{\frac{z_{\alpha}^{2}}{2} \nabla B+z_{\alpha} \nabla A\right\}_{t}-\nabla\left[\frac{\eta^{2}}{2} B_{t}+\eta A_{t}\right] $

(5)

$ \begin{array}{l} \boldsymbol{V}_{2}=\nabla\left\{\left(z_{\alpha}-\eta\right)\left(u_{\alpha} \cdot \nabla\right) A+\right. \\ \left.\frac{1}{2}\left(z_{\alpha}^{2}-\eta^{2}\right)\left(u_{\alpha} \cdot \nabla\right) B+\frac{1}{2}(A+\eta B)^{2}\right\} \end{array} $

(6)

V₃为竖直涡量场的O(μ²)分布：

$ \boldsymbol{V}_{3}=\omega_{0} \boldsymbol{i}^{\rm{z}} \times \overline{\boldsymbol{u}}_{2}+\omega_{2} \boldsymbol{i}^{\rm{z}} \times \boldsymbol{u}_{\alpha} $

(7)

其中：

$ \omega_{0}=\left(\nabla \times \boldsymbol{u}_{\alpha}\right) \cdot \boldsymbol{i}^{\rm{z}}=v_{\alpha, x}-u_{\alpha, y} $

(8)

$ \begin{array}{c} \omega_{2}=\left(\nabla \times \overline{\boldsymbol{u}}_{2}\right) \cdot \boldsymbol{i}^{\rm{z}}= \\ z_{\alpha, x}\left(A_{y}+z_{\alpha} B_{y}\right)-z_{\alpha, y}\left(A_{x}+z_{\alpha} B_{x}\right) \end{array} $

(9)

式中：i^z为竖直方向的单位向量；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个折点处的水深一致并且关于港口长度垂直平分线对称。图中h₁、h₀、h_p和l分别代表港口口门处水深、港口底墙处水深、平台处水深和平台的长度。在本文中，港口口门处的水深和外海水深均为12 m，底墙处的水深为3 m，平台处的水深为7.5 m；各组实验的底部剖面形状各不相同。港内的平均水深为：

$ \bar{h}=\left\{\begin{array}{l} \frac{3 h_{1}-2 h+h_{0}}{4}, & 1 \text { 个折点 } \\ \frac{3 h_{1}-2 h_{\mathrm{p}}+h_{0}}{4}, &2 \text { 个折点 } \end{array}\right. $

(10)

本文共模拟了6组数值实验，为A组~F组，A组与B组为凹型底，其折点处的水深h分别为10.5 m和9.0 m，港内平均水深分别为9 m和8.25 m。C组为平台型底，平台处水深与平均水深均为7.5 m，平台长度l为400 m。D组与E组为凸型底，折点处水深分别为6 m和4.5 m，平均水深分别为6.75 m和6.25 m。F组为对照组，地形为斜直线型地形，平均水深与平台型地形一样，为7.5 m。上述不同地形实验的入射波幅均一致，并以每次增加0.01 m的方式从0.01 m逐渐增加到0.07 m，共7组。图 2为数值实验所用的数值波浪水槽示意图，因为本文所使用的波浪水槽是关于沿港口长度方向的中心线对称的，为节省计算时间，本文只模拟了一半的区域。数值波浪水槽的尺寸为5 500 m×120 m，151个测点沿着港口长度方向均匀的分布在港口中心线上，相邻测点间的间隔为10 m。在口门处设置测点G₀₁，在港口底墙处设置测点G₁₅₁，本文所用的数值波浪水槽均为全反射边界。造波机位于港口口门1 500 m的外海中。每次数值实验模拟的初始条件均为将入射N波的波前置于港口口门处。整个计算域内网格尺寸Δx和Δy均设为1.0 m。网格节点和矩形单元总数的总数分别为335 561和330 000。总模拟时间为500 s。

3 结果和分析 3.1 波面时间序列

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

3.2 港内最大波面

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

3.3 港内波能

由于正交模态分解法只有在港内波浪非线性较弱的时候才能精确分离各共振模态的响应幅值^[12]。故本文所分析的入射波幅为0.01~0.07 m，间隔为0.01 m。图 5为波幅A₀=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重新定义为模拟和拟合的最大爬高和下冲的最大相对误差，即：

$ \begin{array}{c} \text { NFE = } \\ \max \left(\left|\frac{\left(A_{\mathrm{u}}\right)_{\text {fitted }}-A_{\mathrm{u}}}{A_{\mathrm{u}}}\right|, \left|\frac{\left(A_{\mathrm{d}}\right)_{\text {fitted }}-A_{\mathrm{d}}}{A_{\mathrm{d}}}\right|\right) \times 100 \% \end{array} $

(11)

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

本文将相对幅值定义为各共振模态响应幅值与相应的入射波波幅之比。图 7显示了在入射波波幅分别为0.1、0.3、0.5 m条件下各组实验最低的40个共振模态的相对幅值分布。通过对图 7的分析可以得出，当入射波幅相同时，港内地形的平均水深对港内波能的分布有一定的影响。总的来说，港内波能分布会随着平均水深的逐渐增加而向高模态偏移，港内波能分布更加均匀。以A₀=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模态。

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

$ E = \sum\limits_{i = 1}^{40} {\frac{1}{2}} A_i^2 $

(12)

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

4 结论

1) 在本文所研究的地形及入射N波波幅范围内，最大爬高随着港内平均水深的增大而减小，随着入射N波波幅的增大而增大，而折点的存在会略微降低最大爬高。

2) 在本文研究的地形以及港内平均水深范围内来说，港内总波能随着平均水深的增加而呈线性增加。

3) 当入射波为N波时，在港内水深相同的件下平台型地形与直线型地形的港内总波能的相对大小与入射波幅有关。当入射波波幅较小的时候，平台型地形的港内总波能要小于直线型，但当入射波幅较大时，前者反而大于后者。

本文的结论仅限于本文所考虑的细长形的港口形状，港内地形和入射N波的波幅变化范围。