基于条件模拟理论的波浪数值模拟

张 梦，柳淑学，李金宣，张昊宸

(大连理工大学 海岸和近海工程国家重点实验室，大连 116024)

条件模拟是地质统计学中的一个概念。在地质统计学中，因研究对象是一种具有结构特征的随机变量，对它进行数学模拟时要求保持一定的空间相关性，即保持实验信息具有相同的协方差函数或变差函数，在地质统计学中称非条件模拟(unconditional simulation)；如果再增加一个更严格的约束条件，将模拟条件化，即令在各观测点处的模拟值等于该点的实测值，称其为条件模拟(conditional simulation)[1]。

近年来，由于一些海上工程规模的庞大和工程问题的复杂性，其对波浪荷载精确程度的要求也越来越高。目前在实验室中进行物理模拟主要是通过模拟频谱，不能模拟波形。很多时候研究人员常常需要研究某一特定短波列或波群对建筑物产生的作用，且这段短波列包含于满足某一特定波谱的随机波列中，这往往需要进行逐波分析以找到满足此条件的随机波列，费时费力。因此如何有效地产生这样的特殊波列对研究波浪对于复杂工程作用问题具有重要的意义。

起初，波浪模拟主要都是通过模拟频谱的方式来实现的。Borgman[2]建立了线性波浪叠加法和线性过滤法，过滤法主要利用白噪声滤波的方式对波浪进行模拟，线性叠加法中初相位随机分布，因此也称为随机相位谱法。后来，人们发现有时是某种特殊的波列对建筑物产生主要的破坏作用[3]，因此需要模拟一些特殊波列。刘思等基于改进的波包谱[4]的经验公式，建议了单向不规则波群的数值模拟方法[5]，但只是模拟了同时满足频谱及群性要求的波浪，并不是严格地模拟了不规则波列。Borgman[6]提出波浪的条件模拟概念，即在满足波浪的统计特征情况下，人为加入一些约束，阐述了如何利用正态分布的条件概率理论，在时域上产生一段既包含特殊短波列，又属于某一波谱的目标波列。但当计算波列的长度增大时，此方法暴露出一个很大的缺点：自相关函数矩阵越来越病态导致计算结果出现严重偏差。Spanos[7]建立的ARMA波浪模拟模型在计算时遇到同样的病态Toeplitz矩阵。Medina和Sanchez-Carratala[8]对数值波浪模拟中的许多计算模型进行了比较，通过在波浪中加入白噪声的方式很好地解决了这个问题。Gime′nez等[9]从上述方法得到启示，并将其运用到条件波浪模拟中，解决了算法中出现的病态矩阵问题。Hudspeth等[10]在二维水槽中对目标波列进行了重现，同时从能量的角度讨论了嵌入波列的参数对模拟产生的影响。

天然海浪波列中存在的极限波浪，是造成海洋工程建筑物破坏的主要因素之一，而波浪聚焦是目前实验室产生极限波浪的主要方法，即让不同方向、不同频率的波浪在传播过程中互相调制，使波能发生集中，形成极端大波。但是这样模拟波浪的过程中畸形波的产生具有随机性，因此如何实现在不规则波的波列中嵌入已知的畸形波过程，在特定地点和时间产生包含已知畸形波列的确定性模拟，亦是进一步提高物理模型试验水平的需要。Oggiano等[11]基于数值二维水槽，针对一极限波浪过程，建立可以复演该过程的造波边界条件，但其只是针对孤立的极限波浪过程，并非针对整个不规则波过程。Buldakov等[12]则对于给定的二维单向特定的极限波浪及其破碎型式进行了确定性的模拟研究，但是其针对的是短序列的波群序列，而非一完整的包含畸形波的不规则波序列。Pierella等[13]对于二维极限波浪，基于流函数理论，将其嵌入一不规则波波列，实现其在物理水槽中的确定性的模拟。

本文基于Borgman的条件模拟波浪理论，通过加入白噪声的方式，解决了计算中因病态矩阵导致的数值不稳定问题，进而建立了二维数值波浪模拟程序，可在满足某一特定波谱的随机波列中，在任意合理时间点，嵌入某一特定波列后，产生的新的随机波列仍然满足上述特定波谱。分别以规则波和聚焦波为嵌入波列，研究了嵌入波的长度、周期、波高和嵌入点等对模拟波浪产生的影响。

1 基于条件模拟的波浪模拟方法

所谓条件波浪模拟，首先通过基于给定波浪参数和波浪谱，模拟生成一初始的随机波列ηu(t)，可称之为非条件模拟；之后在指定的时间点嵌入一段特殊的波列ηe(t)，从而形成新的同时仍满足给定的波浪谱和统计特征参数的随机波列ηc(t)。

对于初始的随机波列ηu(t)，可采用非条件波浪模拟得到。即假定海浪可看作一平稳随机过程，由无数多个不同的余弦波随机叠加而成[3]

(1)

(2)

把代表M个区间内波能的各余弦波动叠加起来，即得不规则波的波面

(3)

为了在初始波列中某一时间点嵌入某一特定的波列，可将初始波列的波面ηu(t)前后分为3段，表示为{ηu-,vu,ηu+}T，其中vu代表将要被嵌入的短波列替代的部分，而目标波列可表示为{ηc-,vu,ηc+}T；若ve表示将要嵌入的波列，即要求

vc=ve

(4)

接下来就是如何求解ηc-、ηc+，使得所形成的新的波列既满足所要求的波浪谱，同时嵌入已知的某特定波列。这里假设波浪的波面满足正态分布，则由多向正态分布的条件概率公式可得

(5)

式中：C12T、C11-1为初始波列自相关函数矩阵的分块矩阵，其表达式可参见Hudspeth等[10]，详细证明可见Borgman[6]。

2 数值模拟及结果分析

2.1 嵌入规则波数值模拟结果分析

表1 初始波列参数Tab.1 Unconditional simulation wave parameters

表2 嵌入规则波参数Tab.2 Embedded regular wave parameters

(6)

式中：ηu为初始波列的过程线；ηc为目标波列的过程线；N为序列长度；Hmax为初始波列的最大波高。需要说明的是，式(6)中误差累加计算不包括嵌入波列的部分。

1-a 组次R1 1-b 组次R2

图2 初始波列的自相关函数曲线图 图3 波面误差随嵌入点的变化(组次R2)Fig.2 Auto-correlation of the unconditional simulation wave Fig.3 Variation of RMS error coefficient with the embedded time point

图4 目标波列与初始波列特征波高的比值随嵌入点的变化(组次R2) 图5 目标波列与初始波列特征周期的比值随嵌入点的变化(组次R2)Fig.4 Comparison of the variation of the ratios between the unconditional simulation and conditional simulation characteristic wave height with the embedded time point Fig.5 Comparison of the ratios between the variation of the unconditional simulation and conditional simulation characteristic wave period with the embedded time point

图6为各组次目标波列与初始波列的波谱比较图，表3为嵌入规则波后波列各统计参数计算结果。具体来说，从组次R1、R2的统计参数和谱对比来看，进一步说明嵌入点的改变对模拟产生的波浪分析结果影响不大；同样，从组次R1、R3的结果对比来看，当嵌入规则波列的长度增大了9倍，但其统计参数和谱的变化亦较小，因此可以认为当嵌入规则波列的周期和波高与初始波列的有效周期和波高相差不大时，嵌入规则波列的长度变化不会对模拟波浪产生明显的影响；从组次R1、R4的结果对比来看，嵌入规则波列的波高增大了1倍，而其统计有效波高参数和谱亦没有产生很大的变化，但是如果嵌入的规则波波高较大或大于初始波列的最大波高时，新的波列的最大波高会大于初始波列的最大波高。

而从组次R1、R5的结果对比来看，嵌入波列的周期会对模拟得到的新的波列产生较大的影响。从图1组次R5可以明显看出，目标波列与初始波列相比，其在嵌入波列周围产生了较大的特异波浪，从表3中也可以看到各统计波高明显增大，从图6组次R1、R5的对比可以看到，嵌入规则波后的波浪谱整体明显变大。也就是说，嵌入规则波的周期与初始波列谱峰周期偏离越大，嵌入特定规则波后波浪的特征与初始波列特征相差就越大。

6-a 组次R1、R2 6-b 组次R1、R3

表3 嵌入规则波时初始波列和目标波列的统计参数表Tab.3 Statistical parameters of unconditional simulation and conditional simulation after embedding regular waves

综合上述讨论可以看出，当嵌入规则波列时，嵌入点、嵌入波列的长度和波高变化对初始波列的统计特征参数和谱的影响不大，可以认为条件模拟是有效的，但是当嵌入规则波的波高较大时，嵌入规则波后的波列的最大波高会增大。事实上，从某种程度上看，嵌入波列的周期变化也没有对统计参数和谱产生较大的影响，但是由于嵌入波列的周期偏离了初始波列的有效周期，统计特性的不同使得初始波列产生类似“排斥”效应，从而在嵌入波浪附近产生了特异波浪。这种含有远远大于随机波列有效波高的特异波的长波列不满足常规波浪特征，因此嵌入波列的周期应与初始波列的周期基本一致。

2.2 病态矩阵问题

在上述计算过程中，有时会出现如图7所示情况。图7中7-a、7-b、7-c分别代表嵌入一个波高为0.1 m，周期分别为2.16 s、2.18 s、2.2 s，长度为一个周期的规则波。在周期改变很小的情况下，图7-b的模拟结果却出现了很大的波动。其原因是因为计算公式(5)中的矩阵C11有时为病态的，使得计算结果不稳定。Spanos[7]通过ARMA算法进行波浪模拟时遇到了相似的问题。根据Spanos和Mignolet[14]的研究，发现其原因是谱在低频率时波浪谱密度值为0引起的。根据此特点，Medina和Sanchez-Carratala[8]通过在波浪谱中人为地加入白噪声，从而很好地解决了这个问题。假定加入的白噪声占总能量的ε，则新的谱为

7-a Te=2.16 s 7-b Te=2.18 s

(7)

式中：fL、fH分别表示频率的上下界，其取决于模拟波浪时在低频和高频处略去的波浪的能量占波浪总能量的比例[3]，本文取1%。虽然此种方法改变了原有谱形，但是Gime′nez等[9]认为无论是实验室产生的波浪还是海面上记录的波浪都存在着噪声，因此当ε取得足够小时，S*可近似等于S。图7-d为取ε=0.000 5所计算得到的新结果，与图7-b相比有了很大改善。本文所涉及到的谱的计算全部采用S*。

2.3 嵌入聚焦波数值模拟结果分析

2.3.1 二维聚焦波的产生方法(相速度法)[15]

聚焦波可通过模拟频谱的方法产生。对于二维波浪，任一点处波浪自由表面可以表示为不同频率和不同幅值的余弦波线性叠加后的结果，即

(8)

式中：aj为第j个组成波的幅值；kj为第j个组成波的波数；ωj为波角频率；φj为j个组成波的初相位；N为组成波的总个数。kj和ωj满足波浪的色散关系

(9)

式中：g和d分别为重力加速度和水深。若假定波浪在指定时刻t=tb聚焦于x=xb位置处，即各组成波在x=xb处波峰叠加，那么需满足

cos(kjx-ωjt-φj)=1

(10)

则各组成波的初相位应满足下式

φj=kjxb-ωjtb

(11)

此时，波浪的波面表达式为

(12)

由式(12)可以看出，对于给定的聚焦时间tb，聚焦点的坐标xb，频率数N，影响波浪聚焦面参数主要为各组成波的幅值aj。通常把A定义为聚焦点处的波浪幅值，则有

(13)

则生成聚焦波浪组成波的幅值可以表示为

(14)

式中：频谱S(f)为本文采用的JONSWAP谱。

2.3.2 数值模拟及结果分析

以聚焦点为中点向两边对称截取足够长度得到欲嵌入的聚焦波列，表4为嵌入的聚焦波参数。同理，分别研究嵌入点、嵌入聚焦波的历时长度、幅值和谱峰周期对模拟波浪产生的影响。

8-a 组次fc1 8-b 组次fc2

9-a 组次fc1、fc2 9-b 组次fc1、fc3

表4 嵌入聚焦波参数Tab.4 Embedded focusing wave parameters

图8为嵌入聚焦波后各组次目标波列与初始波列的波面比较图，从结果来看，其主要特点与嵌入规则波时类似，即嵌入聚焦波列后目标波列波形只在嵌入波列周围发生变化。图9为嵌入聚焦波后各组次目标波列与初始波列的波谱比较图，表5为嵌入聚焦波后波列各统计参数变化表。从各组次对比结果来看，嵌入聚焦波时，其影响因素也与嵌入规则波时相似，即嵌入点、嵌入聚焦波的长度、幅值不会对模拟产生明显的影响，而嵌入聚焦波的谱峰周期会对模拟波浪产生较大的影响。

表5 嵌入聚焦波时初始波列和目标波列的统计参数表Tab.5 Statistical parameters of unconditional simulation and conditional simulation after embedding focusing waves

这里也可以看到与嵌入规则波不同的地方。从fc1和fc3组成结果的对比可以看到，聚焦波列的长度变大反而使得其统计参数和谱变化更小了，这是因为聚焦波列的长度越长，统计特性更加与初始波列趋于一致，初始波列产生的“排斥”因此就越小。从周期的影响上来看，聚焦波谱峰周期的变化比规则波周期的变化产生的影响更大，同样的嵌入聚焦波的周期为1.5倍的初始波列有效周期，嵌入聚焦波时产生的特异波是嵌入规则波时的1.7倍。

3 结论

本文基于Borgman的条件波浪模拟理论，建立了相应的数值波浪模拟程序，可实现模拟生成既满足某一特定波浪谱，又包含了某一特殊短波列的随机波浪，并通过加入白噪声的方式解决了计算中产生的数值不稳定问题。

嵌入波浪参数对于嵌入波浪后模拟波浪的影响研究结果表明，在进行条件波浪模拟时，嵌入波列的周期是主要影响因素，其与初始波列的有效周期越接近，则模拟效果越好，因此需要限定嵌入波列周期在一定区间范围内。而在限定周期范围内，嵌入点、嵌入波列的波高和长度的变化不会对模拟产生明显的影响，但是如果嵌入波列的波高较大，嵌入波列后波浪的最大波高会增大。