单自由度地震能量计算的程序matlab.docxVIP

funct ion responsenew %基于戏线性滞冋模型的单自由度体系的地震能虽分析程序 %made by Yanan Li Faculty of Infrastructure Enginoering 弔Dalian University of Technology %质量57041kg,阻尼36612 N s/叫初始刚度2350000N/m,刚度折减系数0. 2,屈服位移0?01m, 采用ELCENTRO波 %参数替换直接在下面修改，然后运行 clc format long; m=57041;% 质量 ug= importdataC* ELCENTRO. txt) ;%地震波 txt 文件 ug=ug/100; P=-m*ug; num=size(P, 1); c 二36612; % 阻尼 k 1=2350000; %初始刚度 k2二kl*0?2; a=zeros(lt num) ;v=zeros(l, num) ;x=zeros(l, num);%加速度 速度 位移 Ei=zeros(l, num);Ek二zeros(1, num);Ed二zeros(1, num); Eh二zeros (1, num)输入能 动能 阻 尼耗能滞回耗能 EI=zeros(lt num); EK二zeros (1, num) ;ED=zeros(l, num) ;EH=zeros(l, num) ;%累积的各种能量 t imc=zeros(l, num); a(l)=P (l)/m; hfl=zeros (1, num); h=0. 02;%地震波采样间隔 Xy=0.01;%屈服位移 pxmax=0; nxmax-0; pd=0;%双线型滞冋模型折线的标识，0农示弹性，1表示正向弹犁性，2表示反向弹性，-1 表示反向弹塑性，-2表示正向弹性 for ii=l:num if pd==0 %弹性阶段 k二kl; if x(ii)Xy pd=l; b=[(a(ii)-a(ii-l))/6/h a(ii~l)/2 v(ii-l) x(iiT)-Xy];% 拐点处理 d=roots(b); for j=l:length(d) if isreal (d(j))==l dth=d(j); end end hp二h-dth; vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth*2/h/2); ap=a(i i-l) + (a (i i)-a(ii-l)) *dth/h; pp=m*ap+c*vp+kl*Xy; kd=k2+3*c/hp+6*m/hp/hp; dtpd=P(i i +1) -pp+m* (6*vp/hp+3*ap) +c* (3*vp^-hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*vp-hp*ap/2; x(ii)=Xy+dtx; v(ii)=vp^dtv; dtf=k2*dtx; hfl(ii)=kl*Xy+dtf; a(ii) = (P(ii)-hfl (i i+1)-c*v(i i))/m; elseif x(ii)-Xy pd=-l; b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(iiT) x(iiT)+Xy];% 拐点处理 d=roots(b); for j=l:length(d) if isreal (d(j))==l dth=d(j); end end hp=h-dth; vp=v(i i 1)?a(i i 1)*dth+ ((a(ii)-a(ii-l))*dth2/h/2); ap a(ii-l) + (a(ii)-a(ii-l))*dth/h; pp=m*ap+c*vp_k1*Xy; kd=k2+3*c/hp+6*m/hp/hp; dtpd=P (i i + 1)-pp+m*(6*vp/hp+3*ap)(3*vp^hp*ap/2); dtx=dtpd/kd; dtv-3*dtx/hp-3*vp-hp*ap/2; x(ii)=-Xy+dtx： v(i i)=vp^dtv; dtf=k2*dtx; hfl(ii)=-kl*Xy+dtf; a(i i) = (P(i i)-hf 1 (i i)-c*v(ii))/m; end end if pd==l %正向弹塑性 k=k2; if v(ii)0 pd=2; b=[(a(ii)-a(ii-l))/2/h a(ii-l) v(ii-l)];% 拐点处理 d二rools(b); for j=l:length(d) if isreal (d(j))==l dth=d (j); end end hp=h-dth; xp=x(iiT)+v(i i-1) *dth+a (i i-1) *dth