CN102982239B - 基于海浪谱的三维海浪模拟方法 - Google Patents

Abstract

本发明公开了一种基于海浪谱的三维海浪模拟方法，包括：确定观测时刻t，并取得该时刻的海浪波高；取观测时刻t前海面相对平静的某一时刻为观测的0时刻，并取得该时刻的海浪波高；在考虑海浪形成时间历程的基础上，推导得到海浪的相关函数，并计算得到海浪的衰减权值函数；将0时刻的海浪波高乘上衰减权值函数叠加在观测时刻t的海浪波高上，以此体现前一时刻海浪对观测时刻海浪的影响；改变叠加海浪的方向谱函数的角度初相位模拟风向变化；重复上述步骤，取多个时刻的海浪叠加，最终得到海浪的模拟。本发明通过上述步骤对海浪进行模拟，模拟出的海浪的效果相比传统方法失真小，符合现实的海浪。

Description

基于海浪谱的三维海浪模拟方法

技术领域

本发明涉及流体数值仿真领域，具体是一种基于海浪谱的三维海浪模拟方法。

背景技术

(1)研究海浪及其对工程的作用有三种途径：一是现场观测研究；二是在实验室内进行模拟研究；三是理论分析研究。由于海浪的复杂多变性，加上现场环境恶劣，进行现场观测需花费大量人力物力；理论研究目前也有较大的局限性，特别是对于不规则波浪，很多问题有赖于室内的模拟研究。随着电子计算机的发展和普及，海浪的数值模拟得到迅速发展，它具有经济方便等优点，日益受到人们的重视和广泛的应用

(2)在进行海洋运动体的运动学和动力学分析、海上工程施工过程规划及海上设备的设计时都需要研究海浪对目标设备的作用及其产生的影响，因此必须要有一个较理想的海浪信号来研究上述问题。海浪的仿真正是针对这种问题而提出的。目前学者们对于海浪模拟已经做了许多工作，例如以流体动力学方程组-Navier-Stokes方程组表达，通过一系列近似和简化，通过数值解法来求解。虽然这类方法相对而言，效果比较真实，适用范围较为广泛。但是计算量往往很大，效率较低。

(3)根据线性叠加法利用海浪谱和随机过程来模拟随时间任意变化的海浪。所谓谱分析就是阐明海浪的能量相对于波浪频率、波浪传播方向或其他独立变量的分布规律。是通过一系列的正弦函数的叠加来描述海浪，正弦函数的幅值、相位等等参数则由选定的海浪谱进行一一的取定。但是上述传统的模拟方法模拟出的海浪波峰波谷变化都比较突兀，这与实际的海浪情况是不符合，有较大的失真。

综上所述，传统模型仿真失真原因为：实际的海浪的波峰波谷应该比较柔和地过渡。传统海浪模拟失真，主要原因是将海浪看成平稳随机过程。实际中的海浪并不是随机平稳的，因为海浪的形成的主要作用因素——风是个很复杂的自然现象，是随时变化的。将海浪看成是随机平稳过程，必然导致把海浪形成过程孤立地考虑，而没考虑海浪的相邻的两个状态间的相互作用对于瞬时海浪的形成的影响，这必然导致模拟出来的海浪与现实不符。

发明内容

本发明针对现有技术中存在的上述不足，提供了一种基于海浪谱的三维海浪模拟方法。

本发明是通过以下技术方案实现的。

一种基于海浪谱的三维海浪模拟方法，包括以下步骤：

第一步，确定观测时刻t，并取得该时刻的海浪波高；

第二步，取观测时刻t前海面相对平静的某一时刻为观测的0时刻，并取得该时刻的海浪波高；

第三步，在考虑海浪形成时间历程的基础上，推导得到海浪的相关函数，并计算得到海浪的衰减权值函数；

第四步，将0时刻的海浪波高乘上衰减权值函数叠加在观测时刻t的海浪波高上，以此体现前一时刻海浪对观测时刻海浪的影响；

第五步，改变叠加海浪的方向谱函数的角度初相位模拟风向变化；

第六步，重复上述步骤，取多个时刻的海浪叠加，最终得到海浪的模拟。

所述第三步中，相关函数为：

$R\left(\mathrm{\τ}\right)=E\left[\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\right]$

$={\∫}_0^{+\∞}\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\mathrm{dt}$

$={\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)e^{-\mathrm{\ξ}{\mathrm{\ω}}_i(t+\mathrm{\τ})}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i){\mathrm{dt}}^{,}$

$={A_n}^2\left({\mathrm{\ω}}_i\right){\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i(2t+\mathrm{\τ})}\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

其中，ξ(t)表示考虑衰减情况下t时刻的瞬时波高，ξ(t+τ)表示考虑衰减情况下t+τ时刻的.瞬时波高，A_n(ω_i)表示振幅。

所述第三步中，衰减权值函数为：

其中：i表示第i次的叠加，n为根据模拟精度要求所取的叠加次数，为海浪在100s就基本衰减为0的情况下第时的相关函数值。

所述第四步中得到的观测时刻t的瞬时波高为：

$\mathrm{\ξ}(x,y,t)={\∫}_0^{t}R\left(\mathrm{\τ}\right)\×\mathrm{\ξ}(x,y,\mathrm{\τ})\mathrm{d\τ}$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{R}_i\left(t_i\right)\×{\mathrm{\ξ}}_i(x,y,t_i)$

$=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j}{\overset{n}{\mathrm{\Σ}}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_{i}x\cos {\mathrm{\θ}}_i+k_{i}y\sin {\mathrm{\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{\mathrm{ij}}),$ 其中，

R(τ)为衰减权函数，表征τ时刻的波高到t时刻时的衰减程度；

R_i(t_i)为t_i时刻的衰减权函数，表征t_i时刻的波高到t时刻时的衰减程度；

k_i为波数：k_i＝ω_i ²/g；

ε_ij为在0-2π之间服从均匀分布的随机初相位值；

ω_i的取法是在0.3-0.8之间分m段，并且在分成的区间内随机选取；

θ_j在-0.5π-0.5π之间均匀取，分n段；

A_n(ω_i，θ_j)为不同频率的简单正弦波的幅值。

所述幅值其中，S(ω_i，θ_j)＝S_n(ω_i)×G(ω_i，θ_j)为方向频谱函数，所述G(ω_i，θ_j)为方向函数，所述S_n(ω_i)为频谱函数。

所述频谱函数其中，g为重力加速度，U为可变的海面风速，α是表征波高及平均周期有关的常数，β是与平均周期和U有关的常数。

所述方向函数 $G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2{\mathrm{\θ}}_j+q\·\cos 4{\mathrm{\θ}}_j),$ 其中， $p=0.5+0.82e^{-\frac{1}{2}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right)}^4},$ $q=0.32e^{-\frac{1}{2}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right)}^4},$ 所述 $\left|\mathrm{\θ}\right|\≤\frac{\mathrm{\π}}{2}.$

本发明的原理为：在波浪的形成过程中，风起主导作用，风对波浪的作用主要沿着风向，其次有因为水体运动而相互挤压而激发的浪。如果某一时刻风向发生了变化(这在实际中的确是常常发生的)，其对于浪的作用的主方向将发生变化，但是上个时刻的状态在重力作用下并不能够即时地恢复到平衡态，从而对该时刻的海浪的形态是存在很大的影响的。目前的海浪模拟都只是在平衡状态的基础上进行平稳随机过程的叠加，从而不能够真实地反映海浪的真实形态。将海浪的形成沿时间轴展开，然后取其中的某一时刻的海浪，此时的海浪是在前一时刻海浪的基础上形成的。这样才能够反映实际的海浪的真实形态。由于前一时刻的海浪在重力作用下有所恢复。因此，在背景技术中传统模拟方法的基础上，采用加权叠加的方法，将上一时刻的海浪波高乘上一个代表衰减的权函数叠加在观测时刻的海浪波高上，以此来体现观测时刻前各个时刻海浪对该时刻海浪的影响。通过改变方向谱函数的角度初相位，模拟风向改变后的波浪，我们可以假设在比较短的一段时间内海风在90°范围内变化。海风的方向的变化可以通过将叠加海浪的初相位在(0，0.5π)内变化来实现，并取多个时刻的海浪相叠加，最终得到海浪的模拟。

求解权函数用到了随机振动理论中相关函数，在概率论中，有均值和方差的概念。但是均值和方差只是描述了随机过程单一时刻(即随机变量)的数学特征。随机振动理论中要描述两个不同时刻的随机变量之间的联系，因此引入了相关函数的概念。相关的定义，即两个变量之间的关系，因此研究两个变量之间的关系的问题即相关问题。相关函数分为1)自相关函数2)互相关函数，本发明主要用到的是自相关函数，其数学上的意义是描述一个随机过程在两个不同时刻之间的线性依赖关系。设随机过程x(t)在两个任意时刻t₁、t₂的随机变量为X₁(t)、X₂(t)。则X(t)的自相关函数定义为：

$R_x(t_1,t_2)=\lim \frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{x}_k\left(t_1\right)x_k\left(t_2\right)$

$R_x(t_1,t_2)=E\left[X\left(t_1\right)X\left(t_2\right)\right]$

$=\∫\underset{-\∞}{\overset{\∞}{\∫}}{x}_1{x}_{2}p(x_1,t_1;x_2,t_2)d{x}_{1}d{x}_2$

式中：p(x₁，t₁；x₂，t₂)是这两个随机变量的二维概率密度函数。

在本发明海浪的模拟中自相关函数的物理意义是：定量的表征随机函数x在t时刻和t+τ时刻的值x(t)和x(t+τ)的差别或相似程度的函数。以下通过合理简化后求出海浪的相关函数，并根据求得的相关函数来取权值函数。

由瞬时波高公式可得函数表达式，在该式中取(x，y)＝(0，0)处沿θ_j＝0方向上的波浪为研究对象来考虑海浪的相关函数。当风向由θ_j＝0改变为其他方向时，θ_j＝0方向上的浪失去风的激励，可以认为是有阻尼自由振动。其瞬时浪高随时间衰减，可得：

$\mathrm{\ξ}(\mathrm{0,0},t)=\underset{i}{\overset{m}{\mathrm{\Σ}}}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)$

海浪中两个频率不同的波是不相关的。故取出其中一个频率成分：

$\mathrm{\ξ}\left(t\right)=e^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)$

求该成分自相关函数：

$R\left(\mathrm{\τ}\right)=E\left[\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\right]$

$={\∫}_0^{+\∞}\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\mathrm{dt}$

$={\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)e^{-\mathrm{\ξ}{\mathrm{\ω}}_i(t+\mathrm{\τ})}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

$={A_n}^2\left({\mathrm{\ω}}_i\right){\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i(2t+\mathrm{\τ})}\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

计算上述积分可得：

$R\left(\mathrm{\τ}\right)=$

$\frac{1}{2}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i\mathrm{\τ}}[\frac{\mathrm{\ξ}}{2{\mathrm{\ω}}_i({\mathrm{\ξ}}^2+1)}\cos ({\mathrm{\ω}}_i\mathrm{\τ}+2{\mathrm{\ϵ}}_i)$

$-\frac{1}{2{\mathrm{\ω}}_i({\mathrm{\ξ}}^2+1)}\sin ({\mathrm{\ω}}_i\mathrm{\τ}+2{\mathrm{\ϵ}}_i)$

$+\frac{1}{2\mathrm{\ξ}{\mathrm{\ω}}_i}\cos \left({\mathrm{\ω}}_i\mathrm{\τ}\right)].$

从x轴上处画垂线交曲线于表示时的相关函数值。取两者的比值为代表衰减的权值函数可得：

$R_i=\frac{R\left(t\right)}{R\left(0\right)}=\frac{R\left(\frac{100}{n}i\right)}{R\left(0\right)}$

其中：

i根据模拟需要精度选取的叠加次数。

表征前以往时刻的波浪对该时刻的浪的影响的相关函数值可以在上述图线进行取值，根据上式计算衰减权值函数值。

本发明提供的基于海浪谱的三维海浪模拟方法，在考虑了海浪形成的时间历程的基础上，选择改变方向函数的初相位来表征风向的变化，通过推导得到海浪的相关函数，并且计算得到海浪的衰减函数值来表征之前不同时刻对于观测时刻海浪的不同影响。在这基础上对海浪进行模拟，模拟出的海浪的效果相比传统方法失真小，符合现实的海浪。附图说明

图1为本发明相关函数图线；

图2为传统方法模拟的海浪；其中，(a)为传统方法模拟的海浪的立体图，(b)为传统方法模拟的海浪x向视图；

图3为本发明取i＝3模拟出的海浪；其中，(a)为i＝3模拟的海浪的立体图，(b)为i＝3模拟的海浪x向视图；

图4为本发明取i＝4模拟出的海浪；其中，(a)为i＝4模拟的海浪的立体图，(b)为i＝4模拟的海浪x向视图；

图5为本发明取i＝5模拟出的海浪。其中，(a)为i＝5模拟的海浪的立体图，(b)为i＝5模拟的海浪x向视图；

图中，1为模拟的海浪。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1

本实施例包括以下步骤：

第一步，确定观测时刻t，并取得该时刻的海浪波高；

第二步，取观测时刻t前海面相对平静的某一时刻为观测的0时刻，并取得该时刻的海浪波高；

第三步，在考虑海浪形成时间历程的基础上，推导得到海浪的相关函数，并计算得到海浪的衰减权值函数；

第四步，将0时刻的海浪波高乘上衰减权值函数叠加在观测时刻t的海浪波高上，以此体现前一时刻海浪对观测时刻海浪的影响；

第五步，改变叠加海浪的方向谱函数的角度初相位模拟风向变化；

第六步，重复上述步骤，取多个时刻的海浪叠加，最终得到海浪的模拟。

所述第三步中，相关函数为：

$R\left(\mathrm{\τ}\right)=E\left[\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\right]$

$={\∫}_0^{+\∞}\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\mathrm{dt}$

$={\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)e^{-\mathrm{\ξ}{\mathrm{\ω}}_i(t+\mathrm{\τ})}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i){\mathrm{dt}}^{,}$

$={A_n}^2\left({\mathrm{\ω}}_i\right){\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i(2t+\mathrm{\τ})}\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

其中，ξ(t)表示考虑衰减情况下t时刻的瞬时波高，ξ(t+τ)表示考虑衰减情况下t+τ时刻的.瞬时波高，A_n(ω_i)表示振幅。

所述第三步中，衰减权值函数为：

其中：i表示第i次的叠加，n为根据模拟精度要求所取的叠加次数。表示的是海浪在100s就基本衰减为0的情况下第时的相关函数值，如图1所示。

所述第四步中得到的观测时刻t的瞬时波高为：

$\mathrm{\ξ}(x,y,t)={\∫}_0^{t}R\left(\mathrm{\τ}\right)\×\mathrm{\ξ}(x,y,\mathrm{\τ})\mathrm{d\τ}$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{R}_i\left(t_i\right)\×{\mathrm{\ξ}}_i(x,y,t_i)$

$=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j}{\overset{n}{\mathrm{\Σ}}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_{i}x\cos {\mathrm{\θ}}_i+k_{i}y\sin {\mathrm{\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{\mathrm{ij}}),$ 其中，

R(τ)为衰减权函数，表征τ时刻的波高到t时刻时的衰减程度；

R_i(t_i)为t_i时刻的衰减权函数，表征t_i时刻的波高到t时刻时的衰减程度；

k_i为波数：k_i＝ω_i ²/g；

ε_ij为在0-2π之间服从均匀分布的随机初相位值；

ω_i的取法是在0.3-0.8之间分m段，并且在分成的区间内随机选取；

θ_j在-0.5π-0.5π之间均匀取，分n段；

A_n(ω_i，θ_j)为不同频率的简单正弦波的幅值。

所述幅值其中，S(ω_i，θ_j)＝S_n(ω_i)×G(ω_i，θ_j)为方向频谱函数，所述G(ω_i，θ_j)为方向函数，所述S_n(ω_i)为频谱函数。

所述频谱函数其中，g为重力加速度，U为可变的海面风速，α是表征波高及平均周期有关的常数，β是与平均周期和U有关的常数。

所述方向函数 $G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2{\mathrm{\θ}}_j+q\·\cos 4{\mathrm{\θ}}_j),$ 其中， $p=0.5+0.82e^{-\frac{1}{2}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right)}^4},$ $q=0.32e^{-\frac{1}{2}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right)}^4},$ 所述 $\left|\mathrm{\θ}\right|\≤\frac{\mathrm{\π}}{2}.$

本实施例的原理为：在波浪的形成过程中，风起主导作用，风对波浪的作用主要沿着风向，其次有因为水体运动而相互挤压而激发的浪。如果某一时刻风向发生了变化(这在实际中的确是常常发生的)，其对于浪的作用的主方向将发生变化，但是上个时刻的状态在重力作用下并不能够即时地恢复到平衡态，从而对该时刻的海浪的形态是存在很大的影响的。目前的海浪模拟都只是在平衡状态的基础上进行平稳随机过程的叠加，从而不能够真实地反映海浪的真实形态。将海浪的形成沿时间轴展开，然后取其中的某一时刻的海浪，此时的海浪是在前一时刻海浪的基础上形成的。这样才能够反映实际的海浪的真实形态。由于前一时刻的海浪在重力作用下有所恢复。因此，在背景技术中传统模拟方法的基础上，采用加权叠加的方法，将上一时刻的海浪波高乘上一个代表衰减的权函数叠加在观测时刻的海浪波高上，以此来体现观测时刻前各个时刻海浪对该时刻海浪的影响。通过改变方向谱函数的角度初相位，模拟风向改变后的波浪，我们可以假设在比较短的一段时间内海风在90°范围内变化。海风的方向的变化可以通过将叠加海浪的初相位在(0，0.5π)内变化来实现，并取多个时刻的海浪相叠加，最终得到海浪的模拟。

求解权函数用到了随机振动理论中相关函数，在概率论中，有均值和方差的概念。但是均值和方差只是描述了随机过程单一时刻(即随机变量)的数学特征。随机振动理论中要描述两个不同时刻的随机变量之间的联系，因此引入了相关函数的概念。相关的定义，即两个变量之间的关系，因此研究两个变量之间的关系的问题即相关问题。相关函数分为1)自相关函数2)互相关函数，本发明主要用到的是自相关函数，其数学上的意义是描述一个随机过程在两个不同时刻之间的线性依赖关系。设随机过程x(t)在两个任意时刻t₁、t₂的随机变量为X₁(t)、X₂(t)。则X(t)的自相关函数定义为：

$R_x(t_1,t_2)=\lim \frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{x}_k\left(t_1\right)x_k\left(t_2\right)$

$R_x(t_1,t_2)=E\left[X\left(t_1\right)X\left(t_2\right)\right]$

$=\∫\underset{-\∞}{\overset{\∞}{\∫}}{x}_1{x}_{2}p(x_1,t_1;x_2,t_2)d{x}_{1}d{x}_2$

式中：p(x₁，t₁；x₂，t₂)是这两个随机变量的二维概率密度函数。

在本发明海浪的模拟中自相关函数的物理意义是：定量的表征随机函数x在t时刻和t+τ时刻的值x(t)和x(t+τ)的差别或相似程度的函数。以下通过合理简化后求出海浪的相关函数，并根据求得的相关函数来取权值函数。

由瞬时波高公式可得函数表达式，在该式中取(x，y)＝(0，0)处沿θ_j＝0方向上的波浪为研究对象来考虑海浪的相关函数。当风向由θ_j＝0改变为其他方向时，θ_j＝0方向上的浪失去风的激励，可以认为是有阻尼自由振动。其瞬时浪高随时间衰减，可得：

$\mathrm{\ξ}(\mathrm{0,0},t)=\underset{i}{\overset{m}{\mathrm{\Σ}}}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)$

海浪中两个频率不同的波是不相关的。故取出其中一个频率成分：

$\mathrm{\ξ}\left(t\right)=e^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)$

求该成分自相关函数：

$R\left(\mathrm{\τ}\right)=E\left[\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\right]$

$={\∫}_0^{+\∞}\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\mathrm{dt}$

$={\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)e^{-\mathrm{\ξ}{\mathrm{\ω}}_i(t+\mathrm{\τ})}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

$={A_n}^2\left({\mathrm{\ω}}_i\right){\∫}_0^{+\∞}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i(2t+\mathrm{\τ})}\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)\cos ({\mathrm{\ω}}_i(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)\mathrm{dt}$

计算上述积分可得：

$R\left(\mathrm{\τ}\right)=$

$\frac{1}{2}{e}^{-\mathrm{\ξ}{\mathrm{\ω}}_i\mathrm{\τ}}[\frac{\mathrm{\ξ}}{2{\mathrm{\ω}}_i({\mathrm{\ξ}}^2+1)}\cos ({\mathrm{\ω}}_i\mathrm{\τ}+2{\mathrm{\ϵ}}_i)$

$-\frac{1}{2{\mathrm{\ω}}_i({\mathrm{\ξ}}^2+1)}\sin ({\mathrm{\ω}}_i\mathrm{\τ}+2{\mathrm{\ϵ}}_i)$

$+\frac{1}{2\mathrm{\ξ}{\mathrm{\ω}}_i}\cos \left({\mathrm{\ω}}_i\mathrm{\τ}\right)]$

用matlab画出相关函数的图线如图1，其中取ξ＝0.05，由图1可见当100s以后基本衰减到零。

图1中曲线与y轴的交点为R(0)，表示的是τ＝0时的相关函数值，它表示的是函数与其本身的相似程度，因为没有比函数本身与函数更相似的函数，固其值是相关函数的最大值。

从x轴上处画垂线交曲线于表示时的相关函数值。取两者的比值为代表衰减的权值函数可得：

$R_i=\frac{R\left(t\right)}{R\left(0\right)}=\frac{R\left(\frac{100}{n}i\right)}{R\left(0\right)}$

其中：

i根据模拟需要精度选取的叠加次数。

表征前以往时刻的波浪对该时刻的浪的影响的相关函数值可以在图1的图线进行取值，根据上式计算衰减权值函数值。

本实施例根据不同的精度要求选取不同的n值模拟出海浪，并且与传统模拟模型模拟的海浪(如图2所示)进行比较，证明新模型的有效性。

实施例2

以我国南海某域海面为例，海浪表征流高为6m，平均周期是12.5s，在实施例1的基础上，可得：

α＝0.000665β＝0.35ω₀＝2πf＝0.503

在风速u＝36.7m/s时

取i＝3，t瞬时波高取：

ξ(xyt)＝R₁ξ₁(xyt)+R₂ξ₂(xyt₁)+R₃ξ₃(xyt₂)

其中：

ξ₁(xyt)是t时刻直接产生的瞬时波高；

ξ₂(xyt₁)为时刻t₁的瞬时波高；

ξ₃(xyt₂)为时刻t₂的波高

根据海浪瞬时波高公式可得到t时刻的瞬时波高：

${\mathrm{\ξ}}_1\left(\mathrm{x\; y\; t}\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j}{\overset{n}{\mathrm{\Σ}}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_{i}x\cos {\mathrm{\θ}}_i+k_{i}y\sin {\mathrm{\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{\mathrm{ij}})$

其中方向函数：

$G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2\mathrm{\θ}+q\·\cos 4\mathrm{\θ})$

根据海浪瞬时波高公式可得到t₁时刻的波高：

${\mathrm{\ξ}}_2\left(\mathrm{x\; y\; t}\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j}{\overset{n}{\mathrm{\Σ}}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_{i}x\cos {\mathrm{\θ}}_i+k_{i}y\sin {\mathrm{\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{\mathrm{ij}})$

其中方向函数：

$G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2({\mathrm{\θ}}_j+\frac{\mathrm{\π}}{6})+q\·\cos 4({\mathrm{\θ}}_j+\frac{\mathrm{\π}}{6}))$

方向函数里余弦的初相位取为

${\mathrm{\ξ}}_3\left(\mathrm{x\; y\; t}\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j}{\overset{n}{\mathrm{\Σ}}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_{i}x\cos {\mathrm{\θ}}_i+k_{i}y\sin {\mathrm{\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{\mathrm{ij}})---\left(18\right)$

其中方向函数：

$G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2({\mathrm{\θ}}_j+\frac{\mathrm{\π}}{3})+q\·\cos 4({\mathrm{\θ}}_j+\frac{\mathrm{\π}}{3}))$

方向函数里余弦的初相位由取为

将t₁、t₂两个时刻的瞬时波高ξ₂、ξ₃乘上表征衰减的权值函数R₂、R₃叠加在t时刻的瞬时波高上，得到t时刻的总的波高。

取i＝3在图3中查取0，t₁、t₂时刻的相关函数值。

根据权值函数公式计算可得： R₁＝1

Matlab模拟可得出的图形如图3，可以看出比图2的波峰波谷缓和了许多。

实施例3

实施例3为实施例2的变化例。

当取i＝4时：

根据权值函数公式计算可得

取 $R_4=\frac{-0.25}{6.382},$ $R_3=\frac{-0.045}{6.382},$

$R_2=\frac{1.607}{6.382},$ R₁＝1

模拟如图4，从图中可以看出模拟出的海浪比传统模拟方法模拟出的海浪缓和。

实施例4

实施例4为实施例2或实施例3的变化例。

当取i＝5时，

根据权值函数公式计算可得

$R_5=\frac{0.27}{6.382},$ $R_4=\frac{-0.6}{6.382},$

$R_3=\frac{1.323}{6.382},$ $R_2=\frac{-2.923}{6.382},$ R₁＝1

模拟如图5，从图中可以看出模拟出的海浪比传统模拟方法模拟出的海浪缓和

为了进一步说明新的模型的有效性，分别以传统模拟图和i＝5时的改进模拟图中最为突兀的部分进行定量比较。曲率表征的是曲线的变化的剧烈程度，可以反应出模拟得到的海浪的变化的突兀情况，通过计算该两个部分的曲率平均值并进行比较，可以定量地验证改进模型的有效性。

在传统模拟图中可以看出最突兀的部分大概是在：x＝(350，450)，y＝(550，650)这一部分区域。在该区域取100×100网格上每个节点处的波高的数值，通过matlab计算这部分10000个节点处曲率的平均值，最终得平均曲率为：1.8751×10⁶。

从i＝5的新模拟图中可以看出最突兀的部分大概在：x＝(550，650)，y＝(550，650)这一部分区域。在该区域取100×100网格上每个节点处的波高的数值，通过matlab计算这部分10000个节点处曲率的平均值，最终得平均曲率为：4.1196×10⁵。

可以看出i＝5的新模型的平均曲率要比传统模型的平均曲率小的多。说明其变化比传统模拟缓和。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (7)

1.一种基于海浪谱的三维海浪模拟方法，其特征在于，包括以下步骤：

第一步，确定观测时刻t，并取得该时刻的海浪波高；

第二步，取观测时刻t前海面相对平静的某一时刻为观测的0时刻，并取得该时刻的海浪波高；

第三步，在考虑海浪形成时间历程的基础上，推导得到海浪的相关函数，并计算得到海浪的衰减权值函数；

第四步，将0时刻的海浪波高乘上衰减权值函数叠加在观测时刻t的海浪波高上，形成观测时刻t的叠加海浪，以此体现观测时刻t前海面相对平静的某一时刻0海浪对观测时刻t海浪的影响；

第五步，改变叠加海浪的方向谱函数的角度初相位模拟风向变化；

第六步，重复上述步骤，取多个观测时刻t的叠加海浪，最终得到海浪的模拟。

2.根据权利要求1所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述第三步中，相关函数为：

$\begin{array}{c}R\left(\mathrm{\τ}\right)=E\[\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})\]\\ ={\∫}_0^{+\mathrm{\∞}}\mathrm{\ξ}\left(t\right)\mathrm{\ξ}(t+\mathrm{\τ})dt\\ ={\∫}_0^{+\mathrm{\∞}}{e}^{-{\mathrm{\ξ\ω}}_{i}t}{A}_n\left({\mathrm{\ω}}_i\right)\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)e^{-{\mathrm{\ξ\ω}}_i(t+\mathrm{\τ})}{A}_n\left({\mathrm{\ω}}_i\right)\×\cos \left({\mathrm{\ω}}_i\right(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)dt\\ ={A_n}^2\left({\mathrm{\ω}}_i\right){\∫}_0^{+\mathrm{\∞}}{e}^{-{\mathrm{\ξ\ω}}_i(2t+\mathrm{\τ})}\cos ({\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_i)\cos \left({\mathrm{\ω}}_i\right(t+\mathrm{\τ})+{\mathrm{\ϵ}}_i)dt\end{array},$

其中，ξ(t)表示考虑衰减情况下t时刻的瞬时波高，ξ(t+τ)表示考虑衰减情况下t+τ时刻的.瞬时波高，A_n(ω_i)表示振幅，ω_i为海浪圆频率，ε_i为初相位。

3.根据权利要求1所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述第三步中，衰减权值函数为：

其中：i表示第i次的叠加，n为根据模拟精度要求所取的叠加次数，为海浪在100s就基本衰减为0的情况下第时的相关函数值，R(t)为t时刻的相关函数值，R(0)为0时刻的相关函数值。

4.根据权利要求1所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述第四步中得到的观测时刻t的叠加海浪的瞬时波高为：

$\begin{array}{c}\mathrm{\ξ}(x,y,t)={\∫}_0^t{R}_i\left(\mathrm{\τ}\right)\×\mathrm{\ξ}(x,y,\mathrm{\τ})d\mathrm{\τ}\\ =\underset{i=1}{\overset{n}{\Σ}}{R}_i\left(t_i\right)\×{\mathrm{\ξ}}_i(x,y,t_i)\\ =\underset{i=1}{\overset{m}{\Σ}}\underset{j}{\overset{n}{\Σ}}{A}_n({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)\×\cos (k_i{\mathrm{xcos\θ}}_i+k_i{\mathrm{ysin\θ}}_i-{\mathrm{\ω}}_{i}t+{\mathrm{\ϵ}}_{ij})\end{array},$ 其中，

R_i(τ)为τ时刻的衰减权函数，表征τ时刻的波高到t时刻时的衰减程度；

R_i(t_i)为t_i时刻的衰减权函数，表征t_i时刻的波高到t时刻时的衰减程度；

k_i为波数：k_i＝ω_i ²/g；

ε_ij为在0-2π之间服从均匀分布的随机初相位值；

ω_i的取法是在0.3-0.8之间分m段，并且在分成的区间内随机选取；

θ_j在-0.5π-0.5π之间均匀取，分n段；

A_n(ω_i，θ_j)为不同频率的简单正弦波的幅值；

ξ_i(x，y，t_i)为t_i时刻的初始相位；

m为ω_i的分段数；n为θ_j的分段数。

5.根据权利要求4所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述幅值其中，S(ω_i，θ_j)＝S_n(ω_i)×G(ω_i，θ_j)为方向频谱函数，所述G(ω_i，θ_j)为方向函数，所述S_n(ω_i)为频谱函数；其中，Δω为ω_i的增量，Δθ为θ_j的增量。

6.根据权利要求5所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述频谱函数其中，g为重力加速度，U为可变的海面风速，α是表征波高及平均周期有关的常数，β是与平均周期和U有关的常数。

7.根据权利要求5所述的基于海浪谱的三维海浪模拟方法，其特征在于，所述方向函数 $G({\mathrm{\ω}}_i,{\mathrm{\θ}}_j)=\frac{1}{\mathrm{\π}}(1+p\·\cos 2{\mathrm{\θ}}_j+q\·\cos 4{\mathrm{\θ}}_j),$ 其中， $p=0.5+0.82e^{-\frac{1}{2}{\left(\frac{{\mathrm{\ω}}_i}{{\mathrm{\ω}}_0}\right)}^4},q=0.32e^{-\frac{1}{2}{\left(\frac{{\mathrm{\ω}}_i}{{\mathrm{\ω}}_0}\right)}^4},$ 所述其中，ω₀为所选海域的平均圆频率，θ_j表征方向的参数。