结构动力学newmark法程序

用matlab编程实现Newmark-β法计算多自由度体系的动力响应
姓名：***
学号：**************
专业：结构工程
用matlab 编程实现Newmark -β法 计算多自由度体系的动力响应
一、Newmark -β法的基本原理
Newmark-β法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。

Newmark-β法假定：
t u u u u
t t t t t t ∆ββ∆∆]}{}){1[(}{}{+++-+= (1-1)
2]}{}){2
1
[(}{}{}{t u u t u
u u t t t t t t ∆γγ∆∆∆+++-++= (1-2) 式中，β和γ是按积分的精度和稳定性要求进行调整的参数。

当β=0.5，γ=0.25时，为常平均加速度法，即假定从t 到t +∆t 时刻的速度不变，取为常数)}{}({2
1
t t t u u ∆++ 。

研究表明，当β≥0.5， γ≥0.25(0.5+β)2时，Newmark-β法是一种无条件稳定的格式。

由式(2-141)和式(2-142)可得到用t t u ∆+}{及t u }{，t u
}{ ，t u }{ 表示的t t u ∆+}{ ，t t u ∆+}{ 表达式，即有
t t
t t t t t u u t u u t u
}){121
(}{1)}{}({1}{2
----=++γ∆γ∆γ∆∆ (1-3) t t t t t t t u t u
u u t u
}{)21(}){1()}{}({}{ ∆γ
β
γβ∆γβ∆∆-+-+-=++ (1-4) 考虑t +∆t 时刻的振动微分方程为：
t t t t t t t t R u K u C u
M ∆∆∆∆++++=++}{}]{[}]{[}]{[ (1-5) 将式(2-143)、式(2-144) 代入(2-145)，得到关于u t +∆t 的方程
t t t t R u K ∆∆++=}{}]{[ (1-6)
式中
][][1
][][2
C t M t
K K ∆γβ∆γ++
= )}{)12(}){1(}{]([)}){121
(}{1}{1](
[}{}{2
t t t t t t t t u t u
u t C u u t u t
M R R ∆γ
β
γβ∆γβγ∆γ∆γ∆-+-++-+++=+
求解式(2-146)可得t t u ∆+}{，然后由式(2-143)和式(2-144)可解出t t u
∆+}{ 和t t u ∆+}{ 。

由此，Newmark-β法的计算步骤如下：
1.初始计算：
（1）形成刚度矩阵[K ]、质量矩阵[M ]和阻尼矩阵[C ]；
（2）给定初始值0}{u ， 0}{u
和0}{u ； （3）选择积分步长∆t 、参数β、γ，并计算积分常数
2
01t ∆γα=
，t ∆γβα=1，t ∆γα12
=，1213-=γα， 14-=
γβα，)2(25-=γ
β
∆αt ，)1(6β∆α-=t ，t ∆βα=7； （4）形成有效刚度矩阵][][][][10C M K K αα++=； 2.对每个时间步的计算：
（1）计算t +∆t 时刻的有效荷载：
)}{}{}{]([)}{}{}{]([}{}{541320t t t t t t t t t t u u
u C u u
u M F F αααααα∆∆++++++=++
（2）求解t +∆t 时刻的位移：
[]t t t
t F u K ∆+∆+=}{}
{
（3）计算t +∆t 时刻的速度和加速度：
t t t t t t t u u u u u
}{}{)}{}({}{320 ααα∆∆---=++ t t t t t t u u u u
∆∆αα++++=}{}{}{}{76 Newmark-β方法是一种无条件稳定的隐式积分格式，时间步长∆t 的大小不影响解
的稳定性，∆t 的选择主要根据解的精度确定。

二、 本文用Newmark -β法计算的基本问题
四层框架结构在顶部受一个简谐荷载01
4=sin(
)t
F F t π的作用，力的作用时间1t =5s ，计算响应的时间为100s ，分2000步完成。

阻尼矩阵由Rayleigh 阻尼构造。

具体数据如下图：
图一：结构基本计算简图三、计算Newmark-β法的源程序
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=2000;
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质点的位移响应图三：4质点的位移响应
图四：1质点的速度响应图五：4质点的速度响应
图六：1质点的加速度响应图七：4质点的加速度响应。