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

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

1, 程序说明

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

.11..11222i i i i i u u u t

u u u u t +-+--=∆-+=∆

体系的运动方程为

...

()()()()mu t cu t ku t P t ++=

联立以上三式, 得 1111222i i i i i i i u u u u u m c ku P t t

+-+--+-++=∆∆ 上式中, 假设i u 和1i u -是已知的, 即i t 和i t 以前时刻的运动已知, 则能够把已知

项移到方程的右边, 整理得

11222222i i i i m c m m c u P k u u t t t t

t +-⎛⎫⎛⎫⎛⎫+=---- ⎪ ⎪ ⎪∆∆∆∆∆⎝⎭⎝⎭⎝⎭ 这样, 就能够计算体系任意时刻的位移, 速度和加速度。

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