matlab单自由度的时程分析程序.doc

clear;clc;% 结构模型初始参数-m=3e3; %质量(单位：kg)k=1e6; %刚度（(单位：N/m)）kesai=0.05; %阻尼比取0.05c=2*kesai*sqrt(k*m); %阻尼系数% 读取地震波数据-acc=textread(D:处理后的smc文件51WCW_90_chnua370295.smc_090501.a,%f,headerlines,56);PGA_Max=max(abs(acc) %最大地面加速度绝对值% Newmark-beta法的基本参数-beta=1/6;gama=0.5; %按线性加速度法计算更接近真实结果，故取此组参数dt=0.02; %地震加速度时程波记录时间间隔b1=1/(beta*dt2); b2=1/(beta*dt); b3=1-1/(2*beta); %计算参数b4=gama/(beta*dt); b5=gama/beta-1; b6=(1-gama/(2*beta)*dt;ke=k+m*b1+c*b4; %等效刚度% 设定结构初始状态为零,生成向量空间存储计算值-u=zeros(100/dt,1); v=zeros(100/dt,1); a=zeros(100/dt,1);% Newmark-beta法的主计算程序-for n=2:100/dt fe=-m*acc(n)+b1*u(n-1)+b2*v(n-1)-b3*a(n-1)*m+b4*u(n-1)+b5*v(n-1)-b6*a(n-1)*c;%等效荷载 u(n)=fe/ke; a(n)=b1*u(n)-u(n-1)-b2*v(n-1)+b3*a(n-1); v(n)=b4*u(n)-u(n-1)-b5*v(n-1)+b6*a(n-1);end% 绘制结构在地震作用下的位移、速度、加速度时程曲线-subplot(3,1,1)t=(0:length(a)-1)*dt;plot(t,a) %加速度时程曲线Acc_Max=max(abs(a)title(Earthquake Response Curve Of Station 51WCW-90,fontsize,15)ylabel(Acc(cm/s2),fontsize,12)subplot(3,1,2)plot(t,v) %速度时程曲线Vel_Max=max(abs(v)ylabel(Vel(cm/s),fontsize,12)subplot(3,1,3)plot(t,u) %位移时程曲线Dis_Max=max(abs(u)xlabel(Time/s,fontsize,12)ylabel(Dis/cm,fontsize,12)% End-程序结束-小弟初次发贴，恳请达人们帮分析一下，不胜感激！其中的循环部分是根据结构动力学书上的写的，感觉问题就出在那部分了，请高人们指点一下线性加速度法是直接数值积分法求解地震反应的方法之一，本文所采用的线性加速度法参考大崎顺彦的地震动的谱分析入门第二版。具体计算公式详见大崎顺彦的地震动的谱分析入门第二版P116-P118。clear% *读入地震记录*fid = fopen(ei.txt);Accelerate,count = fscanf(fid,%g); %count 读入的记录的量time=0:0.02:(count-1)*0.02;% *线性加速度法计算各反应*%初始化各储存向量Displace=zeros(1,count); %相对位移Velocity=zeros(1,count); %相对速度AbsAcce=zeros(1,count);%绝对加速度Damp=0.05; %结构阻尼比取为0.05Tc=0.0:0.05:10; %结构自振周期Dt=0.02; %地震记录的步长%记录计算得到的反应，MDis为最大相对位移，MVel为最大相对速度%MAcc某阻尼时最大绝对加速度，用于画图MDis=zeros(1,length(Tc);MVel=zeros(1,length(Tc);MAcc=zeros(1,length(Tc);t=1; %在下一个循环中控制不同的结构自振周期for T=0.0:0.05:10Frcy=2*pi/T ; %结构自振频率DamFrcy=Frcy*sqrt(1-Damp*Damp);%计算公式化简e_t=exp(-Damp*Frcy*Dt);s=sin(DamFrcy*Dt);c=cos(DamFrcy*Dt);A=zeros(2,2);