doubletank.rar_matlab ship roll_减摇_减摇水舱_制导_船舶

% 双水舱被动减摇系统减摇原理仿真 clear; clc; L=168 ; %(船舶设计水线长) B=28; %(船舶设计水线宽) Hw=3.25; %(波浪有义波高) Tw=6.0; %(波浪特征周期) Mctr=1; %(波谱选择控制参数) zg=10; %(船舶重心坐标） pi=3.14; g=9.8; p=1; %(海水密度) D=18378*g; %(船的排水量) ho=2.28; %(横稳心高) Vo=0.12; %(船的无因次阻尼系数) %Vo=0.042*(1+5.0*v) %(船的无因次阻尼系数) Ix=(D/(12*g))*(B^2+4*zg*zg); %(船舶通过纵轴的惯量) Jx=0.29*Ix; %（附连质量对x轴的惯性矩，根据‘船舶耐波性’第77页公式计算） I1=Ix+Jx; %Wo=sqrt((D*ho)/I1); %(船舶固有圆频率) Wo=0.42; for i=1:1:200; %(海浪遭遇频率) %for i=1:1:50 %(海浪遭遇频率) W(i)=i/200 ; At(i)=W(i)./Wo ; %(相对圆频率) TT=2*3.14/Wo; %(船舶固有周期) %%%%%%%%%%%%%%%%%%%%%%%%%% 裸船体横摇运动频率响应分析 Sak(i)=(1-At(i).^2).^2+4*(Vo.^2).*(At(i).^2); Sak(i)=1./sqrt(Sak(i)) ; %裸船体幅频响应 %Sek(i)=-2.0*Vo*At(i)/(1-At(i).^2); %Sek(i)=atan(Sek(i)); %裸船体相频响应 %%%%%%%%%%%%%%%%%%%%%%%%计算海浪谱函数 A=8.1e-3*9.8.^2; %(Mctr=1 为ITTC单参数谱) Ba=3.11/Hw.^2 ; Sw(i)=A*exp(-Ba/W(i)^4)./W(i)^5;%裸船体横摇运动统计计算 Slou(i)=(Sak(i)*W(i)^2/g)^2*Sw(i) ; %船体横摇谱函数 %Xtb= ; %横摇谱函数修正系数 %Slou(i)=(Sak(i)*W(i)^2/g)^2*Sw(i) ; end %plot(W,Sak,'r') Mlou=sum(Slou); Mlou=Mlou-(Slou(1)+Slou(100))/2.0; Dw=(W(100)-W(1))/100; Mlou=Dw*Mlou; Sf31=2.0*sqrt(Mlou)*57.3 %有义摇幅 Sfp=1.25*sqrt(Mlou)*57.3 %平均摇幅 Sfmax=2.55*sqrt(Mlou)*57.3 %最大摇幅 %%%%%%%%%%%%%%%%%%%%%%%%%% U型减摇水舱1设计参数 Lt1=5.5; %（一个水舱通道长度，沿ox轴方向） bt1=2; %（水舱边舱宽度） c1=13; %(水舱中心与边舱中心距离) H1=2.5; %（边舱的空气空间高度） ho1=2.5; %（边舱通道水侧高度） Ht1=5; hd1=0.55; %（水侧底部通道高度） Hz1=10; %（水舱布置高度） So1=Lt1*bt1; %（边舱截面积） S1=hd1*Lt1; %（水舱底部通道截面积） kg1=2.35; %（水舱通道截面积放大系数） Mt1=ho1+So1/S1*c1; %（水舱内水柱相当长度） Wt1=sqrt(g/Mt1) %（水舱1的频率） h1=ho1-hd1/2; %（水舱水平放置时，底部连线通道中心到边舱自由液面的垂直距离） Kcr1=-(Hz1-zg+hd1/2); %（船舶横摇中心到水舱底部连通道中心线的垂直距离,水舱在横摇中心以上时Kcr为负） b1=sqrt(2*(-Kcr1*c1+(ho1-hd1/2)*c1)); %(b平方的开平方) Jst1=p*So1*c1*b1^2; %（船舶与水舱1的耦合横摇惯性矩） Jt1=p*So1*c1^2*2*Mt1; %（水舱1内的液体的横摇质量惯性矩） fJt1=2*p*So1*c1^2*((S1/(So1*c1)+h1/c1^2)*Kcr1^2-h1^2/c1^2*Kcr1+(c1/3)*(S1/So1)+h1+h1^3/(3*c1^2));%（船舶由于计入水舱1而附加的横摇质量惯性矩） %fJt1=0; %%%%%%%%%%%%%%%%%%%%%%%%%% U型减摇水舱2设计参数 Wt2=0.44 %（取此频率为水舱2的固有频率） bt2=2; c2=13; %（船舶中心到水舱2边舱中心的距离） Ast=3.5; %（设水舱容量为3.5度） M=D*ho*Ast*3.14/180; %（所需的最大减摇力矩） Ht2=5; %（H2+ho2,水舱高） ho2=2.5; %(ho2=Ht2/2) So21=(M-p*g*So1*Ht1*c1)/(p*g*Ht2*c2);%（水舱2的截面积） Lt2=11; %(水舱2通道沿ox轴长度，还有待确定） So2=Lt2*bt2; S2=So2*c2/(g/Wt2^2-ho2); %（连通道截面积） hd2=S2/Lt2 Mt2=ho2+So2/S2*c2; %（水舱内水柱相当长度） Hz2=10; %（有待确定） h2=ho2-hd2/2; %（水舱水平放置时，底部连线通道中心到边舱自由液面的垂直距离） Kcr2=-(Hz2-zg+hd2/2); %（船舶横摇中心到水舱底部连通道中心线的垂直距离，此时Kcr2应为正） %Kcr2=15 b2=sqrt(2*(Kcr2*c2+(ho2-hd2/2)*c2)) ; %(取水舱2距基线2米) Jst2=p*So2*c2*b2^2; %（船舶与水舱2的耦合横摇惯性矩） Jt2=p*So2*c2^2*2*Mt2; %（水舱2内的液体的横摇质量惯性矩） fJt2=2*p*So2*c2^2*((S2/(So2*c2)+h2/c2^2)*Kcr2^2-h2^2/c2^2*Kcr2+c2/3*S2/So2+h2+h2^3/(3*c2^2));%（船舶由于计入水舱1而附加的横摇质量惯性矩） %fJt2=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(水舱1的计算参数） X1=Jst1/Jt1 %（与水舱1的位置布置有关的参数） Vt1=0.13 %（水舱1的阻尼系数） K1=2*p*g*So1*c1.^2/(I1+fJt1+fJt2)%(表征水舱1减摇能力的参数) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%（水舱2的计算参数） X2=Jst2/Jt2 %（与水舱2的位置布置有关的参数） Vt2=0.1 %（水舱2的阻尼系数） K2=2*p*g*So2*c2.^2/(I1+fJt1+fJt2)%(表征水舱2减摇能力的参数) %K2=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(计算计入水舱内液体质量的船舶横摇质量惯性矩后船舶的横摇频率） Kcr11=abs(Kcr1); KK1=abs(Kcr1-h1/2); Kcr22=abs(Kcr2); KK2=abs(Kcr2-h2/2); dKs=2*p*g*(c1*S1*Kcr11+h1*So1*KK1)+2*p*g*(c2*S2*Kcr22+h2*So2*KK2);%(由水舱内液体引起的横摇恢复系数增量） Ws1=sqrt((D*ho-dKs)/(I1+fJt1+fJt2)) %（计入水舱内液体质量的船舶横摇质量惯性矩后船舶的横摇频率） Ws=Wo; %%%%%%%%%%%%%%%%%%%%%（计算减摇水舱流体的重量） Mz=p*g*(2*bt1*h1+(B-2*bt1)*hd1)*Lt1+p*g*(2*bt2*h2+(B-2*bt2)*hd2)*Lt2;%（减摇水舱流体重量） dD=Mz/D %（减摇水舱流体重量占船舶排水量的百分比） %%%%%%%%%%%%%%%%%%%%%%（计算减摇水舱自由液面惯性矩和由此引起的横稳心高损失量，并计算安装减摇水舱后的船舶横稳心高） Mfs=p*g*Lt1/12*(B^3-(B-2*bt1)^3)+p*g*Lt2/12*(B^3-(B-2*bt2)^3) ; %（减摇水舱自由液面惯性矩） dh=Mfs/D %（横稳心高损失量） hh=ho-dh %（修正后的船舶横稳心高） for i=1:1:200 %(海浪遭遇频率) W(i)=i/200 ; A1(i)=(Wt1^2-W(i).^2)*Ws^2; B1(i)=2*Vt1*W(i).*Ws^2; C11(i)=W(i).^4-(Ws^2+Wt1^2+4*Vo*Vt1).*W(i).^2+Ws^2*Wt1^2; C12(i)=(-X1^2*K1/Wt1^2)*W(i).^4+2*X1*K1*W(i).^2-K1*Wt1^2; C1(i)=C11(i)+C12(i); D1(i)=-(2*Vo+2*Vt1)*W(i).^3+(2*Vo*Wt1^2+2*Vt1*Ws^2)*W(i); A2(i)=(Wt2^2-W(i).^2)*Ws^2; B2(i)=2*Vt2*W(i)*Ws^2; C21(i)=W(i).^4-(Ws^2+Wt2^2+4*Vo*Vt2)*W(i).^2+Ws^2*Wt2^2; C22(i)=(-X2^2*K2/Wt2^2)*W(i).^4+2*X2*K2*W(i).^2-K2*Wt2^2; C2(i)=C21(i)+C22(i); D2(i)=-(2*Vo+2*Vt2)*W(i).^3+(2*Vo*Wt2^2+2*Vt2*Ws^2)*W(i); M1(i)=A1(i).^2+B1(i).^2; M2(i)=(C1(i)+C22(i).*(C11(i).*C21(i)+D1(i).*D2(i))./(C21(i).^2+D2(i).^2)).^2; M3(i)=(D1(i)+C22(i).*(C21(i).*D1(i)-C11(i).*D2(i))./(C21(i).^2+D2(i).^2)).^2; M4(i)=sqrt(M1(i)./(M2(i)+M3(i))); %M(i)=sqrt((A2.^2*B2.^2)/((C2+C12*(C11*C21+D1*D2)/(C11.^2+D1.^2)).^2+(D2+C12*(C11*D2-C21*D1)/(C11.^2+D1.^2)).^2)); %计算海浪谱函数 A=8.1e-3*9.8^2; %(Mctr=1 为ITTC单参数谱) Ba=3.11/Hw^2 ; Sw(i)=A*exp(-Ba/W(i)^4)/W(i)^5; Sww(i)=(M4(i)*W(i)^2/g)^2*Sw(i);%双水舱之船体横摇谱函数 end Sssc=sum(Sww); Sssc=Sssc-(Sww(1)+Sww(100))/2.0; Dw=(W(100)-W(1))/100; Sssc=Dw*Sssc; F31=2.0*sqrt(Sssc)*57.3 %有水舱之船体横摇运动统计计算 Fp=1.25*sqrt(Sssc)*57.3 Fmax=2.55*sqrt(Sssc)*57.3 M4=M4*0.9; figure(1) plot(W,M4,'r',W,Sak) grid %%%双水舱系统的传递函数%%%%%%%% Vs=Vo; Num=[Ws^2 Ws^2*(2*Vt1+2*Vt2) Ws^2*(Wt1^2+4*Vt1*Vt2+Wt2^2) Ws^2*(2*Vt1*Wt2^2+2*Vt2*Wt1^2) Ws^2*Wt1^2*Wt2^2]; Den=[(1-K1*X1^2/Wt1^2-K2*X2^2/Wt2^2) (2*Vt1+2*Vt2+2*Vs-2*K1*Vt2*X1^2/Wt1^2-2*K2*Vt1*X2^2/Wt2^2)... (Wt1^2+Wt2^2+Ws^2+4*Vt1*Vt2+4*Vt1*Vs+4*Vt2*Vs-X1^2*Wt2^2*K1/Wt1^2-2*K1*X1-X2^2*Wt1^2*K2/Wt2^2-2*K2*X2)... (2*Vt1*Wt2^2 +2*Vs*Wt2^2+2*Vt2*Wt1^2+4*Vs*Vt1*Vt2+2*Vt2*Ws^2+2*Vs*Wt1^2+2*Vt1*Ws^2-4*X1*Vt2*K1-4*X2*Vt1*K2)... (Wt1^2*Wt2^2+4*Wt2^2*Vs*Vt1+Ws^2*Wt2^2+4*Vs*Vt2*Wt1^2+4*Vt1*Vt2*Ws^2+Ws^2*Wt1^2-2*K1*X1*Wt2^2-K1*Wt1^2-2*K2*X2*Wt1^2-K2*Wt2^2)... (2*Vs*Wt1^2*Wt2^2+2*Vt1*Ws^2*Wt2^2+2*Vt2*Ws^2*Wt1^2-2*K1*Vt2*Wt1^2-2*K2*Vt1*Wt2^2)... (Ws^2*Wt1^2*Wt2^2-Wt1^2*Wt2^2*K1-Wt2^2*K2*Wt1^2)]; fprintf('双水舱系统的传递函数：') Doubletanktf=tf(Num,Den)