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 -β法的源程序

clc clear

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=2*m*0.05*sqrt(k/m) ;

t1=5; %力的作用时间 nt=2000; %分2000步完成 dt=0.01; %时间步长 alfa=0.25; %γ=0.25 beta=0.5; %β=0.5

a0=1/alfa/dt/dt; % 2

01t

∆γα=

a1=beta/alfa/dt; %t

∆γβ

α=

1 a2=1/alfa/dt; %t

∆γα1

2=

a3=1/2/alfa-1; %121

3-=

γ

α a4=beta/alfa-1; %14-=

γ

βα a5=dt/2*(beta/alfa-2); % )2(25-=

γ

β

∆αt a6=dt*(1-beta); %)1(6β∆α-=t a7=dt*beta; % t ∆βα=7

d=zeros(4,nt); %初位移为0 v=zeros(4,nt); % 初速度为0 a=zeros(4,nt); % 初加速度为0 for i=2:nt