www.gusucode.com > matlab非线性混凝土模型计算分析源码程序 > main_a20.m
%该程序为主程序 clear;%清除工作空间变量; clc;format compact;format long g; %per=[0:0.2:1.0 1.5:0.5:2 3:1:6];%设定自振周期范围,分别采用不同的密度; per=0:0.2:6; sp=length(per); %采用不同的屈服强度系数; ebx(1)=0.167; ebx(2)=0.2; ebx(3)=0.25; ebx(4)=0.333; ebx(5)=0.5; for iy=1:5 eb=ebx(iy); %=加载地震波=------------------------------------------------------------------- %next one 192 fida(1)=fopen('e:/berkleywave/northridgenew/nor_php270_a.txt','r+'); fidv(1)=fopen('e:/berkleywave/northridgenew/nor_php270_v.txt','r+'); %next one 219 fida(2)=fopen('e:/berkleywave/northridgenew/nor_sar000_a.txt','r+'); fidv(2)=fopen('e:/berkleywave/northridgenew/nor_sar000_v.txt','r+'); %next one 220 fida(3)=fopen('e:/berkleywave/northridgenew/nor_sar270_a.txt','r+'); fidv(3)=fopen('e:/berkleywave/northridgenew/nor_sar270_v.txt','r+'); %next one 227 fida(4)=fopen('e:/berkleywave/northridgenew/nor_scs052_a.txt','r+'); fidv(4)=fopen('e:/berkleywave/northridgenew/nor_scs052_v.txt','r+'); %=加载地震波=------------------------------------------------------------------- um=zeros(1,sp); rm=zeros(1,sp); nw=4; ztd=0;zacv=0; for i=1:nw tag=0.;tvg=0.;ag=0.;vg=0.; for di=1:5 aline=fgetl(fida(i)); vline=fgetl(fidv(i)); if(di==4) dt=fscanf(fida(i),'%f',[1 inf]); end end ag=fscanf(fida(i),'%f',[1 inf])*9.81; vg=fscanf(fidv(i),'%f',[1 inf])/100; status=fclose(fida(i));status=fclose(fidv(i)); st=length(ag);% at=dt*(st-1); %总时间为at; t=0:dt:at; cmag=max(abs(ag)); %ag=ag*2.2/cmag;%%%为和肖结果比较,特添加; ag=ag*0.4*9.81/cmag;;%%%和CHAI比较,假定为0.4g,调整mvg; vg=vg*0.4*9.81/cmag; mag=max(abs(ag));mvg=max(abs(vg)); acv=mag/9.81/mvg; zacv=zacv+acv; tc=2*pi*2.0/2.5*mvg/mag;%计算场地特征周期Tc;cv=2.0;ca=2.5; %%%&***&***能量分析应采用相同震级地面运动加速度进行标准化; %%%计算强震时间td;td=t0.95-t0.05; ia=zeros(1,st); for m=2:st ia(m)=ia(m-1)+pi/2/9.81*(ag(m-1)^2+ag(m)^2)*dt/2; end iac1=0.1; iac2=0.1; for n=1:st ia1=abs(ia(n)/ia(st)-0.05); ia2=abs(ia(n)/ia(st)-0.95); if(ia1<iac1) iac1=ia1; t1=(n-1)*dt; end if(ia2<iac2) iac2=ia2; t2=(n-1)*dt; end end td=t2-t1; ztd=ztd+td; %%%做线弹性计算; para1=1; [rd1,ad1,rv1,av1,ra1,aa1]=response_l(para1,per,ag,dt,at); % figure(1) % mesh(per,t,ad1) %绘制出绝对位移时频反应谱图; % grid on % xlabel('周期 / 秒'); % ylabel('时间 / 秒'); [row,rank]=size(aa1); sa=zeros(1,rank); for ri=1:rank for rj=1:row if (abs(aa1(rj,ri))>sa(ri)) sa(ri)=aa1(rj,ri);%弹性加速度反应谱sa(最大绝对加速度反应); end end end %%%做弹塑性计算; para1=1; [u1,nrd1,nad1,nrv1,nav1,nra1,naa1,rep]=response_n(para1,sa,per,ag,dt,at,eb); % figure(2) % mesh(per,t,nad1) %绘制出弹塑性位移时频反应谱; % grid on % xlabel('周期 / 秒'); % ylabel('时间 / 秒'); um1=zeros(1,sp); for np=1:sp for nt=1:st if (abs(u1(nt,np))>um1(np)) um1(np)=abs(u1(nt,np)); end end end um=um+um1; rm=rm+rep; end um=um/nw; rm=rm/nw; ptd=ztd/nw;pacv=zacv/nw; %%%输出td和a/v数据到硬盘; save('e:\north_concrete\a20\pingjinzhi.txt','ptd','pacv','-ASCII') if(iy==1) figure(1) plot(per(2:sp),um(2:sp),'kd-.') %绘制出弹塑性位移延性反应谱; grid on xlabel('周期T(s)'); ylabel('位移延性系数u'); title('弹塑性位移延性反应谱'); hold on; figure(2) %%%输出数据到硬盘; save('e:\north_concrete\a20\em_a20_1.txt','rm','-ASCII') plot(per,rm,'kd-.') %绘制出累积滞回耗能谱; grid on xlabel('周期T(s)'); ylabel('滞回耗能Eh(J/kg)'); title('弹塑性滞回耗能谱'); hold on; elseif(iy==2) figure(1) plot(per(2:sp),um(2:sp),'kh-.') %绘制出弹塑性位移延性反应谱; hold on; figure(2) %%%输出数据到硬盘; save('e:\north_concrete\a20\em_a20_2.txt','rm','-ASCII') plot(per,rm,'kh-.') %绘制出累积滞回耗能谱; hold on; elseif(iy==3) figure(1) plot(per(2:sp),um(2:sp),'ko-.') %绘制出弹塑性位移延性反应谱; hold on; figure(2) %%%输出数据到硬盘; save('e:\north_concrete\a20\em_a20_3.txt','rm','-ASCII') plot(per,rm,'ko-.') %绘制出累积滞回耗能谱; hold on; elseif(iy==4) figure(1) plot(per(2:sp),um(2:sp),'kp-.') %绘制出弹塑性位移延性反应谱; hold on; figure(2) %%%输出数据到硬盘; save('e:\north_concrete\a20\em_a20_4.txt','rm','-ASCII') plot(per,rm,'kp-.') %绘制出累积滞回耗能谱; hold on; else figure(1) plot(per(2:sp),um(2:sp),'ks-.') %绘制出弹塑性位移延性反应谱; hold on; figure(2) %%%输出数据到硬盘; save('e:\north_concrete\a20\em_a20_5.txt','rm','-ASCII') plot(per,rm,'ks-.') %绘制出累积滞回耗能谱; hold on; end end