MATLAB+鱼雷航行力学.doc

clear all clc global B G M XC YC ZC Jx Jy Jz Lu S L T global Cx global Mxbeta Mxdr Mxdd Mxwx Mxwy Mxp global Cyalfa Cyde Cywz global Mzalfa Mzde Mzwz global Czbeta Czdr Czwy global Mybeta Mydr Mywy global K11 K22 K26 K33 K44 K55 K66 K35 %输入鱼雷流体动力因数 Cx=0.141; Mxbeta=0.00152;Mxdr=-0.000319;Mxdd=-0.0812;Mxwx=-0.0044;Mxwy=0.0008;Mxp=0; Cyalfa=2.32;Cyde=0.51;Cywz=1.17; Mybeta=0.69;Mydr=-0.11;Mywx=0;Mywy=-0.61; Czbeta=-2.32;Czdr=-0.2;Czwy=-1.17; Mzalfa=0.69;Mzde=-0.28;Mzwz=-0.61; K11=0.0222;K22=1.1096;K33=K22;K44=0.1406; K55=3.8129;K66=K55;K26=-0.363;K35=-K26; % 输入鱼雷总体参数 Lu=1019.2; %kg/m^3密度 M=1840; %kg G=M*9.8; %n L=7.738; %m B=14671; %n浮力 S=0.224; %m^2 XC=0.02; %m YC=-0.005; %m ZC=0.008; %m T=10672; %n推力 Jx=289.3; %kg*m^2 Jy=6771.8; %kg*m^2 Jz=Jy; %kg*m^2 % 计算量纲附加质量 V=B/(Lu*9.8); %m^3 L11=K11*Lu*V; %kg L22=K22*Lu*V; %kg L33=K33*Lu*V; %kg L44=K44*Lu*V^(5/3); %kg*m^2 L55=K55*Lu*V^(5/3); %kg*m^2 L66=K66*Lu*V^(5/3); %kg*m^2 L26=K26*Lu*V^(4/3); %kg*m L35=K35*Lu*V^(4/3); %kg*m % 下面计算纵向运动的操纵性 % 计算平衡攻角和平衡舵角 deltaG=G-B; X_C=XC/L;Y_C=YC/L;Z_C=ZC/L; delta_G=deltaG/G;mu=(2*M)/(Lu*S*L); v0=sqrt(2*T/(Lu*S*Cx)); alfa0=(deltaG*Mzalfa-G*X_C*Cyde)/(Cyalfa*Mzde-Cyde*Mzalfa)*(2/(Lu*S*v0^2)); de0=(G*X_C*Cyalfa-deltaG*Mzalfa)/(Cyalfa*Mzde-Cyde*Mzalfa)*(2/(Lu*S*v0^2)); % 计算纵向无量纲位置力臂、旋转力臂与稳定性裕度 L_alfa=Mzalfa/Cyalfa; L_wz=(mu*X_C-Mzwz)/(mu-Cywz); Gwz=1-L_alfa/L_wz; % 计算纵向扰动方程的因数 a21=(M*XC+L26)/(0.5*Lu*v0^2*S);a22=(mu-Cywz)*L/v0;a23=-(M+L22)/(0.5*Lu*v0*S); a24=-Cyalfa;a25=Cyde; a31=(Jz+L66)/(0.5*Lu*v0^2*S*L);a32=(mu*X_C-Mzwz)*L/v0; a33=-(M*XC+L26)/(0.5*Lu*v0*S*L);a34=-Mzalfa;a35=Mzde; % 计算纵向运动扰动方程的状态空间模型 h=a23*a31-a21*a33; A1=[(a21*a34-a31*a24)/h,(a21*a32-a31*a22)/h,0,0; (a33*a24-a23*a34)/h,(a33*a22-a23*a32)/h,0,0; 0,1,0,0; -v0,0,v0,0]; B1=[(a31*a25-a21*a35)/h;(a23*a35-a33*a25)/h;0;0]; %由纵向状态空间模型计算纵向各参数的传递函数 [num_alfa,den_alfa]=ss2tf(A1,B1,[1,0,0,0],0); [num_wz,den_wz]=ss2tf(A1,B1,[0,1,0,0],0); [num_sita,den_sita]=ss2tf(A1,B1,[0,0,1,0],0); [num_y0,den_y0]=ss2tf(