%这是精确解程序
clear all
close all
m=0.5;
c=0.05;
k=1;
x0=0;
v0=0;
tf=30;
delt=0.001;
w=5;
f0=2;
wn=sqrt(k/m);
z=c/(2*m*wn);
lan=w/wn;
wd=wn*sqrt(1-z^2);
A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);
t=0:tf/1000:tf;
phi=atan2(2*z*lan,1-lan^2);
phi1=atan2(v0+z*wn*x0,x0*wd);
B=f0/(k*sqrt((1-lan^2)^2+(2*lan*z)^2));
x=A*exp(-z*wn*t).*sin(wd*t+phi1)+B*sin(w*t+phi);
plot(t,x);grid
xlabel('时间(s)')
ylabel('位移')
title('位移与时间的关系')
%这是Newmark法的程序
close all
clear all
tf=30;
delt=0.01;
fid1=fopen('disp1','wt');
m=0.5;
c=0.05;
k=1;
x0=0;
v0=0;
bita=1/6;
md=inv(m+delt/2*c+bita*delt^2*k);
ml=inv(m);
for t=0:delt:tf
f=2*sin(5*t);
if t==0
xdd0=ml*(f-k*x0-c*v0);
else
xdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0));
xd=v0+delt/2*(xdd0+xdd);
x=x0+delt*v0+(1/2-bita)*delt^2*xdd0+bita*delt^2*xdd;
xdd0=xdd;
v0=xd;
x0=x;
fprintf(fid1,'%10.4f',x0);
end
end
fid2=fopen('disp1','rt');
n=tf/delt;
x=fscanf(fid2,'%f');
t=1:n;
figure('numbertitle','off','name','自由度1的位移','pos',[450 180 400 420]);
plot(t,x),grid on,xlabel('时间*0.1秒'),title('自由度1的位移与时间的关系'); |