matlab计算三层混凝土框架结构的地震波弹塑性时程分析

wypzf_8 · 发表于 2013-3-3 19:47

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

本帖最后由 wypzf_8 于 2013-3-3 20:32 编辑

您好，本贴介绍matlab计算三层混凝土框架结构的地震波弹塑性时程分析，资料来源于《土木工程常用软件分析与应用+MATLAB+SAP2000》，关于三层混凝土框架的信息，本贴不在具体介绍，具体详见附件pdf文档，相应的m文件我也附上，此程序虽能计算出结果，但我对结果持怀疑态度，主要是三折线的滞回规则没有完整的体现，在弹塑性阶段，程序的逻辑判断正确，但是到了状态数为5到6,6到3时，我发现不符合规则。如在400ga的Ell波时，一层的恢复力曲线为
400gal一层恢复力模型.jpg

希望高手能够指点迷津，关键问题是程序逻辑判断没有发现错误。
尝试尝试把400更改为600、800，会发现三折线规则没正确体现。
m文件见附件

三折线弹塑性时程分析.rar (20.2 KB, 下载次数: 23)
欢迎交流QQ：664652476 ，邮箱wypzf_8@163.com

wypzf_8 · 发表于 2013-3-3 20:26

主程序

% 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=[1.690 1.581 1.412]*1e+4;
k0=[1.512 2.012 2.012]*1e+7;
h=[4000 7300 10600];
k3stif=[k0' k0'*0.4 k0'*0.1];
cn=length(m0);
m=diag(m0);
xc=[6.3 4.9 4.2]*1e-3; %层间开裂位移
xy=[21.8 18.9 17.2]*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);
[k4stif,fp,fn]=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;
[kju]=matrixju(kk,cn);
[c0]=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
[wyi1]=cjzhuanhuan(xcjtp,cn);
[sdu1]=cjzhuanhuan(sducjtp,cn);
[jsdu1]=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;
[xcj]=zhuanhuan(wyi,cn);
[sducj]=zhuanhuan(sdu,cn);
[jsducj]=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=[pdd pd]; %阶段符号矩阵
[kju]=matrixju(kk,cn); %刚度矩阵的聚合
[c0]=strdamp(m,kju,kxi1,kxi2,cn); %求解结构的阻尼矩阵
plastic2=[plastic(2:cn);0];
plastic3=plastic-plastic2;
plasticmt=[plasticmt plastic3];
dtplastic=plastic3-plasticmt(:,i);
[wyi1]=cjzhuanhuan(xcj,cn);
[sdu1]=cjzhuanhuan(sducj,cn);
[jsdu1]=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; %新位移
[xcj]=zhuanhuan(wyi,cn);
[sducj]=zhuanhuan(sdu,cn);
jsdu=-m\(m*unit*ag2(i)+c0*sdu+kju*wyi+plastic3);%加速度修正
[jsducj]=zhuanhuan(jsdu,cn);
fcj=kk.*xcj+plastic;
xcjq=[xcjq xcj];
sducjq=[sducjq sducj];
jsducjq=[jsducjq jsducj];
fcjq=[fcjq fcj];
wyimt=[wyimt wyi*1000];
sdumt=[sdumt sdu];
jsdumt=[jsdumt jsdu];
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('三层的恢复力模型')

复制代码

wypzf_8 · 发表于 2013-3-3 20:27

子程序zhuanhuan

%子程序：相对于地面的反应值转换成层间相对反应值
%change absolute value into inter-story value
function[inavx]=zhuanhuan(avx,cn)
inavx(1)=avx(1);
for i=2:cn
inavx(i)=avx(i)-avx(i-1);
end
inavx=inavx';

复制代码

wypzf_8 · 发表于 2013-3-3 20:28

子程序strpara

%caculate the fourth stiffness and initial value of 'fp and fn"
function[k4stif,fp,fn]=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=[k3stif k4th];%%k4stif的第1，2，4，列分别对应3线性的 k1,k2,k4.

复制代码

wypzf_8 · 发表于 2013-3-3 20:30

子程序strdamp

%子程序：求解结构的阻尼矩阵
% solve damping of structure
function [c0]=strdamp(m,kju,kxi1,kxi2,cn)
[zhx,www]=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;

复制代码

wypzf_8 · 发表于 2013-3-3 20:31

子程序matrixju

%子程序：形成刚度矩阵或阻尼矩阵的聚合矩阵
% matrix aggregation of system
function[kcju]=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);

复制代码

wypzf_8 · 发表于 2013-3-3 20:32

子程序cjzhuanhuan

%子程序;层间相对反应值转换成相对地面的反应值
%inter-story value change into absolute value
function[avx]=cjzhuanhuan(inavx,cn)
avx(1)=inavx(1);
for i=2:cn
avx(i)=inavx(i)+avx(i-1);
end
avx=avx';

复制代码

wypzf_8 · 发表于 2013-3-3 20:43

关于时程分析的资料，见本楼附件

xiaobao200 · 发表于 2013-8-24 14:49

谢谢楼主分享~

怪搞悠嘻 · 发表于 2013-9-29 14:54

谢谢分享
！！！！！

hxh3956537 · 发表于 2014-6-10 16:47

看不到pdf

chybeyond · 发表于 2014-6-10 16:58

hxh3956537 发表于 2014-6-10 16:47
看不到pdf

http://forum.vibunion.com/thread-132101-1-1.html

hxh3956537 · 发表于 2014-6-10 17:03

求楼主放PDF啊

hxh3956537 · 发表于 2014-6-10 17:22

chybeyond 发表于 2014-6-10 16:58
http://forum.vibunion.com/thread-132101-1-1.html

多谢了！做任务去了

yangfengdehaizi · 发表于 2014-10-8 09:32

权限不够，PDF看不到啊。

账号		自动登录	找回密码
密码			我要加入

[编程技巧] matlab计算三层混凝土框架结构的地震波弹塑性时程分析

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐: