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

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

当β=0.5，γ=0.25时，为常平均加速度法，即假定从t 到t +∆t 时刻的速度不变，取为常数)}{}({21t 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 tu}){121(}{1)}{}({1}{2----=++γ∆γ∆γ∆∆ (1-3) t t t t t t t u t uu u t u}{)21(}){1()}{}({}{ ∆γβγβ∆γβ∆∆-+-+-=++ (1-4) 考虑t +∆t 时刻的振动微分方程为：t t t t t t t t R u K u C uM ∆∆∆∆++++=++}{}]{[}]{[}]{[ (1-5) 将式(2-143)、式(2-144) 代入(2-145)，得到关于u t +∆t 的方程t t t t R u K ∆∆++=}{}]{[ (1-6)式中][][1][][2C t M tK K ∆γβ∆γ++= )}{)12(}){1(}{]([)}){121(}{1}{1]([}{}{2t t t t t t t t u t uu t C u u t u tM 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 、参数β、γ，并计算积分常数201t ∆γα=，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 uu C u uu M F F αααααα∆∆++++++=++（2）求解t +∆t 时刻的位移：[]t t tt 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 -β法计算的基本问题四层框架结构在顶部受一个简谐荷载014=sin()tF F t π的作用，力的作用时间1t =5s ，计算响应的时间为100s ，分2000步完成。

阻尼矩阵由Rayleigh 阻尼构造。

具体数据如下图：图一：结构基本计算简图三、 计算Newmark -β法的源程序clc clearm=[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.5a0=1/alfa/dt/dt; % 201t∆γα=a1=beta/alfa/dt; %t∆γβα=1 a2=1/alfa/dt; %t∆γα12=a3=1/2/alfa-1; %1213-=γα 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 ∆βα=7d=zeros(4,nt); %初位移为0 v=zeros(4,nt); % 初速度为0 a=zeros(4,nt); % 初加速度为0 for i=2:ntt=(i-1)*dt;if (t<t1),f=[300*sin(6*pi*t)-50*cos(3*pi*t);0;0;0]; %力作用时间内对结构进行加载elsef=[0;0;0;0]; %力作用时间外结构不受力endke=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));% 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时刻的速度endT=[0:dt:19.99]; %离散系统dt为采样周期19.99为终端时间close allfigure %控制窗口数量plot(T,[d]) %绘制位移函数图像title('¸÷ÖʵãλÒÆ×Üͼ') %添加标题为各质点位移总图legend('ÖʵãÒ»','Öʵã¶þ','ÖʵãÈý','ÖʵãËÄ') %添加图例的标注xlabel('ʱ¼ä(s)') %对x轴进行标注为时间（s）ylabel('λÒÆ(m)') %对y轴进行标注为位移（m）grid %显示画图中的个网线figureplot(T,[a]) %绘制加速度函数图像title('¸÷Öʵã¼ÓËÙ¶È×Üͼ') %添加标题为各质点加速度总图legend('ÖʵãÒ»','Öʵã¶þ','ÖʵãÈý','ÖʵãËÄ')xlabel('ʱ¼ä(s)') %对x轴进行标注为时间（s）ylabel('¼ÓËÙ¶È(m/s^2)') %对y轴进行标注为加速度（m/s2）grid%figureplot(T,[v]) %绘制速度函数图像title('¸÷ÖʵãËÙ¶È×Üͼ') %添加标题为各质点速度总图legend('ÖʵãÒ»','Öʵã¶þ','ÖʵãÈý','ÖʵãËÄ')xlabel('ʱ¼ä(s)') %对x轴进行标注为时间（s）ylabel('ËÙ¶È(m/s^2)') %对y轴进行标注为加速度（m/s2）grid四、计算结果截图最后程序分别计算出四个质点的位移、速度、加速度响应。

现将部分截图如下：1、位移响应：图二：1质点的位移响应图三：4质点的位移响应2、速度响应图四：1质点的速度响应图五：4质点的速度响应3、加速度响应图六：1质点的加速度响应图七：4质点的加速度响应。