单自由度简谐振动算例

西风独自凉

例子：见张亚辉、林家浩《结构动力学基础》P34例2-11
解析解P21公式2-77
function vtb22(m,c,k,x0,v0,tf,w,f0)
%单自由度系统的谐迫振动
wn=sqrt(k/m); %频率
z=c/2/m/wn; %阻尼比
wd=wn*sqrt(1-z^2);
e=z*wn;
A1=-2*f0*e*w/(m*((w^2-wn^2)^2+4*e^2*w^2));
A2=f0*(wn^2-w^2)/(m*((w^2-wn^2)^2+4*e^2*w^2));
t=0:tf/1000:tf;
x=exp(-z*wn*t).*((v0+e*x0-e*A1)/wd*sin(wd*t)+x0*cos(wd*t)-A2*w*sin(wd*t)/wd-A1*cos(wd*t))+A1*cos(w*t)+A2*sin(w*t);
plot(t,x)
grid

输入：vtb22(1,0.2*pi,4*pi*pi,0.2,0,10,pi,10)
纽马克法程序1
%ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
function vtb300(tf,delt)
%tf为持续时间，delt为时间步长
fid1=fopen('disp','wt');
m=[1];
c=[0.2*pi ];
k=[4*pi*pi];
x0=[0.2]';
v0=[0]';
bita=1/6;
md=inv(m+delt/2*c+bita*delt^2*k);
m1=inv(m);
[E,F]=eig(m1*k);
diag(sqrt(F));
for t=0:delt:tf
     f=[10*sin(pi*t)]';
     if t==0; xdd0=m1*(f-k*x0-c*v0);
     else
     xdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0));%加速度
     x=md*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);
     xd=v0+delt/2*(xdd0+xdd);%速度
     xdd0=xdd; v0=xd;x0=x;
     fprintf(fid1,'%10.4f',x0);
     t;
     end
end
fid2=fopen('disp','rt');
n=tf/delt;
x=fscanf(fid2,'%f',[1,n]);
x1=x(1,:)';
t=delt:delt:tf;
plot(t,x1,'-')
xlabel('t[sec]')
ylabel('x(m)')


输入：vtb300(10,0.01)

纽马克法程序2
function newmark1(t_max,dt)
%t_max为持续时间，dt为时间步长
M=[1];
C=[0.2*pi];
K=[4*pi*pi];
u=[0.2]';
v=[0]';
a=[0]';
t(1)=0;               %时间
x(:,1)=u;             %位移
x1(:,1)=v;            %速度
x2(:,1)=a;            %加速度
%newmark参数
gama=0.5;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt;
%等效刚度矩阵
K1=K+a0*M+a1*C;
i=1;
t(1)=0;
while t(i)<t_max
      f=10*sin(pi*t(i));
      Q=[f]';
      q(:,i+1)=Q+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i));
      x(:,i+1)=inv(K1)*q(:,i+1);
      x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
      x1(:,i+1)=a1*(x(:,i+1)-x(:,i))-a4*x1(:,i)-a5*x2(:,i);
      i=i+1;
      t(i)=t(i-1)+dt;
end  

x1=x(1:t_max/dt)';
t=dt:dt:t_max;
plot(t,x1,'-')
xlabel('t[sec]')
ylabel('x(m)')
grid
输入：newmark1(10,0.01)

小星星

谢谢楼主！！！！！！