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

用matlab编程实现Newmark - B法计算 多自由度体系的动力响应 用matlab编程实现Newmark - B法 计算多自由度体系的动力响应 、Newmark - B法的基本原理 能很好的适应非线性的Newmark-法是一种逐步积分的方法，避免了任何叠加的应用, 能很好的适应非线性的 反应分析。 Newmark-法假定： {u}t t {u}t [(1 ){u}t {u}t t] t(1-1) 1 2 {u}t t {u}t {u} t [(2 ){u}t {u}t t] t (1-2) 式中，和是按积分的精度和稳定性要求进行调整的参数。当 =，=时，为常平均加速 一 1 度法，即假定从t到t+t时刻的速度不变，取为常数-({u}t {u}t t)。研究表明，当育 +)2时，Newmark-法是一种无条件稳定的格式。 由式(2-141)和式(2-142)可得到用{u}t t及{u}t，{u}t，{u}t表示的{u}t t，{u}t t表达式， 即有 1 1 1 {u}t t 〒(蚀 t {u}t)-t{u}t(丁 1){u}t(1-3) {u}t t ({u}t t {u}t) (1 -){u}t (1 ) t{u}t(1-4) t 2 考虑t+ t时刻的振动微分方程为： [M]{u}t t [C]{u}t t [K]{u}t t {R}t t(1-5) 将式(2-143)、式(2-144)代入(2-145)，得到关于ut+t的方程 [K]{u}t t {R}t t(1-6) 式中 求解式(2-146河得{u}t t，然后由式(2-143)和式(2-144)可解出{u}t t和{u}t t 由此，Newmark-法的计算步骤如下： 10i 1 0 i 初始计算： （1）形成刚度矩阵［KI、质量矩阵［M］和阻尼矩阵［q； (3)(4)给定初始值{U}o ,选择积分步长t、形成有效刚度矩阵{U}o (3) (4) 给定初始值{U}o , 选择积分步长t、 形成有效刚度矩阵 {U}o 和{U}o ； 参数、 ，并计算积分常数 [K] [K] 2)， o[M] t(1 ) ， 7 t ； 1【C］； 对每个时间步的计算： （1） 计算t+ t时刻的有效荷载： （2） 求解t+ t时刻的位移： （3） 计算t+ t时刻的速度和加速度： Newmark-方法是一种无条件稳定的隐式积分格式，时间步长t的大小不影响解的稳定 性，t的选择主要根据解的精度确定。 本文用Newmark - B法计算的基本问题 四层框架结构在顶部受一个简谐荷载 F二F°sin（n）的作用，力的作用时间t1 =5s， t1 计算响应的时间为100s，分2000步完成。阻尼矩阵由Rayleigh阻尼构造。 具体数据如下图： 图一：结构基本计算简图 三、 计算Newmark - B法的源程序 clc clear m=［1,2,3,4］; m=diag（m）;%计算质量矩阵 k=［800-80000; -8002400-16000; 0-; 00-］;%计算刚度矩阵 c=2*m**sqrt（k/m）; t1=5;%力的作用时间 nt=2000;% 分 2000 步完成 dt=;%时间步长 alfa=;% = beta=;% = aO=1/alfa/dt/dt;% a1=beta/alfa/dt;% a2=1/alfa/dt;% 2 1 a3=1/2/alfa-1;% 3 1 2 a4=beta/alfa-1;% 4 一 1 a5=dt/2*(beta/alfa-2);%5 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 ,n t);% 初加速度为 0 fori=2:nt t=(i-1)*dt; if(tt1), f=[300*si n(6*pi*t)-50*cos(3*pi*t);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)); %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)+