《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 263|回复: 0
收起左侧

[提问] 三自由度剪力结构,采用瑞利阻尼,强迫振动分析的代码。哪个大神帮忙看看哪里错了

[复制链接]
发表于 2018-5-15 17:14:44 | 显示全部楼层 |阅读模式
m=[180000,270000,270000];
m=diag(m);
k=[98000000 -98000000 0;-98000000 294000000 -196000000;0 -196000000 441000000];
c=0.9308*m+0.002294*k;
f0=2000000;
theta=250;
t1=5;
tm=20;
dt=0.01;
nt=tm/dt;
T=[0:dt:tm];
alfa=0.5;
beta=0.25;
a0=1/alfa/dt/dt;
a1=beta/alfa/dt;
a2=1/alfa/dt;
a3=1/2/alfa-1;
a4=beta/alfa-1;
a5=dt/2*(beta/alfa-2);
a6=dt*(1-beta);
a7=dt*beta;
d=zeros(3,nt);
v=zeros(3,nt);
a=zeros(3,nt);
for i=2:nt
t=(i-1)*dt;
if (t<t1), f=[0;f0*sin(theta*t);0]; else f=[0;0;0];
end
ke=k+a0*m+a1*c;
fe=f+m*(a0*d(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+c*(a1*d(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));
d(:,i)=inv(ke)*fe;
a(:,i)=a0*(d(:,i)-d(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);
v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);
end
end
close all
figure
plot(T,d(1,:)*1000 ,'--b')
ylabel('{\itd}_1/mm')
figure
plot(T,d(2,:)*1000,'--b')
ylabel('{\itd}_2/mm')
figure
plot(T,d(3,:)*1000,'--b')
ylabel('{\itd}_3/mm')
figure
plot(T,v(1,:),'--g')
ylabel('{\itv}_1/m/s')
figure
plot(T,v(2,:),'--g')
ylabel('{\itv}_2/m/s')
figure
plot(T,v(3,:),'--g')
ylabel('{\itv}_3/m/s')
figure
plot(T,a(1,:),'--r')
ylabel('{\ita}_1/m^2/s')
figure
plot(T,a(2,:),'--r')
ylabel('{\ita}_2/m^2/s'
figure
plot(T,a(3,:),'--r')
ylabel('{\ita}_3/m^2/s')
您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 陕ICP备08102094号 

GMT+8, 2018-10-19 01:34 , Processed in 0.078088 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表