海浪谱仿真

金河镇有黄金

已于 2023-10-08 19:13:05 修改

阅读量659

点赞数

文章标签：算法机器学习人工智能

于 2023-09-26 19:46:48 首次发布

本文链接：https://blog.csdn.net/lmmij/article/details/133322396

版权

海浪模型：
$\times cos(k \times \sqrt{(x-x_0)^2 + (y-y_0)^2 }- \omega t )$

%% -----------------PM谱
clear;close all;
x = 3:0.01:300;
N = length(x);
alpha = 0.0081;
beta = 0.74;
g = 9.8;
S = zeros(1,N);
for i =1:N
    k = 1/x(i);
    S(i) = alpha/4/k^4 * exp(-beta * g^2/k^2/10^4);
    S1(i) = alpha/4/k^4 * exp(-beta * g^2/k^2/7.5^4);
end
figure;plot(2*pi./x,S);
hold on;plot(2*pi./x,S1);hold off;

%% ---------JONSWAP谱
clear;close all;
L = 50000; %风区
U = 10; %风速
g = 9.8;
alpha = 0.076 * (g*L/U^2)^(-0.022);
beta = 0.74;
gamma = 3.3;
x = (100:10:L)/U;
N = length(x);
S = zeros(1,N);
for i = 1:N
    omega = 2*pi /x(i);
    omega_p = 0.12;
    if omega<omega_p
        sigma = 0.07;
    else 
        sigma = 0.09;
    end
    S(i) = alpha*g^2/omega^5 * exp(-1.25 * (omega_p/omega)^4) *gamma^(exp(-(1-(omega/omega_p))^2/2/sigma^2));

end
figure;plot(2*pi./x,S);

%%
%-----------------------Elfouhaily海浪杂波仿真
clear;
close all;
g = 9.8;
km = 370;
U_10 = 10;    %10米高的速度
x = 0.1:0.01:3000;
N = length(x);
SL = zeros(1,N);
SH = zeros(1,N);
for i = 1:N
    k = 1/x(i);   %空间波数1/m
    c = sqrt( g*(1 + k^2/km^2)/k);
    if x(i) <=1
        gamma = 1.7;
    else
        if (x(i)<= 5 && x(i)>1)
            gamma = 1.7 + 6*log10(omu);
        else
            gamma = 2.7 * omu ^0.57;
        end
    end
    k_p = g/U_10^2*sqrt(2*0.74/3);
    c_p = sqrt( g*(1 + k_p^2/km^2)/k_p);
    omu = U_10/c_p;
    if omu<5
        delta = 0.08 * (1+4/omu^3);
    else
        delta = 0.16;
    end
    G = exp(-(sqrt(k/k_p)-1)^2/2/delta^2);
    F_p = gamma^G * exp(-1.25 * (k_p/k)^2) *exp(-omu/10 *(sqrt(k/k_p)-1));
    alpha = 0.006 * sqrt(omu);
    SL(i) = alpha/2 *c_p /c *F_p/k^4;

    uf = 2;
    c_m = 0.23;
    k_m = 2*g/c_m ^2;
    k = 1/x(i);
    c = sqrt( g*(1 + k^2/km^2)/k);
    if uf <c_m
        alpha_m = 0.01*(1+log(uf/c_m));
    else 
        alpha_m =  0.01*(1+3 * log(uf/c_m));
    end
    F_m = exp(((k/k_m)-1)^2/4);
    SH(i) = alpha_m/2 *c_m/c*F_m;
end

figure;plot(log10(2*pi./x),log10(abs(SL)));
hold on;plot(log10(2*pi./x),log10(abs(SH+SL)));
ylim([-14,0]);