%% 随机波速度和荷载计算
% 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
评论0