中心差分法计算单自由度体系动力反应样本.doc

资料内容仅供您学习参考，如有不当或者侵权，请联系改正或者删除。
中心差分法计算单自由度体系动力反应

1, 程序说明
中心差分法基于有限元查分代替位移对时间的求导( 即速度和加速度) 。如果采用等时间步长, , 则速度和加速度的中心查分近似为

体系的运动方程为

联立以上三式, 得

上式中, 假设和是已知的, 即和以前时刻的运动已知, 则能够把已知项移到方程的右边, 整理得

这样, 就能够计算体系任意时刻的位移, 速度和加速度。

2, 程序框图

3, 程序清单
%计算等效刚度和中心差分法计算公式中的系数
clear,
m=17.5e3;k=875500;c=35000;aa=input(' 请选择时间步长 1 or 2 or 3 \n 1: dt=0.02 ; 2 : dt=0.3 ; 3: dt=其它\n');
if aa==1
dt=0.02;
end
if aa==2
dt=0.3;
end
if aa==3
dt=input('请输入时间步长\n dt= ')
end
t=0:dt:1.2;
n=fix(1.2/dt+1);kr=m/(dt * dt) + c/(2 * dt); a=k-2 * m/(dt * dt); b=m/(dt * dt)-c/(2 * dt);
%求力p
p1=0:1:40;p2=39:-1:0;one=ones(1,40);p3=(one<0);p=1000*[p1,p2,p3];
for i=1:n
if t(i)<=0.4,p(i)=100000*t(i);
end
if t(i)>0.4&&t(i)<=0.8,p(i)=8*t(i);
end
if t(i)>0.8,p(i)=0;
end
end

%以下求位移u

u(1,1)=0;u(1,2)=0;
for i=2:(n-1)
u(1,i+1) =( p(1,i) - a * u(1,i) - b * u(1,i-1))/kr;
end
%以下求速度ui, 加速度uii

ui(1,1)=0;uii(1,1)=0;
for i=2:(n-1)
ui(1,i)=(u(1,i+1)-u(1,i-1))/(2 * dt);
uii(1,i)=(u(1,i+1)-2 * u(1,i)+u(1,i-1))/(dt*dt);
end
ui(1,n)=2*ui(1,n-1)-ui(1,n-2);uii(1,n)=2*uii(1,n-1)-uii(1,n-2);
figure(1)
plot(t,u,'b'),xlabel('时间( s) '),ylabel('位移( m) '),title('位移时程曲线 ');
figure(2)
plot(t,ui,'r'),xlabel('时间( s) '),ylabel('速度( m/s) '),title('速度时程曲线 ');
figure(3)
plot(t,uii,'k'),xlabel('时间( s) '),ylabel('加速度(m/s^2)'),title('加速度时程曲线 ');

4, 输入数据
体系质量 m=17.5t
刚度 k=875.5kN/m
阻尼 c=35kNs/m
荷载形式

5, 计算结果

dt=0.02s 计算结果

dt=0.3s 计算结果

6, 稳定性
结构自振周期

由于

因此选t=0.02s满足稳定性条件。