港湾振荡是指当近岸海域的低频波浪传播到港湾内，如果其频率和港湾的本征频率相一致，那么就会激起港湾内剧烈的水体共振的现象。港湾振荡可以诱发港内泊船的剧烈运动，显著影响系泊安全和货物的装卸^[1]。海啸波已被证实能够诱发显著的港湾振荡^[2]，且由其诱发的瞬变港湾振荡往往具有破坏性^[3]。

以往对于港湾振荡的研究主要集中于由平稳的波浪所诱发的稳态振荡^[4-5]，瞬变港湾振荡的研究工作起步较晚且研究的也较少。使用理论和实验方法，Lepelletier于20世纪80年代开始了对海啸等非线性瞬变长波诱发的瞬变共振的研究^[6]。Dong等使用物理模型实验并结合数值模拟对由滑坡产生的冲击波在矩形港内的共振反应进行了研究^[7]。Wang等基于Boussinesq模型对一个具有常底坡的矩形港内由水底运动诱发的瞬变港湾振荡进行了系统地研究^[8]。近期，作者使用Boussinesq模型模拟了孤立波入射诱发的狭长型港湾内的瞬变振荡现象，并采用正交模态分解法系统研究了港口底坡和入射孤立波波高的变化对港内相对波能密度分布的影响^[9]。

目前，海啸诱发的瞬变共振研究多限于使用孤立波来进行研究^{[7, 9]}，但从很多地震海啸的实际观测中发现多数海啸波是由一个大波峰和一个大波谷组成，形状像“N”，这类海啸波被称为N波，且根据波形的不同，N波可进一步分为波峰在前的N波、波谷在前的N波、等边N波和一般化的N波^[10-11]。

本文使用FUNWAVE-TVD模型和正交模态分解法进行了数值实验。通过对波峰在前和波谷在前的两种类型的等边N波诱发的狭长矩形港内的瞬变港湾振荡现象进行的数值模拟，验证了正交模态分解法的准确性。研究了入射N波形式和入射波波幅对港内相对波能分布的影响。

1 数值模型和分析方法 1.1 FUNWAVE-TVD模型

Boussinesq模型近些年来已被广泛地应用于模拟波浪在近岸的传播过程。它可以精确地模拟波浪的色散性、非线性、浅水变形、破碎和波成流等一系列的现象^[12]。本文使用FUNWAVE-TVD模型^[13]来实施数值实验。该数值模型是以Chen^[14]的完全非线性Boussineq方程为基础开发的。Kennedy等的动参考层技术被进一步应用于控制方程中，使得模型能够更好地模拟冲流带和沿海洪水的动边界现象^[15]。该模型使用一个非常均衡的守恒形式来组织控制方程，并使用一个高阶激波捕捉全变差下降格式来对控制方程进行数值离散。该技术允许模型能够不依赖于经验公式而很好地模拟波浪破碎和相关的波能耗散现象。此外，该模型的时间离散使用保持强稳定性的三阶Runge-Kutta方法，使得模型具备了自适应时间步长的功能。此外，模型使用了具有非阻隔通信的消息传递接口技术，使其可进行并行计算，显著提高了数值计算的效率。FUNWAVE-TVD模型已被证实能够很好模拟包括波浪浅水变形、反射、衍射、破碎、在岸滩上的爬高、下冲和波成流等一系列的近岸波浪水动力现象^[13]。

1.2 改进的模态分析法

正交模态分解法由Sobey^[16]首先提出，该方法可用于计算港口本征频率、共振模态形状和定量分离由风暴潮和海啸诱发的瞬变港湾振荡的各共振模态的响应幅值。该方法包含两个计算步骤：1) 以线性长波方程为起点，把控制方程转化为一个Sturm-Liouville问题，并最终离散为矩阵方程，矩阵方程的特征值和特征向量即为港口的本征频率和对应的共振模态形状；2) 通过一个多维优化问题对由风暴潮和海啸诱发的瞬变港湾振荡的各共振模态的响应幅值进行分离，第一步中计算得到的本征频率和共振模态形状将作为已知参数用于第二步的计算中。近期，Gao等^[17]使用镜像原理对正交模态分解法进行了改进，显著提高了其预测港口本征频率和共振模态形状的精度，并用此方法研究了孤立波入射条件下港内发生瞬变共振时港口底坡和入射孤立波波高的变化对港内相对波能密度分布的影响^[9]。本文首次将文献[17]中改进后的正交模态分解法用于N波诱发的瞬变港湾共振问题的研究，定量分离了港内各共振模态的响应幅值，并进一步系统研究了入射N波类型和入射波波幅的变化对于港内相对波能分布的影响。

2 数值实验布置

图 1展示了本文所用港口的平面布置情况。港口长度为L=1 000 m，宽度为W=20 m。港内外水深均为h=10 m。海岸线和港内边墙均设置为全反射直墙。笛卡尔坐标系统(o, x, y, z)的原点取在静水面，z轴垂直向上为正。

本文中，外海入射波采用Tadepalli和Synolakis^[10-11]提出的波峰在前和波谷在前的两种类型的等边N波，它们的表达式为

$\begin{array}{l} \eta \left( {x,0} \right) = \pm \frac{{\sqrt[3]{3}}}{2}{A_0}{\rm{sec}}{{\rm{h}}^2}\left( {k\left( {x - {x_0}} \right)} \right).\\ \quad \quad \quad \quad \tanh \left( {k\left( {x - {x_0}} \right)} \right) \end{array}$

(1)

其中

$k = \frac{3}{{2h}}\sqrt {\frac{{{A_0}}}{h}\sqrt {\frac{3}{4}} } $

(2)

式中：η表示自由波面，A₀表示入射波波幅，x₀表示初始入射N波的中心位置，方程等号右端的正号对应于波峰在前的等边N波，负号对应于波谷在前的等边N波。位于控制方程速度参考面处的初始水质点速度可用如下的线性表达式表示：

$u\left( {x,0} \right) = \sqrt {\frac{g}{h}} \eta \left( {x,0} \right)$

(3)

为了研究不同的入射N波类型对于瞬变港湾振荡的影响，本文设计了两组系列实验，即A组和B组。在A组中，外海初始入射波均为波峰在前的等边N波；而在B组中，外海初始入射波均采用波谷在前的等边N波。同时，为了研究入射N波波幅对于瞬变港湾振荡的影响，在这两组系列实验中，入射波波幅A₀均从0.005 m逐渐增大到0.30 m。所有数值实验的波浪参数见表 1。

图 2显示了入射波波幅A₀为0.005 m和0.30 m时波峰在前和波谷在前的等边N波的波形图。如图所示，当入射波幅相同时，波峰在前和波谷在前的等边N波波形的相位相差了180°。图中的圆点表示对应的N波的“波前”。波前被定义为其所在处的自由波面η_f=±0.05A₀。本文中，等边N波的波长L₀被定义为波前和中心位置x₀之间的两倍距离。当A₀=0.005 m时，波长L₀=1 703.6 m；而当A₀=0.30 m时，波长显著减小为L₀=220.0 m。

图 3为数值实验中所用的数值波浪水槽。数值波浪水槽尺寸为4 000 m×50 m。

因为图 1所示港口关于x轴对称，为了节省计算网格和时间，仅模拟港口和外海的一半区域。港内等间距布置了101个测点(G₀₁~G₁₀₁)，相邻两个测点的间距为10 m。测点G₀₁布置在港口口门处(x=0)，而测点G₁₀₁布置在港口底墙处(x=1 000 m)；x和y方向网格尺寸均为Δx=Δy=1 m。网格节点和单元数分别为204 051和200 000个；所有边界均设置为全反射；作为初始条件，入射N波的波前均设置在港口门位置处；所有的数值实验均模拟到300 s。

3 模拟结果和港内波能分析

在正交模态分解法的使用过程中，需要预先指定所需考虑的共振模态数量M^[16-17]。本文对于所有的数值实验均设置M=40。图 4显示了在实验A02和B02中使用FUNWAVE-TVD模型模拟得到的自由波面和对应的使用正交模态分解法拟合得到的自由波面的比较。图中，A_us和A_uf分别表示模拟和拟合得到的入射N波在港内的最大爬高值，而A_ds和A_df分别表示模拟和拟合得到的入射N波在港内的最大下冲值。对这两个算例，拟合的自由波面和模拟的自由波面吻合得很好。对于实验A02，入射波波幅为A₀=0.01 m，模拟的自由波面(图 4(a))在位于x=1 000 m的港口底墙处、在t=147.6 s的时刻出现了最大爬高值A_us=0.032 22 m；而在相同的位置和时刻，拟合的自由波面(图 4(b))出现A_uf=0.032 16 m的最大爬高值。类似地，模拟的自由波面在港口底墙处、t=178.7 s时出现了最大下冲值A_ds=0.033 19 m；而在相同的位置和时刻，拟合的自由波面也出现A_df=0.033 19 m的最大下冲值。对于实验B02，入射波波幅也为A₀=0.01 m，模拟的自由波面(图 4(c))在港口底墙处、t=148.0 s时出现了最大下冲值A_ds=0.032 21 m，在t=177.8 s时出现了最大爬高值A_us=0.032 24 m；而在相同的位置和对应的时刻，拟合的自由波面(图 4(d))分别出现了A_df=0.032 17 m和A_uf=0.033 20 m的最大下冲值和最大爬高值。本文定义了正交模态分解法在各组数值实验中用于拟合由数值模拟得到的自由波面的数值拟合误差(NFE)：

${\rm{NFE}} = \max \left\{ {\left| {\frac{{{A_{{\rm{uf}}}} - {A_{{\rm{us}}}}}}{{{A_{{\rm{us}}}}}}} \right|,\left| {\frac{{{A_{{\rm{df}}}} - {A_{{\rm{ds}}}}}}{{{A_{{\rm{ds}}}}}}} \right|} \right\} \times 100\% $

(4)

该参数可以直观地反映正交模态分解法分离不同共振模态响应幅值的精确性。对于实验A02和B02，正交模态分解法的数值拟合误差分别为0.19%和0.12%。表 2中呈现了所有数值实验中正交模态分解法的数值拟合误差。可以发现，所有实验的数值拟合误差均小于4.40%。这表明在所有的数值实验中，正交模态分解法分离得到的各共振模态的响应幅值是精确可靠的。

本文将各组数值实验中港内的各共振模态响应幅值与相应的入射波波幅的比值定义为各模态的相对幅值，即：

${{\bar A}_i} = \frac{{{A_i}}}{{{A_0}}}\left( {i = 1,2, \cdots ,40} \right)$

(5)

式中：A_i和A_i分别表示第i模态的幅值与相对幅值。

图 5显示了所有的数值实验中最低的40个共振模态的相对幅值分布。可以看出，相对幅值在各模态上的分布受入射N波波幅的影响很大。以A组为例(图 5(a))。实验A01中，入射波波幅为A₀=0.005 m，港内的波能主要集中在最低的四个模态上。随着入射波波幅的增大，分布于更高模态上的波能的比重逐渐增加，港内的波能分布逐渐趋于均匀。在实验A08中，入射波波幅为A₀=0.30 m，港内大部分的波能分布于更多的共振模态上(主要在最低的15个模态上)，港内波能分布趋于均匀。此外，在实验A01中占据最大相对波能的模态为第二模态，而随着入射波波幅的增大，占据最大相对波能的共振模态逐渐向较高的模态转移，在实验A08中具有最大相对波能的模态已转移到第七模态。对于B组(图 5(b))，同样能够观察到类似的现象。随着入射波波幅的增加，港内波能向高模态转移且波能分布逐渐趋于均匀的原因可能与入射N波波长的减小有关。

为了比较不同入射N波类型对港内相对波能分布的影响，图 6呈现了入射波波幅分别为0.005、0.15和0.30 m时波峰在前和波谷在前的两类等边N波诱发瞬变港湾振荡时港内相对幅值分布的比较。当A₀=0.005 m时(图 6(a))，两类不同的N波产生的港内相对幅值分布几乎完全相同。第二模态具有最大的相对幅值(相对波能)，随着共振模态的增大，相对幅值迅速地单调减小。当A₀增大到0.15 m时(图 6(b))，两类不同的N波产生的港内相对幅值分布出现了明显的不同。对于实验A05，港内最大相对幅值位于第六模态，为0.377；而在实验B05中，最大相对幅值位于第八模态，其值明显低于前者，仅为0.290，并且，总体来说，波能在港内的分布要比前者更为均匀。此外，还可以发现，在实验A05中，相对幅值分布除在第六模态存在最大值，在第十八模态还存在一个明显的次峰，而此现象在实验B05中并未出现。当A₀进一步增大到0.30 m时(图 6(c))，两类不同的N波产生的相对幅值分布的差别同样非常明显。在实验B08中，最大相对幅值仍然明显低于实验A08中的值，波能在港内的分布仍较实验A08中更为均匀。实验A08中相对幅值的次峰相比于实验A05中更加明显，而且实验B08中也出现了较为明显的次峰。

图 7显示了在实验A01和B01中测点G₁₀₁测得的自由波面时间序列和对应的小波谱。小波谱定性地表明了相对波能密度在时域和频域的分布。这两组实验中，虽然测点G₁₀₁测得的自由波面时间序列相位相差了180°，但是它们对应的小波谱几乎是完全相同的，即表明港内相对波能分布情况几乎完全相同。小波谱表明第二模态占有最大的波能。尽管第一模态占有的波能也很显著，但是要明显低于第二模态。第三模态占有的波能要比前两个模态的波能小得多。这些现象和图 6(a)中呈现的结果是一致的。

图 8显示了在实验A08和B08中测点G₁₀₁测得的自由波面时间序列和对应的小波谱。这两组实验中，测点G₁₀₁处的自由波面时间序列和对应的小波谱均存在明显不同。在实验A08中，波浪在港内的最大的爬高值为0.582 9 m，而实验B08中，波浪最大爬高值显著大于前者，为1.059 m。对于波浪在港内的最大下冲值，在实验A08和B08中，分别为-0.627 9 m和-0.488 6 m。在实验A08中，小波谱显示最主要的波能分布在频率f=0.03 Hz附近，另外在f=0.09 Hz附近的波能也较为明显，虽然相比于f=0.03 Hz附近的波能要小得多。而在实验B08中，小波谱显示主要的波能分布在f=0.06 Hz附近，且相比于实验A08来说，港内波能在频域的分布也较为均匀。这些现象和图 6(c)中呈现的结果是一致的。

为了进一步定量地表示港内相对波能在不同共振模态上分布的均匀程度，本文计算了各组实验中最低的40个共振模态的相对幅值分布的标准差系数：

${C_V} = \frac{\sigma }{\mu }$

(7)

其中：

$\begin{array}{l} \quad \quad \mu = \frac{1}{{40}}\sum\limits_{i = 1}^{40} {{{\bar A}_i}} \\ \sigma = \sqrt {\frac{1}{{40}}\sum\limits_{i = 1}^{40} {{{\left( {{{\bar A}_i} - \mu } \right)}^2}} } \end{array}$

标准差系数可以直观地反映各共振模态的相对幅值相对于它们的平均值的离散程度。显然，标准差系数越小，港内相对波能在各模态的分布越均匀；反之，标准差系数越大，港内相对波能在各模态的分布就越不均匀。图 9展示了所有数值实验中最低的40个共振模态的相对幅值分布的标准差系数C_V值。从图中可以发现，当入射波波幅较小时，C_V的值较大，并且A组和B组中相应的C_V值差异很小。这表明在此条件下港内相对波能在不同模态上分布的均匀性较差，但是不同的入射N波形式对港内相对波能分布的影响很有限。当入射波波幅逐渐增大时，A组和B组中的C_V值均单调减少，并且当入射波波幅一定时，B组中的C_V值要明显小于A组中相应的C_V值。这表明当入射波波幅逐渐增大时，港内的相对波能在各共振模态的分布逐渐趋于均匀，并且波谷在前的等边N波所产生的港内相对波能分布比波峰在前的等边N波所产生的港内相对波能分布要更为均匀。该图中所呈现的这些现象同图 5~8中所呈现的现象是一致的。

4 结论

1) 当入射波波幅较小时，波峰在前和波谷在前的等边N波诱发的港内相对波能分布几乎完全相同，并且波能几乎都集中在最低的几个共振模态上，分布于更高模态上的波能极其有限。

2) 随着入射波波幅的增大，分布于更高模态上的波能的比重增加，港内波能的分布趋于均匀，并且占有最大波能的共振模态逐渐由较低的模态向较高的模态转移。

3) 相比于波峰在前的等边N波，波谷在前的等边N波诱发的港内波能分布得更加均匀。最后，需要说明的是，以上结论仅针对本文所使用的港口、入射等边N波类型和入射波波幅范围。