matlab计算三层混凝土框架结构的地震波弹塑性时程分析
本帖最后由 wypzf_8 于 2013-3-3 20:32 编辑您好，本贴介绍matlab计算三层混凝土框架结构的地震波弹塑性时程分析，资料来源于《土木工程常用软件分析与应用+MATLAB+SAP2000》，关于三层混凝土框架的信息，本贴不在具体介绍，具体详见附件pdf文档，相应的m文件我也附上，此程序虽能计算出结果，但我对结果持怀疑态度，主要是三折线的滞回规则没有完整的体现，在弹塑性阶段，程序的逻辑判断正确，但是到了状态数为5到6,6到3时，我发现不符合规则。如在400ga的Ell波时，一层的恢复力曲线为
希望高手能够指点迷津，关键问题是程序逻辑判断没有发现错误。
尝试尝试把400更改为600、800，会发现三折线规则没正确体现。
m文件见附件
欢迎交流QQ：664652476 ，邮箱wypzf_8@163.com
主程序% This is a program for inelastic analysis of 3-layers structure
% structure parameters
disp('--- Welcome to use this program ! ---');
disp('--- This is a program for inelastic analysis of 3-layers structure ! ---');
clear,close all,
m0=*1e+4;
k0=*1e+7;
h=;
k3stif=;
cn=length(m0);
m=diag(m0);
xc=*1e-3;%层间开裂位移
xy=*1e-3;%层间屈服位移
kxi1=0.05;kxi2=0.07;%阻尼比
%earthquake parameters
load elcentroEW.m;
seismicwaves=elcentroEW(:,2);
t=elcentroEW(:,1);
ct=1.4;% in the wilson-ct method
dt=0.02;
ndzh=length(elcentroEW(:,2))-1;
tao=ct*dt;
xs=400/max(abs(seismicwaves)); %尝试把400更改为600、800，发现三折线规则没正确体现
ag=xs*0.01*seismicwaves;
ag1=ag(1:ndzh);
ag2=ag(2:ndzh+1);
agtao=ct*(ag2-ag1);
%initial value
chsh=zeros(cn,1);
pd=chsh;pdd=chsh;
xcj=chsh;sducj=chsh;jsducj=chsh;fcj=chsh;
xcjq=chsh;sducjq=chsh;jsducjq=chsh;fcjq=chsh;
wyi=chsh;sdu=chsh;jsdu=chsh;
wyimt=chsh;sdumt=chsh;jsdumt=chsh;
plastic=chsh;plasticmt=chsh;
unit=ones(cn,1);
xp=xy;xn=-xy;
kk=k3stif(:,1);
=strpara(cn,k3stif,xc,xy);%solve the structure parameters
for i=1:ndzh %对地震波点数的循环计算
pdt=zeros(cn,1);
for j=1:cn %对结构自由度的循环计算
if pd(j)==0
kk(j)=k4stif(j,1);
plastic(j)=0;
fcj(j)=kk(j)*xcj(j)+plastic(j);
if xcj(j)>xc(j)
pd(j)=2;
pdt(j)=1;
bl=(xc(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,2);
plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)<-xc(j)
pd(j)=-2;
pdt(j)=1;
bl=(-xc(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,2);
plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==2
kk(j)=k4stif(j,2);
plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if xcj(j)>xy(j)
pd(j)=3;
pdt(j)=1;
bl=(xy(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif (xcj(j)<xy(j))&&(sducj(j)<0)
x1a(j)=abs(xcj(j));
pd(j)=1;
kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
plastic(j)=0;
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==-2
kk(j)=k4stif(j,2);
plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if xcj(j)<-xy(j)
pd(j)=-3;
pdt(j)=1;
bl=(-xy(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif (xcj(j)>-xy(j))&&(sducj(j)>0);
x1a(j)=abs(xcj(j));
pd(j)=1;
kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
plastic(j)=0;
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==1
kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
plastic(j)=0;
fcj(j)=kk(j)*xcj(j)+plastic(j);
if xcj(j)>x1a(j)
pd(j)=2;
pdt(j)=1;
bl=(x1a(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,2);
plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)<-x1a(j)
pd(j)=-2;
pdt(j)=1;
bl=(-x1a(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,2);
plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==3
kk(j)=k4stif(j,3);
plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if sducj(j)<0
xp(j)=xcj(j);
fp(j)=fcj(j);
pd(j)=4;
kk(j)=k4stif(j,4);
plastic(j)=fp(j)-kk(j)*xp(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==-3
kk(j)=k4stif(j,3);
plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if sducj(j)>0
xn(j)=xcj(j);
fn(j)=fcj(j);
pd(j)=-4;
kk(j)=k4stif(j,4);
plastic(j)=fn(j)-kk(j)*xn(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==4
kk(j)=k4stif(j,4);
plastic(j)=fp(j)-kk(j)*xp(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if fcj(j)<0
pd(j)=-6;
pdt(j)=1;
bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
xa(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
fa(j)=0;
kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
plastic(j)=fa(j)-kk(j)*xa(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)>xp(j)
pd(j)=3;
pdt(j)=1;
bl=(xp(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==-4
kk(j)=k4stif(j,4);
plastic(j)=fn(j)-kk(j)*xn(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if fcj(j)>0
pd(j)=6;
pdt(j)=1;
bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
xb(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
fb(j)=0;
kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
plastic(j)=fb(j)-kk(j)*xb(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)<xn(j)
pd(j)=-3;
pdt(j)=1;
bl=(xn(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==-6
kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
plastic(j)=fa(j)-kk(j)*xa(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if sducj(j)>0
xe(j)=xcj(j);
fe(j)=fcj(j);
pd(j)=-5;
kk(j)=k4stif(j,4);
plastic(j)=fe(j)-kk(j)*xe(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)<xn(j)
pd(j)=-3;
pdt(j)=1;
bl=(xn(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd==6
kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
plastic(j)=fb(j)-kk(j)*xb(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if sducj(j)<0
xf(j)=xcj(j);
ff(j)=fcj(j);
pd(j)=5;
kk(j)=k4stif(j,4);
plastic(j)=ff(j)-kk(j)*xf(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif xcj(j)>xp(j)
pd(j)=3;
pdt(j)=1;
bl=(xp(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
kk(j)=k4stif(j,3);
plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==-5
kk(j)=k4stif(j,4);
plastic(j)=fe(j)-kk(j)*xe(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if fcj(j)>0
pd(j)=6;
pdt(j)=1;
bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
xb(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
fb(j)=0;
kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
plastic(j)=fb(j)-kk(j)*xb(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif (sducj(j)<0)&&(fcj(j)<=0)
xa(j)=xcj(j);
fa(j)=fcj(j);
pd(j)=-6;
kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
plastic(j)=fa(j)-kk(j)*xa(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pd(j)==5
kk(j)=k4stif(j,4);
plastic(j)=ff(j)-kk(j)*xf(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
if fcj(j)<0
pd(j)=-6;
pdt(j)=1;
bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
xa(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
fa(j)=0;
kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
plastic(j)=fa(j)-kk(j)*xa(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
elseif (sducj(j)>0)&&(fcj(j)>=0)
xb(j)=xcj(j);
fb(j)=fcj(j);
pd(j)=6;
kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
plastic(j)=fb(j)-kk(j)*xb(j);
fcj(j)=kk(j)*xcj(j)+plastic(j);
end
end
if pdt(j)==1
dt1=(1-bl)*dt;
tao1=ct*dt1;
=matrixju(kk,cn);
=strdamp(m,kju,kxi1,kxi2,cn);
dtp1=(1-bl)*dtplastic;
for s=1:cn
xcjtp(s)=xcjq(s,i-1)+bl*(xcj(s)-xcjq(s,i-1));
sducjtp(s)=sducjq(s,i-1)+bl*(sducj(s)-sducjq(s,i-1));
jsducjtp(s)=jsducjq(s,i-1)+bl*(jsducj(s)-jsducjq(s,i-1));
end
=cjzhuanhuan(xcjtp,cn);
=cjzhuanhuan(sducjtp,cn);
=cjzhuanhuan(jsducjtp,cn);
kxin1=kju+(3/tao1)*c0+(6/tao1^2)*m;
dpxin1=-m*unit*agtao(i-1)*(1-bl)+m*(6/(tao1)*sdu1+3*jsdu1)+c0*(3*sdu1+tao1/2*jsdu1)-ct*dtp1;
dxtao1=kxin1\dpxin1;
dtjsdu1=6*dxtao1/(ct*(tao1^2))-6*sdu1/(ct*tao1)-(3/ct)*jsdu1;
jsdu=jsdu1+dtjsdu1;
dtsdu1=(dt1/2)*(jsdu+jsdu1);
sdu=sdu1+dtsdu1;
dtwyi1=dt1*sdu1+(dt1^2/3)*jsdu1+(dt1^2/6)*jsdu;
wyi=wyi1+dtwyi1;
=zhuanhuan(wyi,cn);
=zhuanhuan(sdu,cn);
=zhuanhuan(jsdu,cn);
xcjq(:,i)=xcj;
sducjq(:,i)=sducj;
jsducjq(:,i)=jsducj;
fcjq(:,i)=fcj;
wyimt(:,i)=wyi*1000;
sdumt(:,i)=sdu;
jsdumt(:,i)=jsdu;
end
end
pdd=; %阶段符号矩阵
=matrixju(kk,cn); %刚度矩阵的聚合
=strdamp(m,kju,kxi1,kxi2,cn);%求解结构的阻尼矩阵
plastic2=;
plastic3=plastic-plastic2;
plasticmt=;
dtplastic=plastic3-plasticmt(:,i);
=cjzhuanhuan(xcj,cn);
=cjzhuanhuan(sducj,cn);
=cjzhuanhuan(jsducj,cn);
kxin=kju+(3/tao)*c0+(6/tao^2)*m;
dpxin=-m*unit*agtao(i)+m*(6/tao*sdu1+3*jsdu1)+c0*(3*sdu1+tao/2*jsdu1)-ct*dtplastic;
dxtao=kxin\dpxin;
dtjsdu=6*dxtao/(ct*(tao^2))-6*sdu1/(ct*tao)-(3/ct)*jsdu1;
jsdu=jsdu1+dtjsdu;%新加速度
dtsdu=(dt/2)*(jsdu+jsdu1);
sdu=sdu1+dtsdu;%新速度
dtwyi=dt*sdu1+(dt^2/3)*jsdu1+(dt^2/6)*jsdu;
wyi=wyi1+dtwyi;%新位移
=zhuanhuan(wyi,cn);
=zhuanhuan(sdu,cn);
jsdu=-m\(m*unit*ag2(i)+c0*sdu+kju*wyi+plastic3);%加速度修正
=zhuanhuan(jsdu,cn);
fcj=kk.*xcj+plastic;
xcjq=;
sducjq=;
jsducjq=;
fcjq=;
wyimt=;
sdumt=;
jsdumt=;
end
figure(1)
subplot(3,1,1);
plot(t,jsdumt(1,:),'k');
subplot(3,1,2);
plot(t,jsdumt(2,:),'k');
subplot(3,1,3);
plot(t,jsdumt(3,:),'k');
figure(2)
plot(xcjq(1,:),fcjq(1,:));%%%相对位移与层间剪力的关系图形。
grid on
title('一层的恢复力模型')
figure(3)
plot(xcjq(2,:),fcjq(2,:));%%%相对位移与层间剪力的关系图形。
grid on
title('二层的恢复力模型')
figure(4)
plot(xcjq(3,:),fcjq(3,:));%%%相对位移与层间剪力的关系图形。
grid on
title('三层的恢复力模型')
子程序zhuanhuan%子程序：相对于地面的反应值转换成层间相对反应值
%change absolute value into inter-story value
function=zhuanhuan(avx,cn)
inavx(1)=avx(1);
for i=2:cn
inavx(i)=avx(i)-avx(i-1);
end
inavx=inavx'; 子程序strpara%caculate the fourth stiffness and initial value of 'fp and fn"
function=strpara(cn,k3stif,xc,xy)
fp=zeros(cn,1);
fn=zeros(cn,1);
for i=1:cn
k4th(i,1)=(k3stif(i,1)*xc(i)+k3stif(i,2)*(xy(i)-xc(i)))/xy(i);%%塑性时卸载的刚度k4。
fp(i)=k3stif(i,1)*xc(i)+k3stif(i,2)*(xy(i)-xc(i)); %%达到最大弹性时所受的力。
fn(i)=-fp(i);
end
k4stif=;%%k4stif的第1，2，4，列分别对应3线性的 k1,k2,k4.
子程序strdamp%子程序：求解结构的阻尼矩阵
% solve damping of structure
function =strdamp(m,kju,kxi1,kxi2,cn)
=eig(kju,m);
www=sqrt(www);
w=sort(diag(www));%vibration frequency
a=2*w(1)*w(2)*(kxi1*w(2)-kxi2*w(1))/(w(2)^2-w(1)^2);
b=2*(kxi2*w(2)-kxi1*w(1))/(w(2)^2-w(1)^2);
c0=a*m+b*kju; 子程序matrixju%子程序：形成刚度矩阵或阻尼矩阵的聚合矩阵
% matrix aggregation of system
function=matrixju(korc,cn)
kcju=zeros(cn);
for i=1:cn-1
kcju(i,i)=korc(i)+korc(i+1);
kcju(i,i+1)=-korc(i+1);
kcju(i+1,i)=-korc(i+1);
end
kcju(cn,cn)=korc(cn);
子程序cjzhuanhuan%子程序;层间相对反应值转换成相对地面的反应值
%inter-story value change into absolute value
function=cjzhuanhuan(inavx,cn)
avx(1)=inavx(1);
for i=2:cn
avx(i)=inavx(i)+avx(i-1);
end
avx=avx'; 关于时程分析的资料，见本楼附件{:3_47:} 谢谢楼主分享~ 谢谢分享
！！！！！ 看不到pdf hxh3956537 发表于 2014-6-10 16:47
看不到pdf
http://forum.vibunion.com/thread-132101-1-1.html 求楼主放PDF啊 chybeyond 发表于 2014-6-10 16:58
http://forum.vibunion.com/thread-132101-1-1.html
多谢了！做任务去了 权限不够，PDF看不到啊。
页:
[1]
2