集装箱船舰船的回转运动模拟附matlab代码.zip

%画回转圈的子程序，在第4章174米的船舶回转性仿真时用到 function dy=funcH(t,y) dy=zeros(7,1); L=320;%垂线间长 B=58;%船宽 D=28;%型深 d=20.8;%设计吃水 df=20;%首吃水 da=20.8;%尾吃水 dt=0;%吃水差 Cb=0.8098;%方形系数 V1=312622;%排水体积 AR=273;%舵面积 HR=15.799;%舵高 N1=1.82667;%舵展弦比 Dp=9.86;%螺旋桨直径 P=7.109;%螺距 N2=0.431;%盘面比 q=1.025;%海水密度 n1=1.103*10^(-6);%粘性系数 n=127/60;%主机转速 u=y(1); v=y(2); r=y(3); V=(u^2+v^2)^(0.5); v1=v/V;u1=u/V;r1=r*L/V; %无量纲化的惯性矩 m=q*Cb*L*B*d; Iz=m*L^2/16; ix=m/100*(0.398+11.98*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*d/B^2); iy=m*(0.882-0.54*Cb*(1-1.6*d/B)-0.156*(1-0.673*Cb)*L/B+... 0.826*L*d/B^2*(1-0.678*d/B)-0.638*L*d/B^2*Cb*(1-0.669*d/B)); iz=m*L^2/10000*(33-76.85*Cb*(1-0.784*Cb)+3.43*L/B*(1-0.784*Cb))^2; %纵向水动力导数 Rn=V*L/n1; Cf=0.075/(log10(Rn)-0.23)^2; Cr= (35*V1/L^3)^(0.5)/100;%查图谱 CAR=0.0004;%查表得到 Ct=Cf+Cr+CAR; S=V1^(2/3)*(3.432+0.305*L/B+0.443*B/d+0.643*Cb); %Xu1=-S/L/d*Ct*u1^2; %Xvv1=0.4*B/L-0.0006*L/d; t1=dt/d; iy1=iy/m; %Xvr1=(1.11*Cb-1.07)*(1+0.208*t1)*iy1; %Xrr1=0.0003*L/d; %横向水动力导数 %Yv1=-(pi*d/L+1.4*Cb*B/L)*(1+(25*Cb*B/L-2.25)*dt/d) %Yr1=(pi*d/(2*L))*(1+0.8*dt/d) % Nv1=-2*d/L*(1-dt/d) % Nr1=-(0.54*2*d/L-(2*d/L)^2)*(1+(34*Cb*B/L-3.4)*dt/d) % Yvv1=(-2.5*(1-Cb)*B/(d-0.5))*(1-(35.7*Cb*B/L-2.5)*dt/d) % Yvr1=-0.45+1.65*(1-Cb*d/B) % Yrr1=(0.343*Cb*d/B-0.07)*(1+(45*Cb*B/L-8.7)*dt/d) % Yvrr1=(-5.95*(1-Cb)*d/B)*(1+(40*(1-Cb)*d/B-2)*dt/d) % Yvvr1=(1.5*Cb*d/B-0.65)*(1+(110*(1-Cb)*d/B-9.7)*dt/d) % Nvv1=(0.96*(1-Cb)*d/B-0.066)*(1+(58*(1-Cb)*d/B-5)*dt/d) % Nrr1=(0.5*Cb*B/L-0.09)*(1-(30*Cb*B/L-2.6)*dt/d) % Nvrr1=(0.5*Cb*B/L-0.05)*(1+100*dt/d*(48*(Cb*B/L)^2-16*(Cb*B/L)+1.3)) % Nvvr1=(-57.5*(Cb*B/L)^2+18.4*Cb*B/L-1.6)*(1+(3*Cb*B/L-1)*dt/d) %Xvr1=(1.11*Cb-1.07)*(1+0.208*t1)*iy1; %Xrr1=0.0003*L/d; %横向水动力导数 Xu1=-0.0022; Xvv1=0.00274; Xvr1=0.0131; Xrr1=0.00058; Yv1=-0.0241; Yr1=0.00424; Nv1=-0.00794; Nr1=-0.00332; Yvv1=0.0023; Yvr1=-0.45+1.65*(1-Cb*d/B); Yrr1=0.00056; Yvrr1=-0.0403; Yvvr1=-0.0099; Nvv1=-0.00115; Nrr1=-0.00027; Nvrr1=0.00808; Nvvr1=-0.00337; %舵力计算 b=0.1*t*pi/180; if b<=y(7)*pi/180; w=b else w=y(7)*pi/180; aH=0.6874-1.3374*Cb+1.889*Cb^2; xH=-(0.4+0.1*Cb)*L; xR=-0.5*L; tR=0.2618+0.0539*Cb-0.1755*Cb^2; wp0=0.7*Cb-0.18; piao=atan(-v/u); piaop=piao+0.5*r*L/V; wp=wp0*exp(-4*piaop^2); up=u*(1-wp);%u表示船舶纵向分速度 nR=(2.1*(HR/Dp)-1.45)^2; k=0.6/nR; s=1-up/(n*P); if w>=0 %w表示舵角 K=1.065; else K=0.935; end nR1=Dp/HR; uR=nR*up*(1+K*nR1*k*(2-(2-k)*s)*s/(1-s)^2)^0.5; yR=1.163314-1.98283*Cb+1.390152*Cb^2; lR=-0.95*L; vR=yR*(v+lR*r);%v表示船舶横向分速度,r表示船舶角速度 UR=((0.85*u)^2+v^2)^0.5; aR=w-yR*(piao-2*xR*r/V); fa=6.13*N1/(2.25+N1); FN=-0.5*q*AR*fa*UR^2*sin(aR); XR=(1-tR)*FN*cos(w); YR=(1+aH)*FN*cos(w); NR=(xR+aH*xH)*FN*cos(w); %螺旋桨的作用 K1=0.0821;K2=-0.027;K3=-0.2768;K4=0.3499;%待定 J=up/n/Dp;tp=0.5*Cb-0.15; KT=K1*J^3+K2*J^2+K3*J+K4; XP=(1-tp)*q*n^2*Dp^4*KT; %解方程 XH=(Xu1+Xvv1*v1*v1+Xvr1*v1*r1+Xrr1*r1*r1)*(0.5*q*V^2*L*d); YH=(Yv1*v1+Yr1*r1+Yvv1*abs(v1)*v1+Yrr1*abs(r1)*r1+Yvr1*abs(v1)*r1+... Yvvr1*v1^2*r1+Yvrr1*v1*r1^2)*(0.5*q*V^2*L*d); NH=(Nv1*v1+Nr1*r1+Nrr1*abs(r1)*r1+Nvv1*abs(v1)*v1+... Nvvr1*(v1)^2*r1+Nvrr1*v1*(r1)^2)*(0.5*q*V^2*L^2*d); xc=2; %dy(1)=(XH+XR+(m+iy)*v*r)/(m+ix); dy(1)=0; dy(2)=(YH+YR-(m+ix)*y(1)*y(3)-(m+ix+iy)*xc*dy(3))/(m+iy); dy(3)=(NH+NR-YH*xc-(m+ix+iy)*(dy(2)+y(1)*y(3))*xc)/(Iz+iz); dy(4)=y(3); dy(5)=u*cos(y(4))-v*sin(y(4)); dy(6)=v*cos(y(4))+u*sin(y(4)); end