罗丽，柳淑学*，李金宣，王磊

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

1 引言

真实的海浪是多向不规则的波浪，其方向分布对波浪传播特性及波浪与结构的相互作用有着重要的影响。为了准确地研究海洋环境及海洋工程中的波浪特性，有必要将实测的波浪历时过程线（波浪时间序列）进行再次模拟和重现，即波浪的确定性模拟。

一些学者对二维不规则波浪的确定性数值或物理模拟取得了较好的结果。Baldock 和Swan[1]基于对实测波浪在空间和时间上进行傅里叶级数展开，在数值模型中求解了具有非线性特性的单向不规则波浪，成功的通过数值模拟方法重现了物理实验的波浪过程线。裴玉国等[2]利用一个基本波列和一个瞬态波列进行线性叠加的方式对实测的畸形波进行了定点模拟生成，模拟产生的波列和实测的波浪时间过程吻合良好。Liang 等[3]基于快速傅里叶变换（FFT）方法，将二维不规则波分解为具有不同初始相位的组成波，数值重现了不规则波过程。崔成等[4]基于单波列叠加模型，在完全非线性的数值波浪水槽对不规则的波列进行了确定性模拟。Finnegan 和Goggins[5]通过FFT 方法将大西洋上实测的一点海浪表示为线性波浪叠加的形式，将此叠加的波面作为二维计算流体动力学（CFD）数值模型的入射边界，并将数值模型输出的波浪过程线分别与滤波分析后的波浪过程线及实测的海浪过程线进行了对比分析，验证了重现波浪的准确性。沈王刚[6]基于非线性波浪分离方法，在二维势流数值波浪水池中确定性重构了具有反射成分的强非线性波浪的波动过程。

然而，上述研究大都针对的是单向不规则波浪，而真实的海浪是多向不规则的。但是关于多向不规则波浪的确定性模拟，目前国内外研究还较少。Naaijen 和Blondel-Couprie[7–8]基于3D FFT 方法对多向不规则海浪在较大区域进行了模拟及预测，但是，海浪的模拟精度随着时间的增长而降低。刘思[9]将多向不规则波群一个测点的波面信息进行FFT 变换，得到各组成波的频率、振幅和相位的信息，然后根据多向不规则波的理论模拟方法[10]确定组成波的方向角，再根据线性叠加理论确定多向不规则波群的数值入射边界条件，对实测波列在数值水池的指定测点处进行了确定性模拟，但其主要致力于波浪要素、波浪群性和方向分布的研究。Draycott 等[11]采用基于时间序列差相位分析（Phase Time Path Difference，PTPD）的方法在圆形水池内重构了多向不规则波浪的方向谱，并对比分析了模拟波列及实验波列，但其研究主要致力于波浪统计参数如波浪方向谱的验证。Luo 等[12]通过等分波浪方向分布的方法确定性模拟了多向不规则波浪，但其组成波方向角与特定组成频率随机对应，波浪模拟存在一定的误差。

本文采用线性单叠加模型，建立了多向不规则波浪的确定性模拟方法。通过将模拟的波浪结果与已知波浪过程线进行对比分析，验证本文建立的多向不规则波浪确定性模拟方法的有效性，并确定了模拟方法的应用条件和适用范围。

2 波浪确定性模拟方法的建立

采用单叠加方法对多向不规则波浪进行确定性模拟，可将波浪的波面η(x,y,t)表示为

式中，N为组成波个数；ai、fi和θi分别为组成波的振幅、频率及传播方向角；ki为组成波的波数；εi为组成波的初始相位，在[0, 2π]内均匀的随机分布。进行波浪确定性模拟时，N取2n−1，2n为对采集波浪进行FFT 时的数据采样长度。

公式(1)中，ai、θi及εi需要基于给定的波浪序列来确定。需要说明，为了按下述方法确定波浪的传播方向，需要3 个以上合适位置的同步时间序列。本文的浪高仪布置形式如图1 所示。

ai和εi的求解及选取参考Luo 等[12]的处理方法。ai取基于所有测点处的时间序列按单向波求得的组成波幅值aim的平均值，即：

式中，M+1 为阵列浪高仪的总个数；S(fim)为通过FFT分析方法[3]求解的第m个浪高仪处第i个组成波的谱密度值。初始相位εi近似取为图1 所示的O 点浪高仪通过FFT 变换后直接求解的波浪初相位。

图1 波浪浪高仪的布置形式Fig. 1 Sketch of the wave gauges

组成波方向角θi的求解采用Borgman[13]和Esteva[14]提出的PTPD 方法，即对3 个非共线的浪高仪同步测得的波浪过程线进行分析，求解组成波的传播方向角。其求解步骤如下：

（1）在M+1 个浪高仪中，根据的组合个数，确定三角形浪高仪对的组合总个数(本文取M=5)。

（2）对于组成波频率为fi的组成波，选取有效的三角形浪高仪对组合，即确保每个浪高仪对的任意两个浪高仪间距L12、L13和L23大于0.01Li且小于等于0.5Li[15–16]，其中，Li为频率fi基于波浪色散方程所计算的波长。

（3）对所有有效的三角形浪高仪对组合，通过互谱分析，计算相对相位Ф12(fi)和Ф13(fi)。

（4）确定每个有效的浪高仪对对应频率为fi的组成波的波浪传播方向角，即

式中，sgn(*)为符号函数；P的计算公式为

（5）最终，确定对应频率为fi的组成波的波浪传播方向θi。可以采用以下两种方法：在所有有效的浪高仪对中，对不同浪高仪组合所计算的方向角进行平均所得的θi（以下记为average-θi）；在所有有效的浪高仪对中，选取对应波浪能量最大的仪器对所计算的方向角θi（以下记为max-Energy-θi）。

3 确定性波列模拟结果分析

3.1 规则波确定性模拟

不规则波及多向波均是由单向规则波叠加得到的，因此单向规则波精确地确定性模拟是不规则波及多向波浪模拟的基础。作为示例，取水深h=0.5m，规则波传播方向θ=π/6，波高H0=0.04m，周期T0=250/256 s=0.977 s (f0=1/T0=1.024 Hz)，为了使得FFT 分析能够准确地分辨规则波浪的频率，取T0/Δt为整数，波浪采样间隔Δt=0.00977 s。计算波浪方向角的浪高仪布置见图1，其外接半径R/L=0.224 (R/L=C代表浪高仪间距为C倍的波长)，2n=4096。通过FFT 变换，可以求得：f0=1.024 Hz，H0=2×0.01994m，ε0=−0.012578 rad，average-θi和max-Energy-θi计算的波浪传播方向角分别为：θ=π×30.173/180 rad 和θ=π×30.215/180 rad。可以看出，通过本文方法计算得到的波浪参数与初始给定的波浪参数基本一致。将模拟的波浪过程线按照求出的H0、ε0、θ计算波高，即

式中，k表示波浪的波数。

作为示例，图2 给出了确定性模拟的波浪过程线(θ=π×30.173/180 rad)和已知理论计算波浪过程线在t=0～10 s 时于2 号和5 号浪高仪位置的对比。可以看出，模拟的波浪过程线和已知波浪过程线几乎完全重合。

图2 2 号（a）和5 号（b）浪高仪在T0=250/256 s，H0=0.04m，θ=π/6 时的波浪过程线Fig. 2 Comparison of the wave surface elevation histories for T0= 250/256 s, H0=0.04m, θ=π/6 of No.2 (a) and No.5 (b) gauges

3.2 理论单向不规则波浪的确定性模拟

进一步，以Ts=1.5 s，Hs=0.04m，θ=0°的理论单向不规则波浪为例，验证本文方法对单向不规则波浪确定性模拟的有效性。波浪模拟时的频谱采用Goda 改进的JONSWAP 谱(谱峰升高因子取为3.3)[17]。图3给出了理论波浪频谱在fi∈(0，4fp)(fp表示波浪的谱峰频率)的分布。波浪模拟时，采用时间间隔Δt=0.01 s，采样长度2n=8192。图4给出了R/Ls=0.06 和0.12 时，基于average-θi和max-Energy-θi的组成波方向角误差θEr沿fi的分布情况（Ls表示有效波波长）。θEr的定义如下：

图3 理论波浪频谱(Ts=1.5 s，Hs=0.04m)Fig. 3 Theoretical wave frequency spectrum (Ts=1.5 s，Hs=0.04m)

图4 R/Ls 为0.06 和0.12 时，average-θi 和max-Energy-θi 的θEr 沿fi 的分布（Ts=1.5 s, Hs=0.04m）Fig. 4 Variation of the θEr for average-θi andmax-Energy-θi with fi for R/Ls=0.06 and R/Ls=0.12 (Ts=1.5 s, Hs=0.04m)

式中，θi_original和θi分别为已知的和计算的组成波方向角。通过图4 可以看出，R/Ls=0.06 和0.12 时的θEr分别在fi>3.2fp和fi>2.3fp时突然增大，这说明计算波浪方向角的fi超出上述范围后，计算的组成波方向角与理论组成波方向角不一致。这种现象产生的主要原因是由于计算θi的浪高仪间距R大于频率fi对应波长的0.5 倍所导致的，这与Fernandes 等[15–16]的结论是一致的。此外，通过对比average-θi和max-Energyθi的θEr沿fi的分布可以发现，整体而言，采用averageθi计算的组成波方向角的精度相对较高，θEr相对较小。

进一步，以average-θi在R/Ls=0.06 时计算的组成波传播方向模拟计算的波浪序列为例，图5 给出了t=20～40 s，波浪模拟空间范围的rr/Ls=0.12、0.20、0.35 时，已知波列和重现波列在5 号浪高仪位置的对比（rr代表指定位置与O 点的空间距离）。可以看出，确定性模拟的波列和原始波列几乎完全重合，仅有微小的振幅高度差及相位偏差，且随着rr/Ls的增大，波浪的模拟误差略有增大。

为了定量评估本文提出的波浪确定性模拟方法，定义误差系数Er为

式中，ηoriginal为已知波浪的过程线；ηdet为模拟所得波浪的过程线；ηmax//original为已知波列的最大幅值。以R/Ls为0.06 和0.12 时确定的θi进行波浪模拟的结果为例，图6a 和图6b 分别给出了基于average-θi和max-Energy-θi进行波浪模拟时的误差Ermax随rr/Ls的变化，其中Ermax为rr/Ls确定后，1～5 号浪高仪位置的波浪确定性模拟的最大误差。可以看出，Ermax随着rr/Ls的增大而增大。当R/Ls=0.06 时，基于average-θi和max-Energy-θi的波浪确定性模拟误差基本一致，但当R/Ls=0.12 时，采用average-θi进行波浪模拟时的误差明显比采用max-Energy-θi时的小。此外，结合图3 及图4 进行分析，可以看出，fi∈(0，3fp～4fp)时各fi对应的θEr大小是影响波浪确定性模拟精度的重要因素。

3.3 理论多向不规则波浪的确定性模拟

理论多向不规则波浪的模拟基于孙忠滨等[18]采用的频率−方向对应模型。波浪模拟的频谱同样采用Goda[17]改进的JONSWAP 谱(谱峰升高因子取为3.3)。波浪方向分布采用Mstuyasu 型方向分布函数[17]，即

式中，s为方向分布集中度参数；θ0为波浪的主波向。

图5 单向不规则波已知波列和模拟波列在不同rr/Ls 时的对比（Ts=1.5 s，Hs=0.04m）Fig. 5 Comparison of the original and simulated uni-directional irregular wave surface elevation histories (Ts=1.5 s，Hs=0.04m)

图6 采用average-θi 和max-Energy-θi 模拟所得波浪的Ermax 随rr/Ls 的变化（Ts=1.5 s，Hs=0.04m）Fig. 6 Variation of Ermax for simulated waves with rr/Ls using average-θi andmax-Energy-θi respectively (Ts=1.5 s, Hs=0.04m)

理论多向不规则波的波浪参数为：Ts=1.5 s，Hs=0.04m，θ0=0°，s=25。波浪模拟的采样时间间隔Δt=0.01 s，模拟的序列长度2n=16384。

与单向不规则波浪的研究类似，图7 以averageθi计算的θEr为例，给出了R/Ls不同时θEr沿fi∈(0，4fp)范围内的分布。明显地，θEr沿fi的误差分布范围随着R/Ls的增大而增大，且θEr主要表现在高频部分。这说明，R/Ls越小，计算的组成波方向角θi在低频部分的误差越小。

图7 R/Ls 不同时，采用average-θi 方法计算的θEr 沿fi 的分布（Ts=1.5 s，Hs=0.04m，s=25）Fig. 7 Variation of θEr with fi for average-θimethod with different R/Ls (Ts=1.5 s, Hs=0.04m, s=25)

图8 rr/Ls=0.35，R/Ls 不同时，3 号浪高仪位置理论波列和确定性模拟波列的对比（Ts=1.5 s，Hs=0.04m，s=25)Fig. 8 Comparison of theoretical and simulated wave surface elevation histories for different R/Ls with rr/Ls=0.35 at No.3 gauge (Ts=1.5 s,Hs=0.04m, s=25)

为了进一步验证R/Ls的变化对多向不规则波浪确定性模拟的影响，图8 以采用average-θi所得波浪模拟结果为例，给出了rr/Ls=0.35，不同R/Ls时，3 号浪高仪位置的理论波列和模拟波列的对比。可以看出，当R/Ls≤0.12 时，模拟波列和理论波列的波动过程是基本一致的，但当R/Ls=0.16 时，两者的差别相对较大。结合图7 的θEr在fi∈(0，4fp)范围内的分析结果可认为，模拟波列同理论波列差别产生的原因是由于当R/Ls=0.12，fi<2.3fp时，绝大多数组成波方向角误差θEr<2%，仅有个别的组成波方向角误差θEr稍大，但也小于10%，保证了绝大多数的波浪能量对应的计算的组成波的传播方向θi是准确的，从而保证了波浪模拟的整体精度。明显地，当R/Ls=0.16 时，在fi<1.95fp内的θEr大于R/Ls>0.12 时的θEr，且绝大多数θEr≈5%；此外，波浪高频成分组成波的θEr也会导致重现波浪在高频部分产生较大的模拟误差，从而导致模拟波列和已知波列整体差别较大。

为了研究多向不规则波确定性模拟的有效范围，图9 以R/Ls=0.06，选取average-θi进行的波浪确定性模拟结果为例，给出了3 号浪高仪位置处理论波列和模拟波列在t=20～40 s，rr/Ls=0.16、0.35、0.60 和0.80时的对比。可以看出，随着波浪模拟空间范围的增大，确定性模拟的波浪精度减小，模拟的波列与理论波浪过程线在时域上存在的振幅差及相位偏差随着rr/Ls的增大而逐渐增大。当rr/Ls=0.80，即确定性模拟的波浪范围为1.6Ls时，波浪确定性模拟的误差相对较大，但模拟波列和理论波列仍然具有一致性。

与单向不规则波浪的波浪模拟精度分析类似，图10以R/Ls=0.06 为例，给出了多向不规则波采用averageθi和max-Energy-θi进 行 波 浪 模 拟 时 分 析 的Ermax随rr/Ls的变化。明显地，当rr/Ls确定后，采用averageθi方法计算的组成波方向角进行波浪确定性模拟时的Ermax较采用max-Energy-θi方法计算的组成波方向角进行波浪模拟时的Ermax小。这与单向不规则波浪所得的结论是一致的。

图10 R/Ls=0.06 时采用计算的average-θi 和max-Energy-θi 模拟所得波浪的Ermax 随rr/Ls 的变化（Ts=1.5 s,Hs=0.04m, s=25）Fig. 10 Variation of Ermax for simulated waves with rr/Ls using average-θi andmax-Energy-θi of R/Ls=0.06 (Ts=1.5 s,Hs=0.04m, s=25)

图11 不同的R/Ls 采用average-θi 模拟所得波浪的Ermax 随rr/Ls 的变化(Ts=1.5 s, Hs=0.04m, s=25)Fig. 11 Variation of Ermax for simulated waves with rr/Ls using average-θi with different R/Ls (Ts=1.5 s, Hs=0.04m, s=25)

图11 以采用average-θi计算的组成波方向角进行确定性模拟的多向不规则波为例，给出了R/Ls不同时，波浪模拟的Ermax随空间范围rr/Ls的变化情况。随着rr/Ls的增大，波浪模拟的Ermax增大；当rr/Ls相同时，Ermax随着R/Ls的增大而增大。明显地，rr/Ls一定时，R/Ls=0.16 时的Ermax与R/Ls=0.12 时的Ermax绝对差值|Ermax_0.16−Ermax_0.12|明显地比|Ermax_0.12−Ermax_0.09|、|Ermax_0.09−Ermax_0.06|及|Ermax_0.06−Ermax_0.03|时 的 大；且R/Ls=0.03～0.12时，确定性模拟波浪的Ermax的差别不是很大。因此，为了保证采用本文方法进行波浪确定性模拟的准确性，采集波浪的浪高仪间距应小于波浪有效波长的0.12 倍。另外，综合图11 和图9b 的结果可以看出，当R/Ls≤0.12 时，在空间范围rr/Ls<0.35 内，即确定性模拟的波浪范围小于0.7Ls时，模拟的波列与原始波列的差别较小，即Ermax≤20%，一般满足波浪对结构物作用的试验要求。

4 结论

基于单叠加模型，通过对给定的合适位置的多向不规则波浪同步波动过程求解，确定各组成的波浪成分ai、ωi、θi及εi，并以单叠加数值计算模型，实现了多向不规则波的确定性模拟。通过定性对比分析模拟波列和已知波列在时域的分布情况及定量分析模拟波列的误差分布情况，可得出如下结论：

（1）本文建立的方法可以在一定范围内确定性的模拟给定的多向不规则波波动过程。但是，波浪的相位差和振幅高度差随着波浪模拟范围的增大而增大，即波浪重现误差随着波浪重现范围的增大而增大。

（2）采用对不同浪高仪组合所计算的方向角进行平均所得的average-θi为组成波方向角进行波浪确定性模拟的结果优于采用对应波浪能量最大仪器对所计算的方向角max-Energy-θi为组成波方向角时的模拟结果。

（3）给定计算组成波传播方向θi的各浪高仪间距应小于波浪有效波长的0.12 倍。

（4）当R/Ls≤0.12 时，在空间范围rr/Ls<0.35 内，即确定性模拟的波浪范围小于0.7Ls时，模拟的波浪与原始波浪的差别较小，即Ermax≤20%，一般满足波浪对于结构物作用的试验要求。