%画回转圈的子程序,在第4章174米的船舶回转性仿真时用到
function dy=USVMMG(t,y)
Uc=0;%水流流速
angle_c=20*pi/180; %流向
dy=zeros(10,1);
M=7.066;
Lc=9.05;%垂线间长
L=10.5;%船长
B=2.88;%船宽
D=0.825;%型深
df=0.64;%首吃水
da=0.66;%尾吃水
d=(df+da)/2;%平均吃水
dt=da-df;%吃水差
q=1.025;%海水密度
V1=M/q;%排水体积
Cb= V1/(B*Lc*d);%方形系数
Cp=0.62;%菱形系数,没有找到该数据,大致根据经验确定值
HR=0.58;%舵高
BR=0.36;%舵宽
N1=HR/BR;%舵展弦比
AR=HR*BR;%舵叶面积
Dp=0.46;%螺旋桨直径
P=0.42;%螺距
N2=0.431;%盘面比
n1=1.103*10^(-6);%粘性系数
nt=y(9);
nmax=y(10);
if nt>nmax
nt=nmax;
end
n=nt/60;%主机转速
ur=y(1);
vr=y(2);
r=y(3);
sigma=y(4);
theta=y(7);
b=y(8);
%水流对船的影响
uc=Uc*cos(angle_c-sigma);
vc=Uc*sin(angle_c-sigma);
u=ur-uc;
v=vr-vc;
V=(u^2+v^2)^(0.5);
if V==0
v1=0;u1=0;r1=0;
else
v1=v/V;u1=u/V;r1=r*L/V;
end
%无量纲化的惯性矩
m=q*Cb*L*B*d;
Iz=m*(1+Cb^4.5)+(L^2+B^2.4)/24;
%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=L^2*0.01*(33-76.85*Cb*(1-0.784*Cb)+3.43*L/B*(1-0.63*Cb));
%纵向水动力导数
Rn=V*L/n1;%雷诺数
Cf=0.075/(log10(Rn)-2.03)^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;
t1=dt/d;
iy1=iy/m;
%Xvr1=(1.11*Cb-1.07)*(1+0.208*t1)*iy1;
%Xrr1=0.0003*L/d;
%横倾力矩计算
%Np=0.12*sqrt((m*(B^2+4*ZG^2)/(12))*H*m*9.8);
%GZFAI=H*sin(FAI);
%ZH=ZG-(4-B/D+0.02*(B/D-5.35)^3)*D;
%KH=-Np*p-m*9.8*GZFAI-YH*ZH;
%横向水动力导数
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-0.27*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);
% Nvfai=0.5*midu*L^2*D*V*(-1.72*K);
% Nrfai=0.5*midu*L^3*D*V*(2.6*(0.54*K-K^2));
% Nfai=0.5*midu*L^2*D*V^2*(-0.008);
Xvv1=0.4*B/L-0.0006*L/d;
Xvr1=(1.11*Cb-0.07)*(1+0.208*t1)*iy1-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=theta+(y(7)/abs(y(7)))*15*t*pi/180;
%if abs(b)<=abs(theta)
w=b;
%else
% b=y(7);
%end
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.5*Cb-0.05; %泰洛公式,单桨
%wp0=0.55*Cb-0.20; %双桨
%wp0=0.7*Cp-0.18;
%piao=atan(-v/u);
if u==0
piao=0;
else
piao=atan(-v/u);
end
piaop=piao+0.5*r1;
wp=wp0*exp(-4*piaop^2);
up=u*(1-wp);%u表示船舶纵向分速度
%nR=(2.1*(HR/Dp)-1.45)^2;
nR=1.08;
k=0.6/nR;
if n==0
s=0;
else
s=1-up/(n*P);
end
if w>=0 %w表示舵角
K=1.065;
else
K=0.935;
end
nR1=Dp/HR;
%Gs=nR1*k*(2-(2-K)*s)*s/(1-s)^2;
%uR=up*(1+K*Gs)^0.5;
if s==1
uR=0;
else
%uR=up*(1+(6.3*s)^1.5)^0.5;
uR=up*(1+K*nR1*k*(2-(2-k)*s)*s/(1-s)^2)^0.5;
end
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;
UR=sqrt(uR^2+vR^2);
%aR=w-yR*(piao-2*xR*r1/L);
if uR==0
aR=0;
else
aR=w-atan(-vR/uR);
end
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;%待定
K1=0.1736;K2=-0.4464;K3=-0.04861;K4=0.4952;%待定
if n==0
J=0;
else
J=up/(n*Dp);
end
%tp=0.5*Cb-0.15;
tp0=0.77*Cp-0.30;
tp=tp0*exp(-4*piaop^2);
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=0;
%风浪流的影响
eta = [y(5);y(6);real(y(4))];
eta(3)=rem(eta(3),2*pi);
W = zeros(3,1);
%windspeed = 15; % 风的绝对速度
%angle_w = 90*pi/180; % 风的绝对角度,相对于左舷为正
%v_cur = 1*0.8; % 流速
%angle_c = 90*pi/180; % 流向角
angle_wave = 90*pi/180; % 浪向角
%F_wind = F_feng1(windspeed, angle_w,eta,u,v,piao);
%F_cur = F_current1(v_cur, angle_c, eta,u,v,piao);
%F_wave = F_lang1(windspeed, angle_wave,eta);
%F_wave=[0;0;0];
%W = force1([0,0,0]', F_wind, F_wave)/1000;
%W=F_wind/1000;
%W = zeros(3,1);
dy(1)=(XH+XP+XR+W(1)+(m+iy)*v*r)/(m+ix)+r*vc;
dy(2)=(YH+YR+W(2)-(m+ix)*y(1)*y(3)-(m+ix+iy)*xc*dy(3))/(m+iy)+r*uc;
dy(3)=(NH+NR+YH+W(3)-(m+ix+iy)*(dy(2)+y(1)*y(3))*xc)/(Iz+iz);
%dy(3)=(1+YH)/(Iz+iz);
dy(4)=y(3);
dy(5)=v*cos(y(4))+u*sin(y(4));
dy(6)=u*cos(y(4))-v*sin(y(4));
dy(8)=(theta-w)/2.5;
%dy(9)=100;