地震工程作业

地震工程课程大作业作业1绘制1940ElCentro，N-S分量地震动的绝对加速度、相对速度和相对位移反应谱。地震动：在PEERGroundMotionDatabase自行下载经典的1940ElCentro，N-S分量。要求：在此模板内完成，A4纸打印。自编程序与软件（Bispec或Seismosigna等）计算反应谱进行对比。提交自编写程序。Matlab程序：clearfid=fopen('E:\Earthquake\Elcentro.txt');[Accelerate,count]=fscanf(fid,'%g');%count读入的记录的量Accelerate=9.8*Accelerate';%单位统一为m和stime=0:0.02:(count-1)*0.02;%单位sDisplace=zeros(1,count);%相对位移Velocity=zeros(1,count);%相对速度AbsAcce=zeros(1,count);%绝对加速度DampA=[0.00,0.02,0.05];%三个阻尼比TA=0.0:0.02:4;%TA=0.000001:0.02:4;结构周期Dt=0.02;%地震记录的步长%记录计算得到的反应，MDis为某阻尼时最大相对位移%MVel为某阻尼时最大相对速度，MAcc某阻尼时最大绝对加速度MDis=zeros(3,length(TA));MVel=zeros(3,length(TA));MAcc=zeros(3,length(TA));j=1;forDamp=[0.00,0.05,0.1]t=1;forT=0.0:0.02:4Frcy=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);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DamFrcy;A(2,1)=-Frcy*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(Frcy^2*Dt);d_3t=Damp/(Frcy^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy^2)*c)-2*d_3t;B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/Frcy)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/Frcy^2)*(DamFrcy*s+Damp*Frcy*c))+1/(Frcy^2*Dt);B(2,2)=e_t*(1/(Frcy^2*Dt)*c+s*Damp/(Frcy*DamFrcy*Dt))-1/(Frcy^2*Dt);fori=1:(count-1)Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i+1);Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i+1);AbsAcce(i+1)=-2*Damp*Frcy*Velocity(i+1)-Frcy^2*Displace(i+1);endMDis(j,t)=max(abs(Displace));MVel(j,t)=max(abs(Velocity));ifT==0.0MAcc(j,t)=max(abs(Accelerate));elseMAcc(j,t)=max(abs(AbsAcce));endDisplace=zeros(1,count);Velocity=zeros(1,count);AbsAcce=zeros(1,count);地震工程课程大作业t=t+1;endj=j+1;endcloseallfigure%绘制位移反应谱plot(TA,MDis(1,:),'-b',TA,MDis(2,:),'-r',TA,MDis(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=0','ζ=0.02','ζ=0.05')gridfigure%绘制速度反应谱plot(TA,MVel(1,:),'-b',TA,MVel(2,:),'-r',TA,MVel(3,:),':k')title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=0','ζ=0.02','ζ=0.05')gridfigure%绘制绝对加速度反应谱plot(TA,MAcc(1,:),'-b',TA,MAcc(2,:),'-r',TA,MAcc(3,:),':k')title('AbsoluteAcceleration')xlabel('Tn(s)')ylabel('absoluteacceleration(m/s^2)')legend('ζ=0','ζ=0.02','ζ=0.05')grid由此得到下图反应谱与由Seismosigna得到的反应谱作出比较相对位移反应谱（上）与Seismosigna计算相对位移反应谱（下）的比较地震工程课程大作业00.511.522.533.5400.10.20.30.40.50.60.7DisplacementTn(s)Displacement(m)ζ=0ζ=0.02ζ=0.05相对速度反应谱（上）与Seismosigna计算相对速度反应谱（下）的比较00.511.522.533.5400.511.522.5VelocityTn(s)velocity(m/s)ζ=0ζ=0.02ζ=0.05地震工程课程大作业绝对加速度反应谱（左图）与Seismosigna计算绝对加反应谱（右图）的比较00.511.522.533.5405101520253035AbsoluteAccelerationTn(s)absoluteacceleration(m/s2)ζ=0ζ=0.02ζ=0.05作业2、地震工程课程大作业使用中心差分方法和Newmark方法分别计算下面单自由度结构的地震反应。题：考虑下面的具有2个自由度的一个简单系统，分析时间步长Δt=0.28sec的情况，求每时刻的加速度、速度、位移求：要求：在此模板内完成，A4纸打印。MATLAB程序：clear%求自振圆频率与振型m=[20;01];k=[6-2;-24];[v,d]=eig(inv(m)*k);w=sqrt(d);fai1=v(:,1)./v(1,1);fai2=v(:,2)./v(1,2);fai=[fai1fai2];mm=fai'*m*fai;kk=fai'*k*fai;%等效质量矩阵及等效刚度矩阵p0=[0;10];%中心差分法s0=[0;0];u0=[0;0];%初速度及初位移均为0%以下过程均采用q、q1、q2的形式分别表示位移、速度、加速度q0=[fai1'*m*s0./(fai1'*m*fai1);fai2'*m*s0./(fai2'*m*fai2)];q01=[fai1'*m*u0./(fai1'*m*fai1);fai2'*m*u0./(fai2'*m*fai2)];P0=fai'*p0;q02=mm^(-1)*(P0-kk*q0);t=0.28;q(:,1)=q0-t*q01+t^2*q02/2;%q(:,1)表示q-1，即初始时刻前一时刻的位移q(:,2)=[0;0];%q(:,2)表示初始时刻的位移00U00U地震工程课程大作业Kh=mm/(t^2);a=mm/(t^2);b=kk-2*mm/(t^2);P=fai'*p0;fori=2:7;Ph(:,i)=P-a*q(:,i-1)-b*q(:,i);q(:,i+1)=Kh^(-1)*Ph(:,i);q1(:,i)=(q(:,i+1)-q(:,i-1))/(2*t);q2(:,i)=(q(:,i+1)-2*q(:,i)+q(:,i-1))/(t^2);endu=fai*q;u1=fai*q1;u2=fai*q2;fori=1:6U(:,i)=u(:,i+1);U1(:,i)=u1(:,i+1);U2(:,i)=u2(:,i+1);endUU1U2%Newmark法（采用线加速度法）gamma=1/2;beta=1/6;Kh1=kk+mm/(beta*t^2);a1=mm/(beta*t);b1=mm/(2*beta);detaP=[00000;00000];q1(:,1)=[0;0];q11(:,1)=[0;0];q12(:,1)=q02;fori=1:5detaPh(:,i)=detaP(:,i)+a1*q11(:,i)+b1*q12(:,i);detaq(:,i)=Kh1^(-1)*detaPh(:,i);detaq1(:,i)=gamma*detaq(:,i)/(beta*t)-gamma*q11(:,i)/beta+t*(1-gamma/(2*beta))*q12(:,i);detaq2(:,i)=detaq(:,i)/(beta*t^2)-q11(:,i)/(beta*t)-q12(:,i)/(2*beta);地震工程课程大作业q1(:,i+1)=q1(:,i)+detaq(:,i);q11(:,i+1)=q11(:,i)+detaq1(:,i);q12(:,i+1)=q12(:,i)+detaq2(:,i);endv=fai*q1;v1=fai*q11;v2=fai*q12;fori=1:6V(:,i)=v(:,i);V1(:,i)=v1(:,i);V2(:,i)=v2(:,i);endVV1V2表1中心差分法Time0.28s0.56s0.84s1.12s1.40s位移（m）00.03070.16750.48711.01700.39201.44512.83384.14415.0152速度（m/s）0.05490.29920.81491.51692.16752.58054.36034.81973.89541.9873加速度(m/s^2)0.39201.35292.33122.68281.96428.4320402812-1.0001-5.6022-8.0268表2Newmark法Time0.28s0.56s0.84s1.12s1.40s位移（m）0.00470.04440.18260.48500.97800.37261.38092.73174.04474.9744速度（m/s）0.05020.27510.75551.42382.07202.59264.42444.98504.17712.3360加速度(m/s^2)0.35861.24762.18392.58962.04038.51884.5654-0.5615-5.2088-7.9416作业3、如图1所示一3层钢筋混凝土框架结构FR3，各层层高均为4.2m。设计荷载：楼面恒荷5.0kN/m2，楼面活荷载2.0kN/m2，屋面恒荷载5.0kN/m2，屋面活荷载2.0kN/m2，结构构件尺寸和配筋见表1，该结构遭受一地震动作用。使用非线性静力分析方法或时程分析方法求取该地震动作用下此结构（图1右图FR3-A轴平面框架。不考虑楼板贡献。）的顶点地震工程课程大作业最大侧向位移值（参考使用的软件：OpenSees、Abaqus、Ansys、Sap2000、Idarc等）