用matlab编程实现法计算多自由度综合体系的动力响应.doc

用matlab编程实现法计算多自由度体系动力响应

姓名：
学号：
班级：
专业：

用matlab编程实现法
计算多自由度体系动力响应
一、法基础原理
Newmark-b法是一个逐步积分方法，避免了任何叠加应用，能很好适应非线性反应分析。
Newmark-b法假定：
(1-1)
(1-2)
式中，b和g是按积分精度和稳定性要求进行调整参数。当b=0.5，g=0.25时，为常平均加速度法，即假定从t到t+Dt时刻速度不变，取为常数。研究表明，当b≥0.5， g≥0.25(0.5+b)2时，Newmark-b法是一个无条件稳定格式。
由式(2-141)和式(2-142)可得到用及，，表示，表示式，即有
(1-3)
(1-4)
考虑t+Dt时刻振动微分方程为：
(1-5)
将式(2-143)、式(2-144) 代入(2-145)，得到相关ut+Dt方程
(1-6)
式中

求解式(2-146)可得，然后由式(2-143)和式(2-144)可解出和。
由此，Newmark-b法计算步骤以下：
1.初始计算：
（1）形成刚度矩阵[K]、质量矩阵[M]和阻尼矩阵[C]；
（2）给定初始值， 和；
（3）选择积分步长Dt、参数b、g，并计算积分常数
，，，，
，，，；
（4）形成有效刚度矩阵；
2.对每个时间步计算：
（1）计算t+Dt时刻有效荷载：

（2）求解t+Dt时刻位移：

（3）计算t+Dt时刻速度和加速度：

Newmark-b方法是一个无条件稳定隐式积分格式，时间步长Dt大小不影响解稳定性，Dt选择关键依据解精度确定。

二、 本文使用方法计算基础问题
四层框架结构在顶部受一个简谐荷载作用，力作用时间=5s，计算响应时间为100s，分步完成。阻尼矩阵由Rayleigh阻尼结构。
具体数据以下图：

图一：结构基础计算简图
三、 计算法源程序

m=[1,2,3,4];
m=diag(m);
k= [800 -800 0 0;
-800 2400 -1600 0;
0 -1600 4800 -3200;
0 0 -3200 8000];
c=0.05*m+0.02*k;
f0=100;
t1=5;
nt=;
dt=0.01;

alfa=0.25;
beta=0.5;
a0=1/alfa/dt/dt;
a1=beta/alfa/dt;
a2=1/alfa/dt;
a3=1/2/alfa-1;
a4=beta/alfa-1;
a5=dt/2*(beta/alfa-2);
a6=dt*(1-beta);
a7=dt*beta;

d=zeros(4,nt);
v=zeros(4,nt);
a=zeros(4,nt);

for i=2:nt
t=(i-1)*dt;
if (t<t1), f=[f0*sin(4*pi*t/t1);0;0;0]; else f=[0;0;0;0]; end
ke=k+a0*m+a1*c;
fe=f+m*(a0*d(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+c*(a1*d(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));

d(:,i)=inv(ke)*fe;
a(:,i)=a0*(d(:,i)-d(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);
v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);
end

四、 计算结果截图
最终程序分别计算出四个质点位移、速度、加速度响应。
现将部分截图以下：

1、位移响应：

图二：1质点位移响应

图三：4质点位移响应

2、速度响应

图四：1质点速度响应

图五：4质点速度响应

3、加速度响应

图六：1质点加速度响应

图七：4质点加速度响应