传递矩阵-matlab程序

%main_critical.m

%该程序使用Riccati传递距阵法计算转子系统的临界转速及振型

%本函数中均采用国际单位制

% 第一步：设置初始条件（调用函数shaft_parameters）

%初始值设置包括：轴段数N，搜索次数M

%输入轴段参数：内径d,外径D，轴段长度l，支撑刚度K，单元质量mm，极转动惯量Jpp[N,M,d,D,l,K,mm,Jpp]=shaft_parameters;

% 第二步：计算单元的5个特征值（调用函数shaft_pra_cal）

%单元的5个特征值：

%m_k：:质量

%Jp_k：极转动惯量

%Jd_k：直径转动惯量

%EI：弹性模量与截面对中性轴的惯性矩的乘积

%rr：剪切影响系数

[m_k,Jp_k,EI,rr]=shaft_pra_cal(N,D,d,l,Jpp,mm);

% 第三步：计算剩余量（调用函数surplus_calculate）,并绘制剩余量图

%剩余量：D1

for i=1:1:M

ptx(i)=0;

pty(i)=0;

end

for ii=1:1:M

wi=ii/1*2+50;

[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,JD_k,l,EI,rr);

D1;

pty(ii)=D1;

ptx(ii)=w1

end

ylabel(‘剩余量’)；

plot(ptx,pty)

xlabel(‘角速度red/s’);

grid on

% 第四步：用二分法求固有频率及振型图

%固有频率：Critical_speed

wi=50;

for i=1:1:4

order=i

[D1,SS,Sn]=surplus_calculate(N,wi,k,m_k,Jp_k,Jd_k,l,EI,rr);

Step=1;

D2=D1;

kkk=1;

while kkk<5000

if D2*D1>0

wi=wi+step;

D2=D1;

[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);

end

if D1*D2<0

wi=wi-step;

step=step/2;

wi=wi+step;

[D1,SS,Sn] =surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);

End

D1;

Wi;

If atep<1/2000

Kkk=5000;

end

end

Critical_speed=wi/2/pi*60

figure;

plot_mode(N,l,SS,Sn)

wi=wi+2;

end

%surplus_calculate,.m

%计算剩余量

%（1）计算传递矩阵

%（2）计算剩余量

function [D1,SS,Sn]= surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); % (1)计算传递矩阵

%＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

%(a)初值设为0

%＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

for i=1:1:N+1

for j=1:1:2

for k=1:1:2

ud11(j,k.i)=0;

ud12(j,k.i)=0;

ud21(j,k.i)=0;

ud22(j,k.i)=0;

end

end

end

for i=1:1:N

for j=1:1:2

for k=1:1:2

us11(j,k.i)=0;

us12(j,k.i)=0;

us21(j,k.i)=0;

us22(j,k.i)=0;

end

end

end

for i=1:1:N

for j=1:1:2

for k=1:1:2

u11(j,k.i)=0;

u12(j,k.i)=0;

u21(j,k.i)=0;

u22(j,k.i)=0;

end

end

end

%＝＝＝＝＝＝＝＝＝＝＝＝

%(b)计算质点上传递矩阵―――点矩阵的一部分！

%＝＝＝＝＝＝＝＝＝＝＝＝

for i=1:1:N+1

ud11(1,1,i)=1; ud11(1,2,i)=0; ud11 (2,1,i)=0; ud11(2,2,i)=1;

ud21(1,1,i)=0; ud21(1,2,i)=0; ud21 (2,1,i)=0; ud21(2,2,i)=0;

ud22(1,1,i)=1; ud22(1,2,i)=0; ud22 (2,1,i)=0; ud22(2,2,i)=1;

end

%＝＝＝＝＝＝＝＝＝＝＝＝

%(c)计算质点上传递矩阵―――点矩阵的一部分！

%＝＝＝＝＝＝＝＝＝＝＝＝

for i=1:1:N+1

ud12(1,1,i)=0;

ud12(1,2,i)=(Jp_k(i)-Jd_k(i))*wi^2; %%%考虑陀螺力矩

ud12(2,1,i)=m_k(i)*wi^2-k(i);

ud12(2,2,i)=0;

end

%＝＝＝＝＝＝＝＝＝＝＝＝

%(d)以下计算的是无质量梁上的传递矩阵―――场矩阵

%计算的锥轴的us是不对的，是随便令的，在后面计算剩余量时，zhui中会把错误的覆盖掉

%＝＝＝＝＝＝＝＝＝＝＝＝

for i=1:1:N

us11(1,1,i)=1; us11(1,2,i)=1(i); us11 (2,1,i)=0; us11(2,2,i)=1;

us12(1,1,i)=0; us12(1,2,i)=0; us12 (2,1,i)=0; us12(2,2,i)=0;

us21(1,1,i)=1(i)^2/(2*EI(i)); us21(1,2,i)=(1(i)^3*(1-rr(i))/(6*EI(i)); us21(2,1,i)=1(i)/EI(i);

us21(2,2,i)=1(i)^2/(2*EI(I));

us22(1,1,i)=1; us22(1,2,i)=1(i); us22 (2,1,i)=0; us22(2,2,i)=1;

end