# 等分频率法模拟随机波列（线性波叠加原理）

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) ,
a i a_i 为第 i i 个组成波的振幅， k i 和 ω i k_i和\omega_i 为第 i i 个组成波的波数和圆频率。 ϵ i \epsilon_i ( 0 , 2 π ) (0,2\pi) 之间的随机数，代表随机相位。假设靶谱的能量大多分布在区间 [ ω L ω H ] [\omega_L\quad\omega_H] ,其他部分可忽略不计。将该区间平分为M个子区间，其间距为 Δ ω i = ω i − ω i − 1 \Delta\omega_i=\omega_i-\omega_{i-1} ，取
ω 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} ,

η ( 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) ,

# 主程序

% 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
收藏
• 一键三连
• 扫一扫，分享海报

04-24 2万+
03-15 4391
09-28 46万+
12-01 1万+
08-27
05-03 3349
10-21 1万+
07-09 2413
08-28