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

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

姓名：

学号： 班级： 专业：

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

一、Newmark-β法的基本原理

Newmark- 法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。 Newmark- 法假定：

}t t {u }t [(1 ){u }t {u }t t] t {u

(1-1)

1

} t [( ){u }t {u }t t] t2 (1-2) {u}t t {u}t {u

2

式中， 和 是按积分的精度和稳定性要求进行调整的参数。当 =0.5， =0.25时，为常平均加速度法，即假定从t到t+ t时刻的速度不变，取为常数1

}t {u }t t)。研究表明，当 ≥0.5， ≥0.25(0.5+ )2时，Newmark- 法是一种({u2

无条件稳定的格式。

}t t，{u }t，{u }t表示的{u }t t由式(2-141)和式(2-142)可得到用{u}t t及{u}t，{u表达式，即有

}t t {u

111

}t (1-3) ({u} {u}) {u} ( 1){ut ttt2

t2 t

}t t {u

}t (1 ) t{u }t (1-4) ({u}t t {u}t) (1 u

t 2

考虑t+ t时刻的振动微分方程为：

}t t [C]{u }t t [K]{u}t t {R}t t (1-5) [M]{u

将式(2-143)、式(2-144) 代入(2-145)，得到关于ut+ t的方程

[]{u}t t {t t (1-6)

式中

[] [K]

1

[M] [C] 2

t t

111

}t){u} {u} ( 1){utt2

t2 t

{ {R}t t [M

}t ( 1) t{u }t) [C{u}t ( 1){u

t 2

}t t和{u }t t。 求解式(2-146)可得{u}t t，然后由式(2-143)和式(2-144)可解出{u由此，Newmark- 法的计算步骤如下：

1.初始计算：

（1）形成刚度矩阵[K]、质量矩阵[M]和阻尼矩阵[C]；

}0和{u }0； （2）给定初始值{u}0， {u

（3）选择积分步长 t、参数 、 ，并计算积分常数

0

111 ，，， 1， 231

t2 t t2

4

t

1， 5 ( 2)， 6 t(1 )， 7 t； 2

（4）形成有效刚度矩阵[] [K] 0[M] 1[C]； 2.对每个时间步的计算：

（1）计算t+ t时刻的有效荷载：

}t 3{u }t){t t {F}t t [M]( 0{u}t 2{u

}t 5{u }t) [C]( 1{u}t 4{u

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

{u}

t t

{t t

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

}t t 0({u}t t {u}t) 2{u }t 3{u }t {u

}t t {u }t 6{u }t 7{u }t t {u

Newmark- 方法是一种无条件稳定的隐式积分格式，时间步长 t的大小不影响解

的稳定性， t的选择主要根据解的精度确定。

二、 本文用Newmark-β法计算的基本问题

四层框架结构在顶部受一个简谐荷载F=F0sin(

4 t

)的作用，力的作用时间t1

t1=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、位移响应：

图二：1质点的位移响应

图三：4质点的位移响应

2、速度响应

图四：1质点的速度响应

图五：4质点的速度响应

3、加速度响应

图六：1质点的加速度响应

图七：4质点的加速度响应