% newmark方法的程序实现
% 算例为书《结构动力学基础》中单自由度系统的强迫振动问题，例2-11中的结果，与计算结果一致
% ddx+2w*0.05dx+w2x=f1/m
clc
clear
%*******************生成刚度阵、质量阵、阻尼阵********************
kk=4*pi^2;
mm=1;
w=sqrt(kk/mm);
m=[1]; %质量矩阵
c=[2*w*0.05]; %阻尼矩阵
k=[w^2]; %刚度矩阵
f1=[0]; %激励力向量
%******************初始位移、速度*****************
d0=[0.2];
d1=[0];
%******************计算DDX0************************
d2=inv(m)*(f1-c*d1-k*d0);
%******************计算积分常数********************
b=0.25; %参数a
r=0.5; %DETA>=0.5
dt=1/100; %detaT
a0=1/(b*(dt)^2); %积分常数
a1=r/(b*dt);
a2=1/(b*dt);
a3=1/(2*b)-1;
a4=r/b-1;
a5=0.5*dt*((r/b)-2);
a6=dt*(1-r);
a7=r*dt; %积分常数
%******************计算等效刚度*******************
pk=k+a0*m+a1*c;
pk=inv(pk);
%******************每一时间部计算*****************
for i=1:900
t(i)=dt*i; %时间间隔
ff=10*sin(pi*t(i))/mm; %t(k+1)时刻载荷已知
ff=ff+m*(a0*d0+a2*d1+a3*d2)+c*(a1*d0+a4*d1+a5*d2); %t(k+1)时刻等效力
d00=pk*ff; %x(k+1)
d22=a0*(d00-d0)-a2*d1-a3*d2; %ddx(k+1)
d11=d1+a6*d2+a7*d22; %dx(k+1)
d0=d00;
d1=d11;
d2=d22;
a(i)=d00(1);
end
plot(t,a)
评论0