Matlab模拟船舶三自由度MMG模型

%% 函数定义 %函数声明行 %利用四阶龙哥库塔求解MMG运动方程 %%主函数 function dy=mmg111(t,y) u=y(1);%X轴向速度 v=y(2);%Y轴向速度 r=y(3);%转首角速度 xg=y(4);%X坐标 yg=y(5);%Y坐标 course=y(6);%航向 rudder=y(7);%舵角 %参数赋值（实船数据） %参数取自《船舶运动数学模型》P162，商船“怀集河”的参数 p=1.025;%水密度 L=114;%船长 Lw=114;%设计水线长度 B=20.5;%船宽 D=8.7;%型深 d=6.7;%设计吃水 Cb=0.694;%方形系数 Cp=0.708;%菱形系数，没有找到该数据，大致根据经验确定值 Vpai=10805.7;%排水量 m=10805.7*p;%船舶质量 M=m/(0.5*p*L^2*d);%无量纲化后的船舶质量 Ar=17.95;%舵面积 Hr=5.49;%舵高 Dp=4.5;%螺旋桨直径 P=3.99;%螺旋桨螺距 zhanbi=1.68;%展铉比 n=126/60;%螺旋桨转数 Ae=0.71;%盘面比 Z=4;%桨叶数 %船舶运动初始状态赋值 V=(v^2+u^2)^0.5;%船舶运动速度 u0=13.5*1852/3600;%船舶航速 vnian=1.18831*10^(-6);%流体运动粘性系数,经验系数 %利用周昭明回归公式求解船舶附加质量Mx、My %求取附加质量Mx、My，因船舶质量取无量纲化后的质量M，因此无需再对Mx、My无量纲化 Mx=M/100*(0.398+11.97*Cb*(1+3.73*d/B)-2.89*Cb*L/B*(1+1.13*d/B)+0.175*Cb*(L/B)^2*(1+0.541*d/B)-1.107*L/B*d/B); My=M*(0.882-0.54*Cb*(1-1.6*d/B)-0.156*L/B*(1-0.673*Cb)+0.826*d/B*L/B*(1-0.678*Cb/B)-0.638*Cb*d/B*L/B*(1-0.669*d/B)); %求取惯性力矩I和附加惯性力矩J，无量纲化后分别为Izz和Jzz I=(1+Cb^4.5)*p*m/24*(L^2+B^2);%经验公式 Izz=I/(0.5*p*L^4*d);%无量纲化后结果 Jzz=0.01*(33-76.85*Cb*(1-0.784*Cb)+3.43*L/B*(1-0.63*Cb))*M/2.6;%经验公式求解（不理解无量纲化的方法） %利用桑地公式求取船体浸湿面积S %桑地公式S=k*(Vpai*Lw)^0.5，k为面积系数，不好选取，故采用单桨商船的泰勒公式 S=(3.432+0.305*L/B+0.443*B/d-0.643*Cb)*Vpai^(2/3); %求取雷诺数 Rn=V*L/vnian;%V为航速，L为船长，vnian为流体运动粘性系数 %求取摩擦阻力系数Cf Cf=0.4631/(log10(Rn))^2;%桑海公式 %求取剩余阻力系数Cr Cr=(0.0067*V^4-0.166*V^3+1.553*V^2-6.3699*V+10.02)*10^(-3)*0.85/0.95;%经验公式 %确定粗糙度补粘系数 Car Car=0.000368; %求取船舶总阻力系数Ct Ct=Cf+Cr+Car; %求取直航阻力Xuu Xuu=-S/(L*d)*Ct; %其它流体力导数的估算 Xvv=0.4*B/L-0.006*L/d; Xvr=(1.11*Cb-0.07)*My; Xrr=0.0003*L/d; %求取纵向粘性流体力XH XH=Xuu*(u/u0)^2+Xvv*(v/u0)^2+Xvr*(v/u0)*(r*L/u0)+Xrr*(r*L/u0)^2; %求直航时的推力减额系数tp0 tp0=0.5*Cp-0.12;%对于单桨标准商船（Cb=0.54~0.84２），采用汉克歇尔公式 tp0per=tp0; %求解漂角piaojiao piaojiao=atan(-v/u); %求解舵处漂角piaojiaoR lr=-0.95; piaojiaoR=piaojiao-lr*r*L/V;%lr为经验系数应取-0.9~-1.0，此处取-0.95 %求解仅有船舶纵向影响下的伴流系数wp0 wp0=0.70*Cp-0.18;%对于单桨标准商船（Cb=0.54~0.842），采用汉克歇尔公式 % 求解斜流影响的伴流系数wp piaojiaoP=piaojiao+0.5*L/u0*r; wp=wp0*exp(-4*piaojiaoP^2); %求解lcb Xc=0;%浮心纵向坐标 lcb=Xc/(100*L); %求解ra ra=(B/d)*(1.3*(1-Cb)-3.1*lcb); % 求解kt kt=0.00023*(ra*L/Dp)-0.028; %推力减额系数tp f=kt*piaojiaoR; tp=tp0per*exp(-4*piaojiaoP^2); %求取螺旋桨推力T up=(1-wp)*u; J=up/(n*Dp); Kt=0.0821*J^3-0.2308*J^2-0.2768*J+0.3499;%螺旋桨推力系数，一般随船已知 T=p*n^2*Dp^4*Kt; %求取推力XP XP=(1-tp)*T/(0.5*p*L*d*V^2);%有效推力 %求取舵的阻力减额系数tr tr=0.2618+0.0539*Cb-0.1755*Cb^2; %求取参数CP yxbi=Dp/Hr; up0=u*(1-wp0); s0=1-up0/(n*P); s=1-(1-wp)*u/(n*P); CP=1/(1+0.6*yxbi*(2-1.4*s)*s/(1-s)^2)^0.5; %求取参数CS if piaojiaoR<=0.5/0.45 CS=0.45*piaojiaoR; else CS=0.5; end %求取参数zx zx=CP*CS; %求取参数舵角rudder rudder0=-2*s0*pi/180; rudder=35*pi/180;%实际舵角 rudderE=35*pi/180;%命令舵角 if abs(rudderE-rudder)/2.5>3*pi/180 rudder=3*pi/180; else rudder=(rudderE-rudder)/2.5; end %求取舵有效来冲角ar ar=rudder-rudder0-zx*piaojiaoR; %求取舵处来流的有效速度Ur if rudder-rudder0>=0 C=1.065 ;% 向右转时 else C=0.935;% 向左转时 end wr0=0.25; wr=wr0*exp(-4*piaojiaoR); k=0.6*(1-wp)/(1-wr); Gs=yxbi*k*(2-(2-k)*s)*s/(1-s)^2; Ur=V*(1-wr)*(1+C*Gs)^0.5; %求舵的升力系数 fa=6.13*zhanbi/(2.25+zhanbi); %求取舵的正压力Fn Fn=-0.5*p*Ar*fa*Ur^2*sin(ar)/(0.5*p*L*d*V^2); %求取舵力XR XR=(1-tr)*Fn*sin(rudder); %求取YH Yv=-pi*d/L*(1.0+0.4*Cb*B/d); Yr=-pi*d/L*(-0.5+2.2*B/L+0.08*B/d); Yvv=0.048265-6.293*(1-Cb)*d/B; Yvr=-0.3791+1.28*(1-Cb)*d/B; Yrr=0.0045-0.445*(1-Cb)*d/B; YH=Yv*v/V+Yr*r*L/V+Yvv*v/V*abs(v/V)+Yvr*abs(v/V)*r*L/V+Yrr*r*L/V*abs(r*L/V); %求取NH Nv=-pi*d/L*(0.5+0.24*d/L); Nr=-pi*d/L*(0.25-0.56*B/L+0.039*B/d); Nr=-pi*d/L*(0.25-0.56*B/L+0.039*B/d); Nvvr=-6.0856+137.4735*(Cb*B/L)-1029.514*(Cb*B/L)^2+2480.6082*(Cb*B/L)^3; Nvrr=-0.0635+0.04414*(Cb*d/B); Nrr=-0.0805+8.6092*(Cb*B/L)^2-36.9816*(Cb*B/L)^3; NH=Nv*v/u0+Nr*r*L/u0+Nvvr*(v/u0)^2*r*L/u0+Nvrr*v/u0*(r*L/u0)^2+Nrr*r*L/u0*abs(r*L/u0); %求取系数aH0 aH0=0.6784-1.3374*Cb+1.8891*Cb^2; if J<=0.3 aH=aH0*J/0.3; else aH=aH0; end %求取YR YR=(1+aH)*Fn*cos(rudder); %求取NR xH=-(0.4+0.1*Cb); NR=(-0.5+aH*xH)*Fn*cos(rudder); %%分函数 y1=((Mx+My)*v/u0*r*L/u0+XH+XP+XR)/(Mx+My) ; %X轴方向速度 y2=(-(Mx+My)*u/V*r*L/V+YH+YR)/(Mx+My);%Y轴方向速度 y3=(NH+NR+YH)/(Izz+Jzz);%角加速度 y4=(u*cos(course)-v*sin(course))/L;%X坐标 y5=(u*sin(course)+v*cos(course))/L;%Y坐标 y6=r;%转首角速度 y7=(rudderE-rudder)/2.5;%舵角变化率 dy=[y1;y2;y3;y4/V;y5/V;y6;y7];