马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
最近在学习传递矩阵法,按照传递矩阵理论推导,f1=S1*e1,其中f1=0,e1≠0,于是必有|S1|=0;前阵子在论坛里面看到其他人写的程序(如下),里面提到初值S=[0,0;0,0],而|S1|=0和S=[0,0;0,0]还是有区别的,想请教一下各位,这里S矩阵初值必须是零矩阵吗?
clear;
l=[1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,0]; %13个结点
m=[2940,5880,5880,5880,5880,5880,5880,5880,5880,5880,5880,5880,2940];
Jp=[0,0,0,0,0,0,0,0,0,0,0,0,0];
Jd=[0,0,0,0,0,0,0,0,0,0,0,0,0]; %不计转动惯量和陀螺力矩
I=[100,100,100,100,100,100,100,100,100,100,100,100,100];
E=4393;
v=[0,0,0,0,0,0,0,0,0,0,0,0,0]; %不计剪切影响
k=[1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9]; %油膜刚度
kb=[2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9]; %轴承座刚度
mb=[3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577]; %参振质量
S=[0,0;0,0]; %Riccati第一矩阵
s=1;
x=[]; %用于记录剩余量
for n=1864:1864 %试算频率
for i=1:13
K=k(i)*(kb(i)-mb(i)*n^2)/(k(i)+kb(i)-mb(i)*n^2) %总刚度
u11=[1,l(i);0,1];
u12=[l(i)*(m(i)*n^2-K),(Jp(i)-Jd(i))*n^2;m(i)*n^2-K,0];
u21=(l(i)/(E*I(i))).*[l(i)/2,l(i)^2*(1-v(i))/6;1,l(i)/2];
u22=[1+l(i)^3*(1-v(i))*(m(i)*n^2-K)/(6*E*I(i)),l(i)+l(i)^2*(Jp(i)-Jd(i))*n^2/(2*E*I(i));l(i)^2*(m(i)*n^2-K)/(2*E*I(i)),1+l(i)*(Jp(i)-Jd(i))*n^2/(E*I(i))];
S=[u11*S+u12]*inv([u21*S+u22])
end
D=det(S);
x=[x,D]; %记录剩余量
end
n=1864:1:1864; %产生曲线的横坐标
grid on
plot(n,x)
|