%画回转圈的子程序，在第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
