newmark法程序

1用matlab编程Newmark-β法一、Newmark-β法原理Newmark-法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。Newmark-法假定：tuuuutttttt]}{}){1[(}{}{(1-1)2]}{}){21[(}{}{}{tuutuuutttttt(1-2)式中，和是按积分的精度和稳定性要求进行调整的参数。当=0.5，=0.25时，为常平均加速度法，即假定从t到t+t时刻的速度不变，取为常数)}{}({21tttuu。研究表明，当≥0.5，≥0.25(0.5+)2时，Newmark-法是一种无条件稳定的格式。由式(2-141)和式(2-142)可得到用ttu}{及tu}{，tu}{，tu}{表示的ttu}{，ttu}{表达式，即有tttttttuutuutu}){121(}{1)}{}({1}{2(1-3)tttttttutuuutu}{)21(}){1()}{}({}{(1-4)考虑t+t时刻的振动微分方程为：ttttttttRuKuCuM}{}]{[}]{[}]{[(1-5)将式(2-143)、式(2-144)代入(2-145)，得到关于ut+t的方程ttttRuK}{}]{[(1-6)式中][][1][][2CtMtKK2)}{)12(}){1(}{]([)}){121(}{1}{1]([}{}{2ttttttttutuutCuututMRR求解式(2-146)可得ttu}{，然后由式(2-143)和式(2-144)可解出ttu}{和ttu}{。由此，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）形成有效刚度矩阵][][][][10CMKK；2.对每个时间步的计算：（1）计算t+t时刻的有效荷载：)}{}{}{]([)}{}{}{]([}{}{541320ttttttttttuuuCuuuMFF（2）求解t+t时刻的位移：ttttFuK}{}{（3）计算t+t时刻的速度和加速度：tttttttuuuuu}{}{)}{}({}{320ttttttuuuu}{}{}{}{76Newmark-方法是一种无条件稳定的隐式积分格式，时间步长t的大小不影响解的稳定性，t的选择主要根据解的精度确定。二、Newmark-β法计算四层框架结构在顶部受一个简谐荷载014=sin()tFFt的作用，力的作用时间1t=5s，计算响应的时间为100s，分2000步完成。阻尼矩阵由Rayleigh阻尼构造。3具体数据如下图：图一：计算简图三、Newmark-β法程序m=[1,2,3,4];m=diag(m);k=[800-80000;-8002400-16000;0-16004800-3200;00-32008000];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;4a3=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);fori=2:ntt=(i-1)*dt;if(tt1),f=[f0*sin(4*pi*t/t1);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));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四、计算结果截图最后程序分别计算出四个质点的位移、速度、加速度响应。现将部分截图如下：51、位移响应：图二：1质点的位移图三：4质点的位移62、速度响应图四：1质点的速度图五：4质点的速度73、加速度响应图六：1质点的加速度响应图七：4质点的加速度响应