随机波浪速度及波浪力计算

%% 随机波速度和荷载计算 % Jonswap波浪谱生成随机波面 % 小振幅波浪理论速度解计算速度 % Morison方程计算单位长度圆柱结构体水平波浪荷载 % Code by JN_Cui 20221009 %% % ~~~~~~~~~~~~~~wave elevation~~~~~~__~~~~~~~~~~~~~ % | velocity u | | | % | acceleration ax | | | % | ===> | | | % | wave force | | | % d drag force Pd | | | % | inertia force Pi | | | % | |__| | % | diameter of cylinder D | % |_____________sea bottom__________________________| %% clear;clc H_s = 4;%有效波高 T_s = 7;%有效波高对应周期 dm = 0.1;%波浪谱频率步长 dt = 0.1;%时间步长 T = 100;%时间长度 dx = 1;%空间步长 X = 1;%空间长度 Cd = 1.2;%阻力系数 Cm = 2;%惯性系数 D = 5;%圆柱结构直径 rho = 1000;%海水密度 [S,Omega,omega_p,T_m] = Improved_Jonswap_spectral(H_s,T_s,dm); d = 2000;%总水深 dz = 5;%垂向步长 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; 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:N waitbar(i/(LT1+LT2)); if getappdata(Wait,'canceling') delete(Wait); 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); ZZ = 0:-dz:-d; u = zeros(length(ZZ),M,NN-1); w = zeros(length(ZZ),M,NN-1); ax = zeros(length(ZZ),M,NN-1); az = zeros(length(ZZ),M,NN-1); Pd = zeros(length(ZZ),M,NN-1); Pi = zeros(length(ZZ),M,NN-1); for j = 1:M waitbar((j+LT1)/(LT1+LT2)); if getappdata(Wait,'canceling') delete(Wait); break end for ii = 1:NN Eta(j,ii) = sum(a.*cos(k.*x(ii)-ww.*t(j)+e).*(a./L<=0.142)); for iii = 1:length(ZZ) Z = ZZ(iii); u(iii,j,ii) = sum(pi*2*a./T_w.*cosh(k.*(Z+d))./sinh(k.*d).*cos(k.*x(ii)-ww.*t(j)+e),'omitnan'); ax(iii,j,ii) = sum(2*pi^2.*2*a./T_w.^2.*cosh(k.*(Z+d))./sinh(k.*d).*sin(k.*x(ii)-ww.*t(j)+e),'omitnan'); w(iii,j,ii) = sum(pi*2*a./T_w.*sinh(k.*(Z+d))./sinh(k.*d).*sin(k.*x(ii)-ww.*t(j)+e),'omitnan'); az(iii,j,ii) = sum(-2*pi^2.*2*a./T_w.^2.*sinh(k.*(Z+d))./sinh(k.*d).*cos(k.*x(ii)-ww.*t(j)+e),'omitnan'); Pd(iii,j,ii) = 1/2*Cd*D*rho*u(iii,j,ii)^2*sign(u(iii,j,ii)); Pi(iii,j,ii) = rho*Cm*pi*D^2/4*ax(iii,j,ii); end end end delete(Wait); subplot(5,1,1) plot(Eta(:,1)) ylabel('\eta (m)'); subplot(5,1,2) plot(u(1,:,1)) ylabel('u (m/s)'); subplot(5,1,3) plot(ax(1,:,1)) ylabel('a_x (m/s^2)'); subplot(5,1,4) plot(Pd(1,:,1)) ylabel('P_d (N/m)'); subplot(5,1,5) plot(Pi(1,:,1)) ylabel('P_i (N/m)'); %% J谱函数 function [S,Omega,omega_p,T_p]=Improved_Jonswap_spectral(H_s,T_p,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 %% 二分法 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 %% J谱拟合 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)); 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