等分频率法模拟随机波列(线性波叠加原理)

随机波浪理论 同时被 2 个专栏收录
5 篇文章 1 订阅
6 篇文章 0 订阅

线性叠加法

(大家反映程序没法运行,原来是之前没有放子程序的缘故,这里统一放上~)
海浪可看做一系列不同周期不同初相位的线性波叠加而成的:
η ( t ) = ∑ i = 1 M a i cos ⁡ ( k i x − ω i t + ϵ i ) \eta(t)=\sum\limits_{i=1}^{M}a_i\cos(k_ix-\omega_it+\epsilon_i) η(t)=i=1Maicos(kixωit+ϵi),
a i a_i ai为第 i i i个组成波的振幅, k i 和 ω i k_i和\omega_i kiωi为第 i i i个组成波的波数和圆频率。 ϵ i \epsilon_i ϵi ( 0 , 2 π ) (0,2\pi) (0,2π)之间的随机数,代表随机相位。假设靶谱的能量大多分布在区间 [ ω L ω H ] [\omega_L\quad\omega_H] [ωLωH],其他部分可忽略不计。将该区间平分为M个子区间,其间距为 Δ ω i = ω i − ω i − 1 \Delta\omega_i=\omega_i-\omega_{i-1} Δωi=ωiωi1,取
ω i ^ = ( ω i − 1 + ω i ) / 2 , a i = 2 S η η ( ω i ^ ) Δ ω i \hat{\omega_i}=(\omega_{i-1}+\omega_i)/2,a_i=\sqrt{2S_{\eta\eta}(\hat{\omega_i})\Delta\omega_i} ωi^=(ωi1+ωi)/2,ai=2Sηη(ωi^)Δωi ,
则将代表M个区间波能的M个线性波叠加起来,可得到海浪波面:
η ( t ) = ∑ i = 1 M 2 S η η ( ω i ^ ) Δ ω i cos ⁡ ( ω i ~ t + ϵ i ) \eta(t)=\sum\limits_{i=1}^{M}\sqrt{2S_{\eta\eta}(\hat{\omega_i})\Delta\omega_i}\cos(\tilde{\omega_i}t+\epsilon_i) η(t)=i=1M2Sηη(ωi^)Δωi cos(ωi~t+ϵi),
式中 ω i ~ \tilde{\omega_i} ωi~为第i个组成波的代表频率。

程序

主程序

% Randam wave simulation
% Designed by: JN-Cui
% Modified on 12/09/2019
%% DEFINITIONS
% alpha - energy scale factor; gama - spectral peak elevation factor;
% omega_m - spectral peak circular frequency; f_m - spectral peak frequency;
% U - wind speed at 10 m above sea surface; H_s - significant wave height;
% g - gravity acceleration;
%% FOR AVERAGE JONSWAP SPECTRAL
% gama=3.3; k=83.7; sigma_a=0.07; sigma_b=0.09;
% alpha=0.076*(X_bar)^(-0.22);
% X_bar=10^(-1)~20^(5); omega_m=22(g/U)*(X_bar)^(-0.33);
% f_m=3.5(g/U)(X_bar)^(-0.33);
%% FUNCTION

function [W_s,t,x,Omega,S]=random_wave_v1_0(H_s,T_s,dm,dt,T,dx,X)
[S,Omega,omega_p,T_m]=Improved_Jonswap_spectral(H_s,T_s,dm);
d=2000;
g=9.8;
N=200;
M=T/dt;
e=rand(1,N)*2*pi;
w=0:max(Omega)/(N-1):max(Omega);
ww(N)=max(w);
LT1=length(w);
LT2=M;
LT3=M;
for i=1:N-1
    ww(i)=unifrnd(w(i),w(i+1));
end
ww(N)=w(N);
Wait=waitbar(0,'程序计算中,请稍后', 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
setappdata(Wait,'canceling',0);
SS_o=zeros(1,length(w));
f=zeros(1,length(w));
T_w=zeros(1,length(w));
K=zeros(1,length(w));
k=zeros(1,length(w));
a=zeros(1,length(w));
L=zeros(1,length(w));
% for i=1:length(w)
for i=1:N
    waitbar(i/(LT1+LT2+LT3));
    if getappdata(Wait,'canceling')
        break
    end
    SS_o(i)=S_Improved_Jonswap_spectral(ww(i),H_s,T_s,dm);
    f(i)=ww(i)/2/pi;
    T_w(i)=1/f(i);
    K(i)=g*T_w(i)^2/(2*pi);
    fun=@(L1) L1/tanh(2*pi/L1*d)-K(i);
    L(i)=Division(fun,0.0001,0,K(i));
    k(i)=(2*pi/L(i));
    a(i)=sqrt(2*SS_o(i)*(w(2)-w(1)));
end
% Time steps
t=zeros(1,M);
for i=1:M
    t(i)=(i-1)*dt;
end
eta=zeros(N,1);
NN=floor(X/dx)+1;
x=0:dx:X;
Eta=zeros(M,NN-1);
for j=1:M
    waitbar((j+LT1)/(LT1+LT2+LT3));
    if getappdata(Wait,'canceling')
        break
    end
    for ii=1:NN
        Eta(j,ii)=sum(a.*cos(k.*x(ii)-ww.*t(j)+e));
    end
end
%%
% Wave Surface
W_s=Eta;
delete(Wait);
end

子程序1

%% DEFINITIONS
% alpha - energy scale factor; gama - spectral peak elevation factor;
% omega_m - spectral peak circular frequency; f_m - spectral peak frequency;
% U - wind speed at 10 m above sea surface; H_s - significant wave height;
% g - gravity acceleration;
%% FOR AVERAGE JONSWAP SPECTRAL
% gama=3.3; k=83.7; sigma_a=0.07; sigma_b=0.09;
% alpha=0.076*(X_bar)^(-0.22);
% X_bar=10^(-1)~20^(5); omega_m=22(g/U)*(X_bar)^(-0.33);
% f_m=3.5(g/U)(X_bar)^(-0.33);
%% IMPUT PARAMETERS
% H_s - significant wave height; T_s -  wave period at 1/3 wave height
% dm - calculation interval of omega
%% FUNCTION
function [S,Omega,omega_p,T_p]=Improved_Jonswap_spectral(H_s,T_p,dm)%the input of dm??
gama=3.3; sigma_a=0.07;sigma_b=0.09;
beta_j=0.06238/(0.23+0.0336*gama-0.185*(1.9+gama)^(-1))*(1.094-0.01915*log(gama));
T_s=T_p*(1-0.132*(gama+0.2)^(-0.559));
f_p=1/T_p;
omega_p=f_p*2*pi;
i=1;
df=dm/2/pi;
S_o=zeros(1,length(0:dm:1/T_s*2*pi*4));
Omega1=zeros(1,length(0:dm:1/T_s*2*pi*4));
for omega=0:dm:1/T_s*2*pi*4 %?显著波高对应频率的四倍  包含绝大部分谱的能量
    if omega<omega_p 
        sigma=sigma_a;
        S_o(i)=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...
            *gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);
    else
        sigma=sigma_b;
        S_o(i)=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...
            *gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);
    end
    Omega1(i)=omega; 
    i=i+1;
end
Omega=Omega1(2:end);
S=S_o(2:end);
end

子程序2

function [S_omega]=S_Improved_Jonswap_spectral(omega,H_s,T_s,dm)
gama=3.3; sigma_a=0.07;sigma_b=0.09;
beta_j=0.06238/(0.23+0.0336*gama-0.185*(1.9+gama)^(-1))*(1.094-0.01915*log(gama));
T_p=T_s/(1-0.132*(gama+0.2)^(-0.559));
% T_p=T_bar/(1-0.532*(gama+2.5)^(-0.569));
f_p=1/T_p;
omega_p=f_p*2*pi;
i=1;
df=dm/2/pi;
S_o=zeros(1,length(0:dm:1/T_s*2*pi*4));
Omega1=zeros(1,length(0:dm:1/T_s*2*pi*4));
    if omega<omega_p
        sigma=sigma_a;
        S_omega=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...
            *gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);
    else
        sigma=sigma_b;
        S_omega=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...
            *gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);
    end
end

子程序3

function output=Division(fun,x,a,b) 
p=-1; 
while (fun(a)*fun(b) <=0) && (abs(a-b)>x) 
    c=(a+b)/2; 
    if fun(c)*fun(b)<=0 
        a=c; 
        p=p+1; 
    else
        p=p+1;
        b=c;
    end
end
output=(a+b)/2; 
end
  • 4
    点赞
  • 10
    评论
  • 10
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

相关推荐
程序员的必经之路! 【限时优惠】 现在下单,还享四重好礼: 1、教学课件免费下载 2、课程案例代码免费下载 3、专属VIP学员群免费答疑 4、下单还送800元编程大礼包 【超实用课程内容】  根据《2019-2020年中国开发者调查报告》显示,超83%的开发者都在使用MySQL数据库。使用量大同时,掌握MySQL早已是运维、DBA的必备技能,甚至部IT开发岗位也要求对数据库使用和原理有深入的了解和掌握。 学习编程,你可能会犹豫选择 C++ 还是 Java;入门数据科学,你可能会纠结于选择 Python 还是 R;但无论如何, MySQL 都是 IT 从业人员不可或缺的技能!   套餐中一共包含2门MySQL数据库必学的核心课程共98课时   课程1:《MySQL数据库从入门到实战应用》   课程2:《高性能MySQL实战课》   【哪些人适合学习这门课程?】  1平时只接触了语言基础,并未学习任何数据库知识的人;  2对MySQL掌握程度薄弱的人,课程可以让你更好发挥MySQL最佳性能; 3想修炼更好的MySQL内功,工作中遇到高并发场景可以游刃有余; 4被面试官打破沙锅问到底的问题问到怀疑人生的应聘者。 【课程主要讲哪些内容?】 课程一:《MySQL数据库从入门到实战应用》 主要从基础篇,SQL语言篇、MySQL进阶篇三个角度展开讲解,帮助大家更加高效的管理MySQL数据库。 课程二:《高性能MySQL实战课》主要从高可用篇、MySQL8.0新特性篇,性能优化篇,面试篇四个角度展开讲解,帮助大家发挥MySQL的最佳性能的优化方,掌握如何处理海量业务数据和高并发请求 【你能收获到什么?】  1.基础再提高,针对MySQL核心知识点学透,用对; 2.能力再提高,日常工作中的代码换新貌,不怕问题; 3.面试再加,巴不得面试官打破沙锅问到底,竞争力MAX。 【课程如何观看?】  1、登录CSDN学院 APP 在我的课程中进行学习; 2、移动端:CSDN 学院APP注意不是CSDN APP哦  本课程为录播课,课程永久有效观看时长 【资料开放】 课件、课程案例代码完全开放给你,你可以根据所学知识,自行修改、优化。  下载方式:电脑登录课程观看页面,点击右侧课件,可进行课程资料的打包下载。
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值