s_waveforce.rar_ittc单参数_s函数 海浪_波倾角_海浪 仿真_海浪波倾角

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 海浪模型，输出波倾角，波高 function [sys,x0,str,ts]=s_waveforce(t,x,u,flag) switch flag, case 0, [sys,x0,str,ts]=mdlInitializeSizes; case 3, sys=mdlOutputs(t,x,u); case {1,2,4,9}, sys=[]; otherwise error(['unhandled flag=',num2str(flag)]); end end function [sys,x0,str,ts]=mdlInitializeSizes sizes=simsizes; sizes.NumContStates =0; sizes.NumDiscStates =0; sizes.NumOutputs =1; %波倾角 sizes.NumInputs =3; % 有义波高，遭遇角， sizes.DirFeedthrough =1; sizes.NumSampleTimes =1; sys=simsizes(sizes); x0=[]; str=[]; ts=[0 0]; end function sys=mdlOutputs(t,x,u) g = 9.81; D=1458*10^3; %排水量 L=98.0; %船长 BB=10.2; %船宽 T=3.1; %吃水深度 h=1; %横稳心高 V=18*1.852*1000/3600; %航速18节，，，变换单位为m/s ,(=9.375m/s) 1节(kn)航速=1海里/时=(1852/3600)m/s =1.852千米(km)/h wo=0.698; %船舶横摇固有频率 (=0.8373) %wo=sqrt(0.4874); IxIx=D*h/wo^2; %惯量和附加惯量 nu=0.1325; %无因次衰减系数 即2*nu=0.3 Nu=nu*wo*IxIx; hh=u(1);%有义波高 g=9.81; %重力加速度 sida=u(2); %遭遇角 alfa=sida*pi/180; %遭遇角 角度变成弧度 *（pi/180） A=8.1*10^(-3)*g^2; B=3.11/hh^2; Et=0; %%%%%%%%%%%%%%%%%%%%%%%%%% Ett=0; %%%%%%%%%%%%%%%%%%%%%%%%%% Att=0; %%%%%%%%%%%%%%%%%%%%%%%%%% N=50; %%%%%%%%%%%%%%%%%%%%%%%%%% persistent EE%确保每个叠加的波的初相位相同 for i=1:N-1 %海浪谱的等分数 if u(3)==0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EE(i)= rand(1)*2*pi; %均匀分布在0--2*pi之间的随机初始相位 end end for i=1:N-1 ww(i)=(B/log(N/i))^(0.25); %%%自然频率 % EE(i)= rand(1)*2*pi; we(i)=ww(i)+ww(i)^2*V*cos(alfa)/g; %遭遇频率 %波面角修正系数 x1=BB*ww(i)^2/(2*pi*g); %%%%%%%%%%%%%% 从哪里得到 ？？？？？？？？？？？？？？？？？？？？ x2=T*ww(i)^2/(2*pi*g); %%%%%%%%%%%%%% 从哪里得到 ？？？？？？？？？？？？？？？？？？？？ if x1<=1.6 kb=-0.03565055*x2^3+0.7970778*x2^2-1.509432*x2+1; % kb 和 kt 是考虑了船宽和船舶吃水影响之后的修正系数 else kb=0.07781226; end if x2<=0.8 kt=0.3996212*x1^3-1.038799*x1^2+0.062684*x1+1; else kt=0.28443; end Kq=kt*kb; we(i)=ww(i)+ww(i)^2*V*cos(alfa)/g; %遭遇频率 Et=Et+sqrt(2*A/(4*N*B))*cos(ww(i)*t+EE(i));%长峰波海浪（随时间的幅值P49） Ett=Ett+sqrt(2*A/(4*N*B))*cos(we(i)*t+EE(i));%具有遭遇频率的波高仿真 Att=Att+(sqrt(2*A/(4*N*B))*cos(we(i)*t+EE(i))*Kq*ww(i)^2*sin(alfa)/g); %有效波倾角 end sys(1)=Att; %输出波倾角弧度 %sys(1)=Att*(180/pi);%输出波倾角角度 %sys(1)=Ett; %输出波高 end