newmark法程序法计算多自由度体系的动力响应

用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{u2式中，和是按积分的精度和稳定性要求进行调整的参数。当=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{u111}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的方程

[K]{u}tt{R}tt (1-6)

式中

}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1111， ，，，1322t2tt4t1，5(2)，6t(1)，7t； 2（4）形成有效刚度矩阵[K][K]0[M]1[C]； 2.对每个时间步的计算： （1）计算t+t时刻的有效荷载： （2）求解t+t时刻的位移： （3）计算t+t时刻的速度和加速度： Newmark-方法是一种无条件稳定的隐式积分格式，时间步长t的大小不影响解的稳定性，t的选择主要根据解的精度确定。 二、 本文用Newmark-β法计算的基本问题 四层框架结构在顶部受一个简谐荷载F=F0sin(4t)的作用，力的作用时间t1t1=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; % 1

02ta1=beta/alfa/dt; %1 t1

a2=1/alfa/dt; %2ta3=1/2/alfa-1; %311 2a4=beta/alfa-1; %41 t(2) 2a5=dt/2*(beta/alfa-2); % 5a6=dt*(1-beta); %6t(1a7=dt*beta; % ) 7t d=zeros(4,nt); %初位移为0 v=zeros(4,nt); % 初速度为0 a=zeros(4,nt); % 初加速度为0 for i=2:nt t=(i-1)*dt; if (t % t+t

时刻的有效荷载

d(:,i)=inv(ke)*fe; %求解t+t时刻的位移 a(:,i)=a0*(d(:,i)-d(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1); %计算t+t时刻的加速度

v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i); %计算t+t时刻的速度

end

T=[0:dt:19.99]; %离散系统dt为采样周期 19.99为终端时间 close all

figure %控制窗口数量 plot(T,[d]) %绘制位移函数图像

title('?÷?êμ???ò?×üí?') %添加标题为各质点位移总图 legend('?êμ?ò?','?êμ??t','?êμ?èy','?êμ???') %添加图例的标注 xlabel('ê±??(s)') %对x轴进行标注为时间（s） ylabel('??ò?(m)') %对y轴进行标注为位移（m） grid %显示画图中的个网线 figure

plot(T,[a]) %绘制加速度函数图像 title('?÷?êμ??ó?ù?è×üí?') %添加标题为各质点加速度总图 legend('?êμ?ò?','?êμ??t','?êμ?èy','?êμ???')

xlabel('ê±??(s)') %对x轴进行标注为时间（s） ylabel('?ó?ù?è(m/s^2)') %对y轴进行标注为加速度（m/s） grid % figure plot(T,[v]) %绘制速度函数图像 title('?÷?êμ??ù?è×üí?') %添加标题为各质点速度总图 legend('?êμ?ò?','?êμ??t','?êμ?èy','?êμ???') xlabel('ê±??(s)') %对x轴进行标注为时间（s） ylabel('?ù?è(m/s^2)') %对y轴进行标注为加速度（m/s） grid 22

四、 计算结果截图 最后程序分别计算出四个质点的位移、速度、加速度响应。 现将部分截图如下： 1、位移响应： 图二：1质点的位移响应 图三：4质点的位移响应 2、速度响应 图四：1质点的速度响应 图五：4质点的速度响应 3、加速度响应

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