【船舶运动仿真】水面船舶三度运动仿真（风浪流模型）【含Matlab源码 3491期】

本文链接：https://blog.csdn.net/weixin_63266434/article/details/135569184

在这里插入图片描述

⛄一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【船舶运动仿真】基于matlab水面船舶三度运动仿真（风浪流模型）【含Matlab源码 3491期】
点击上面蓝色字体，直接付费下载，即可。

备注：
点击上面蓝色字体付费专栏Matlab物理应用（初级版），扫描上面二维码，付费29.9元订阅海神之光博客付费专栏Matlab物理应用（初级版），凭支付凭证，私信博主，可免费获得1份本博客上传CSDN资源代码（有效期为订阅日起，三天内有效）；
点击CSDN资源下载链接：1份本博客上传CSDN资源代码

⛄二、部分源代码

%%%%%%%%%%%% 20160515 wds %%%%%%%%%%%%%
clc
clear
% system data
global A_D; A_D = pi/180;
% ship data
m =4.3e53.85417^3; %4.3e5kg质量
x_g = -0.01373.85417; %-0.0137m重心X点坐标
Iz = 5.4e63.85417^5; %kg.s^2.m2转动惯量Iz
rho = 1025; %kg.m^-3 海水密度
L = 483.85417; %48m船长

nu = zeros(3,1); % velocity
eta = zeros(3,1); % position
% tao=[200,0,0]‘;
% tao=[0,200,0]’;
tao=[0,0,200]';

%无因次水动力系数
X_u_neg = 0;
Y_v_neg = -1.566e-2;
Y_r_neg = 1.922e-3;
N_v_neg = -1.78961988e-3;
N_r_neg = -2.42e-3;

X_du_neg = -9.469e-4;
Y_dv_neg = -8.316e-3;
N_dr_neg = -3.647e-4;
Y_dr_neg = -8.078e-4;

% dimensional hydrodynamic coefficients 有因次化
X_du = X_du_neg0.5rhoL^3;
Y_dv = Y_dv_neg0.5rhoL^3;
N_dr = N_dr_neg0.5rhoL^5;
Y_dr = Y_dr_neg0.5rhoL^4;

X_u = X_u_neg0.5rhoL^2;
Y_v = Y_v_neg0.5rhoL^2;
Y_r = Y_r_neg0.5rhoL^3;
N_v = N_v_neg0.5rhoL^3;
N_r = N_r_neg0.5rho*L^4;

% system inertia matrix
M = zeros(3,3);
M(1,1) = m-X_du;
M(2,2) = m-Y_dv;
M(2,3) = mx_g-Y_dr;
M(3,2) = mx_g-Y_dr;
M(3,3) = Iz-N_dr;

% damping matrix
D = zeros(3,3);
D(1,1) = -X_u;
D(2,2) = -Y_v;
D(2,3) = -Y_r;
D(3,2) = -N_v;
D(3,3) = -N_r;

h=0.3;
t=0:300;

for i=1:length(t)
%%transfer matrix
J_eta=zeros(3,3);
J_eta(1,1)=cos(eta(3));
J_eta(1,2)=-sin(eta(3));
J_eta(2,1)=sin(eta(3));
J_eta(2,2)=cos(eta(3));
J_eta(3,3)=1;

%科里奥利向心力矩阵   低速时C(v)v的二次项可以忽略，较高航速下不再成立
C=zeros(3,3);
C(1,3)=-(m-Y_dv)*nu(2,1)-(m*x_g-Y_dr)*nu(3,1);
C(2,3)=(m-X_du)*nu(1,1);
C(3,1)=(m-Y_dv)*nu(2,1)+(m*x_g-Y_dr)*nu(3,1);
C(3,2)=-(m-X_du)*nu(1,1);
% 环境干扰 W 
W =  zeros(3,1);
windspeed = 10;             % 风的绝对速度
angle_w = 225*pi/180;       % 风的绝对角度，相对于左舷为正 
v_cur =  1*0.8;             % 流速
angle_c = 10*pi/180;        % 流向角
angle_wave = 15*pi/180;     % 浪向角

F_wind = F_feng(windspeed, angle_w, eta);        
F_cur =  F_current(v_cur, angle_c, eta);   
F_wave= F_lang(windspeed,angle_wave,eta);
W = force(F_cur, F_wind, F_wave);

%%%迭代解算
k11=J_etanu;
k21=-inv(M)Dnu+inv(M)tao;
k12=J_eta(nu+h/2k21);
k22=-inv(M)D(nu+h/2k21);
k13=J_eta(nu+h/2k22);
k23=-inv(M)D(nu+h/2k22);
k14=J_eta*(nu+hk23);
k24=-inv(M)D(nu+hk23);
eta=eta+h/6*(k11+2k12+2k13+k14);
nu=nu+h/6*(k21+2k22+2k23+k24);

eta_draw(:,i)=eta;
nu_draw(:,i)=nu;

end
figure(1)
subplot(321);
plot(t,eta_draw(1,:),‘k-’,‘linewidth’,1); ylabel(‘北面(m)’);
hold on;grid on;
subplot(322);
plot(t,nu_draw(1,:),‘k-’,‘linewidth’,1);ylabel(‘纵向速度’);
hold on;grid on;
subplot(323);
plot(t,eta_draw(2,:),‘k-’,‘linewidth’,1);ylabel(‘东面(m)’);
hold on;grid on;
subplot(324);
plot(t,nu_draw(2,:),‘k-’,‘linewidth’,1);ylabel(‘横向速度’);
hold on;grid on;
subplot(325);
plot(t,eta_draw(3,:),‘k-’,‘linewidth’,1);ylabel(‘艏向角’);
hold on;grid on;
subplot(326);
plot(t,nu_draw(3,:),‘k-’,‘linewidth’,1);ylabel(‘艏向角速度’);
hold on;grid on;