Course Stability Codes for Ship Towing System

clear
load coefficient.dat %towing boat parameter
load bargeskeg.dat   %towed boat para
rho=998;
U=3.6; % towing speed

prope1=coefficient;
L1=40;%length of towed body
d1=2.2;%draft


ellT=-0.2*L1;%lt distance between com and towing point (on tug)
m1=494.7*998;% tug mass
%
Xuu1=prope1(1)*(0.5*rho*L1*d1);
Yv1=prope1(2)*(0.5*rho*L1*d1*U);
Yr1=prope1(3)*(0.5*rho*L1*L1*d1*U);
Nv1=prope1(4)*(0.5*rho*L1*L1*d1*U);
Nr1=prope1(5)*(0.5*rho*L1*L1*L1*d1*U);
mx1=prope1(6)*(0.5*rho*L1*L1*d1);
my1=prope1(7)*(0.5*rho*L1*L1*d1);
Iz1=0.02*(0.5*rho*L1*L1*L1*L1*d1);
Jz1=prope1(8)*(0.5*rho*L1*L1*L1*L1*d1);
Yd=prope1(9)*(0.5*rho*L1*d1*U*U);
Nd=prope1(10)*(0.5*rho*L1*d1*U*U);

%towed boat
L2=60.96;%length
d2=2.74;%draft
prope2=bargeskeg;
nm2=prope2(1);
nmx2=prope2(2);
nmy2=prope2(3);
nIz2=prope2(4);
nJz2=prope2(5);
nXuu2=prope2(6);
nYv2=prope2(7);
nYr2=prope2(8);
nNv2=prope2(9);
nNr2=prope2(10);

%nondimension-->dimension
m2=prope2(1)*(0.5*rho*L2*L2*d2);
mx2=prope2(2)*(0.5*rho*L2*L2*d2);
my2=prope2(3)*(0.5*rho*L2*L2*d2);
Iz2=prope2(4)*(0.5*rho*L2*L2*L2*L2*d2);
Jz2=prope2(5)*(0.5*rho*L2*L2*L2*L2*d2);
Xuu2=prope2(6)*(0.5*rho*L2*d2);
Yv2=prope2(7)*(0.5*rho*L2*d2*U);
Yr2=prope2(8)*(0.5*rho*L2*L2*d2*U);
Nv2=prope2(9)*(0.5*rho*L2*L2*d2*U);
Nr2=prope2(10)*(0.5*rho*L2*L2*L2*d2*U);

towRope=linspace(0.2,5,500);% nondimensional
comPos=linspace(0.1,1,length(towRope)); % nondimensional
CSbody=zeros(length(towRope),length(comPos));
nCSbody=CSbody;
CSsystem=CSbody;
for i=1:length(towRope)
    for j=1:length(comPos)
        nell1=towRope(i);
        nell2=comPos(j);
        ell1=towRope(i)*L2;
        ell2=comPos(j)*L2;
        
% only considering dynamics of towed body 
        a41=-ell2*U^2*Xuu2/(Iz2+Jz2);
        a42=-Nv2*ell1/(Iz2+Jz2);
        a43=(ell2*U^2*Xuu2-U*Nv2)/(Iz2+Jz2);
        a44=(-Nv2*ell2+Nr2)/(Iz2+Jz2);
        tmp=-ell2/ell1;
        a21=tmp*a41+ U^2*Xuu2/ell1/(m2+my2);
        a22=tmp*a42+ Yv2/(m2+my2);
        a23=tmp*a43+ (-Xuu2*U^2+U*Yv2)/ell1/(m2+my2);
        a24=tmp*a44+ (Yv2*ell2-Yr2-U*(my2-mx2))/ell1/(m2+my2);

        
        A=[0 1 0 0;
           a21 a22 a23 a24;
            0 0 0 1;
            a41 a42 a43 a44];
        [T,AA]=balance(A);
        if(all(real(eig(AA))<0))
            CSbody(i,j)=1;
        else
            CSbody(i,j)=0;
        end
        
        % nondimen analysis
D0=nXuu2*(nell2*nYv2-nNv2);
        D1=nXuu2*(nNr2+(nell1+nell2)*nell2*nYv2-nell1*nNv2-nell2*(nYr2-nmx2+nmy2+nNv2));        D2=nell1*((nm2+nmx2-nYr2)*nNv2+nYv2*nNr2)-nXuu2*(nell2*(nm2+nmy2)*(nell1+nell2)+nIz2+nJz2);
        D3=-nell1*((nm2+nmy2)*nNr2+(nIz2+nJz2)*nYv2);
        D4=nell1*(nm2+nmy2)*(nIz2+nJz2);
        D5=D3*D2*D1-D3^2*D0-D4*D1^2;
        vec=[D0;D1;D2;D3;D4;D5];
        if(all(vec>0))
            nCSbody(i,j)=1;
        end
        
    

% one tug and one towed body

        b11=Yv1/(m1+my1);
        b12=(Yr1-(m1+mx1)*U)/(m1+mx1);
        b13=0;
        b14=0;
        b15=-Xuu2*U^2/(m1+my1);
        b16=Xuu2*U^2/(m1+my1);
        b17=0;
        b21=Nv1/(Iz1+Jz1);
        b22=Nr1/(Iz1+Jz1);
        b23=0;
        b24=0;
        b25=-Xuu2*ellT*U^2/(Iz1+Jz1);
        b26=Xuu2*ellT*U^2/(Iz1+Jz1);
        b27=0;
        b41=Nv2/(Iz2+Jz2);
        b42=Nv2*ellT/(Iz2+Jz2);
        b43=-Nv2*ell1/(Iz2+Jz2);
        b44=(Nr2-Nv2*ell2)/(Iz2+Jz2);
        b45=U*Nv2/(Iz2+Jz2);
        b46=-ell2*Xuu2*U^2/(Iz2+Jz2);
        b47=ell2*Xuu2*U^2/(Iz2+Jz2)-Nv2*U/(Iz2+Jz2);

        tmp1=-ell2/ell1;
        tmp2=ellT/ell1;
        tmp3=1/ell1;
        b31=tmp1*b41+tmp2*b21+tmp3*b11+Yv2/(-ell1*(m2+my2));
        b32=tmp1*b42+tmp2*b22+tmp3*b12+(Yv2*ellT-(m2+my2)*U)/(-ell1*(m2+my2));
        b33=tmp1*b43+tmp2*b23+tmp3*b13+(-Yv2*ell1)/(-ell1*(m2+my2));
        b34=tmp1*b44+tmp2*b24+tmp3*b14+(Yr2-Yv2*ell2+(my2-mx2)*U)/(-ell1*(m2+my2));
        b35=tmp1*b45+tmp2*b25+tmp3*b15+Yv2*U/(-ell1*(m2+my2));
        b36=tmp1*b46+tmp2*b26+tmp3*b16+(-Xuu2*U^2)/(-ell1*(m2+my2));
        b37=tmp1*b47+tmp2*b27+tmp3*b17+(Xuu2*U^2-Yv2*U)/(-ell1*(m2+my2));

        b51=0;b52=1;b53=0;b54=0;b55=0;b56=0;b57=0;
        b61=0;b62=0;b63=1;b64=0;b65=0;b66=0;b67=0;
        b71=0;b72=0;b73=0;b74=1;b75=0;b76=0;b77=0;
        B=[b11 b12 b13 b14 b15 b16 b17;
            b21 b22 b23 b24 b25 b26 b27;
            b31 b32 b33 b34 b35 b36 b37;
            b41 b42 b43 b44 b45 b46 b47;
            b51 b52 b53 b54 b55 b56 b57;
            b61 b62 b63 b64 b65 b66 b67;
            b71 b72 b73 b74 b75 b76 b77];
%         real(eig(B))
        [T,BB]=balance(B,'noperm');
        if(all(real(eig(BB))<1e-6))
            CSsystem(i,j)=1;
        else
            CSsystem(i,j)=0;
        end
    end
end

figure(1);
contourf(towRope,comPos,CSbody);
xlabel('l1');
ylabel('l2');
colormap('gray');
title('tug fixed');


% figure(2);
% contourf(towRope,comPos,nCSbody);
% xlabel('l1');
% ylabel('l2');
% colormap('gray');


figure(3);
contourf(towRope,comPos,CSsystem);
xlabel('l1');
ylabel('l2');
colormap('gray');
title('tug-ship connection');

Course stability results when tug is fixed.
Course stability results when considering the flexibility of both tug and ship.

The coefficients used in the code can be download from here.

Leave a Comment

Your email address will not be published. Required fields are marked *