地震工程作业哈工大MATLAB程序

1、地震工程大作业 哈工大 李金平 弹性反应谱 原创性声明，主程序由本人独立编写完成，支持各种检验。所选地震动RSN2345CHICHI.AT2 peer索取号：2345:位移反应谱局部放大求比值绝对加速度反应谱局部放大求比值进行比较速度反应谱局部放大观察差异求比值进行比较。结论：对比：拟合效果非常好，短周期 0sT0.2s误差较大，其中相对速度谱的误差较大，接近30%，加速度和位移误差都在3%以内，周期0.2sT10s各反应谱的误差都保持在1%以内，精度相同，猜想，两种计算程序应该用的同一种计算方法（newmark法）。理解：本次计算的反应谱是在给定的地震动下，具有相同阻尼比（5%）不同周期的单

2、自由度结构的线弹性反应幅值，得到的速度、加速度、移幅值随周期变化的三条曲线。a1(i)=max(abs(a+xg); v1(i)=max(abs(v);x1(i)=max(abs(x);用法：今后给定一个结构，我们可以计算其周期T1，然后在相应的反应谱图表里面找到对应的谱值如Sa 、Sv、 Sd，即为我们所要求得的该结构在指定地震动作用下最大动力反应，将动力分析简化成为静力计算，简单方面。（弹性范围内）主程序 DZZY.mclearclcM=1;%fid=fopen(RSN2345CHICHI.txt);%读取地震动加速度记录xg=9.8*fscanf(fid,%f);fclose(fid);

3、n=length(xg);F=-M*xg; %生成地震力for i=1:1000 tn(i)=0.01*i; K=2*pi*2*pi/tn(i)/tn(i); lamda=2*pi/tn(i);dt=0.005;t=(0:0.005:(n-1)*0.005);x0=0;v0=0;a0=0;ksi=0.01*5;DAMPER=2*ksi*lamda*M;x,v,a=newmarkb(M,K,DAMPER,1,F,x0,v0,a0,dt,n);a1(i)=max(abs(a+xg);v1(i)=max(abs(v);x1(i)=max(abs(x);endtest=importdata(a2.tx

4、t); %读取Seismosignal绝对加速度a2=test(:,2)*9.8;testv=importdata(v2.txt); %读取Seismosignal速度v2=testv(:,2);testd=importdata(x2.txt); %读取Seismosignal相对位移x2=testd(:,2);tn=0.01:0.01:10;figure (1) plot(tn,x1*1000,r-)xlabel(周期T (s)ylabel(位移 (mm)hold onplot(tn,x2*1000,linewidth,2)legend(matlab相对位移,Seismosignal相对位移

5、)% % figure (2)plot(tn,a1,r-,linewidth,1)xlabel(周期T (s)ylabel(Sa (m/s2)hold onplot(tn,a2,k,linewidth,2)legend(matlab-Sa,Seismosignal-Sa)figure (3)plot(tn,v1)xlabel(周期T (s)ylabel(速度 (m/s)hold onplot(tn,v2,linewidth,2)legend(matlab-Sv,Seismosignal-Sv)hold offfigure (4)hold onplot(t,xg)legend(地震动加速度)xl

6、abel(时间 (s)ylabel(加速度 (m/s2)for j=1:1000 %求比值ca(j)=a2(j)/a1(j);cv(j)=v2(j)/v1(j);cx(j)=x2(j)/x1(j);end figure (6)plot(tn,ca)legend(Seismosignal加速度/MATLAB加速度比值)xlabel(周期T (s) ylabel(比值)figure (7)plot(tn,cv)legend(Seismosignal速度/MATLAB速度 比值)xlabel(周期T (s) ylabel(比值)figure (8)plot(tn,cx)legend(Seismosi

7、gnal位移/MATLAB位移 比值)xlabel(周期T (s) ylabel(比值)子程序 newmarkb.mfunction x,v,a=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)% newmark-beta method% obtain the response of the dynamic system% x,v,a=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)% M - mass matrix% K - stiffness matrix% C - damping matrix% N - DO

8、F% P - loads% x0 - initial displacement% v0 - initial velocity% a0 - initial acceleration% dt - interval% RecordLength - number of sampling pointsx=zeros(N,RecordLength);v=zeros(N,RecordLength);a=zeros(N,RecordLength); x(:,1)=x0;v(:,1)=v0;a(:,1)=a0;deta=0.50;alpha=0.25;a0=1/alpha/dt2;a1=deta/alpha/dt;a2=1/alpha/dt;a3=1/2/alpha-1;a4=deta/alpha-1;a5=dt*(deta/alpha-2)/2;a6=dt*(1-deta);a7=deta*dt;K_=K+a0*M+a1*C;iK=inv(K_); for i=1:RecordLength-1 P_(:,i+1)=P(:,i+1)+M*(a0*x(:,i)+a2*v(:,i)+a3*a(:,i)+C*(a1*x(:,i