MVH 发表于 2005-7-28 17:46

maker系列-有阻尼自由振动程序

一个简单的例子,可以对照相关振动和动力学的书,非常容易理解,望对初学者有一定的帮助,针对临界阻尼,欠阻尼,过阻尼的情况

function maker1(m,c,k,x0,v0,time)
%maker1用来计算单自由度系统响应
%此程序针对两种输入形势(m,c,k,x0,v0,tf)和(zeta,wn,x0,v0,tf)
%m-为质量,c-阻尼,k-刚度,x0-初位移,v0-除速度;time-仿真时间
%zeta-阻尼系数,wn-固有频率
%chinamaker@dytrol.com
%2003.11.26
clc
%针对输入形势做出判断
if nargin==5 %输出参数的个数
z=m;
wn=c;
time=v0;
v0=x0;
x0=k;
m==1;
c=2*z*wn;
k=wn^2;
end

k=1;
wn=sqrt(k/m);%固有频率
z=c/2/m/wn;%相对阻尼系数
wd=wn*sqrt(1-z^2);%系统的实际振动频率
fprintf('固有频率%.3g.rad/s.\n',wn);
fprintf('阻尼系数%.3g.\n',z);
fprintf('振动频率%.3g.\n',wd);
t=0:time/1000:time;
if z<1 %欠阻尼情况
Ampli=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);%系统振幅
phi=atan2(x0*wd,v0+z*wn*x0);%系统相位角
x=Ampli*exp(-z*wn*t).*sin(wd*t+phi);%系统振动时程曲线
fprintf('Ampli=%.3g\n',Ampli);
fprintf('Ampli2=%,3g\n',phi);
elseif z==1 %临界阻尼
a1=x0;
a2=v0+wn*x0;
fprintf('a1=%.3g.\n',a1);
fprintf('a2=%.3g.\n',a2);
x=(a1+a2*t).*exp(-wn*t);
else %过阻尼
a1=(-v0+(-zsqrt(z^2-1))*wn*x0)*wn*x0/2/wn/sqrt(z^2-1);
a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);
fprintf('a1=%.3g.\n',a1);
fprintf('a2=%.3g.\n',a2);
x=exp(-wn*t).*(a1*exp(wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t));
end
plot(t,x);
grid;
xlabel('Time(s)');
ylabel('Amplitude');

运行结果:
>>maker1(1,0.05,1,1,1,100)
固有频率1.rad/s.
阻尼系数0.025.
振动频率1.
Ampli=1.43
phi=0.773

vib 发表于 2006-12-7 20:40

输入形式是(zeta,wn,x0,v0,tf)是什莫含义?程序里对应法nargin==5的else在哪里?

suffer 发表于 2006-12-9 09:43

输入形式是(zeta,wn,x0,v0,tf)是什莫含义?

注释中不是有说明吗?

suffer 发表于 2006-12-9 09:45

程序里对应法nargin==5的else在哪里?

不需要else,如果输入参数不是5个,就不需要做特殊处理,直接跳过

if nargin==5 %输出参数的个数
z=m;
wn=c;
time=v0;
v0=x0;
x0=k;
m==1;
c=2*z*wn;
k=wn^2;
end

这一段
页: [1]
查看完整版本: maker系列-有阻尼自由振动程序