求弹塑性时程分析程序（MATLAB）
本帖最后由 hub115 于 2011-12-7 14:44 编辑如题
徐赵东《matlab语言在建筑工程中的应用中》一书中油一段弹塑性时程分析程序，本人条了很久也没有调通，不知道有没有人调通过的，发出来大家借鉴一下！
急切希望大家帮忙看看，这个东西很多人能用到，调通以后对我是个帮助，对其他想用的人也是帮助。 楼主你好，我也看了这本书，不过弹塑性那段的程序没有看懂，不知楼主可否注解一下，我倒是能把他的程序调出来，下面给出他的 两个程序
楼主可将我的程序中的elcen.dat换成你的EI centro.dat这样就可以运行了弹性%%%%%%%%%%%%% MAIN PROGRAM %%%%%%%%%%%%%
% response analysis of structure--- elastic time history analysis method
% structure parameter
tic;
h=; % 各层层高, mm
m=*1e+3;
k0=*1e+5;
cn=length(m); % 总层数
% earthquake parameter
load elcen.dat
dzhbo=elcen(:,2); %dzhbo指地震波
ct=1.4; % Wilson-θ法的θ值
dt=0.01;
xs=200/max(abs(dzhbo)); % 调整地震输入加速度幅值
ag=dzhbo*0.01*xs; % 地震波, 400 步, 步长为0.02s, 单位m/s2
ndzh=2000;
ag1=ag(1:ndzh);
ag2=ag(2:ndzh+1);
agtao=ct*(ag2-ag1);
% initial value
chsh=zeros(cn,1);
wyi1=chsh; %位移
sdu1=chsh; %速度
jsdu1=chsh;%加速度
wyimt=chsh;
sdumt=chsh;
jsdumt=chsh;
unit=ones(cn,1);
m=diag(m);
=matrixju(k0,cn);%子程序形成刚度矩阵
=eig(ik,m);
d=sqrt(d);
w=sort(diag(d)); %将频率从大到小排列
a=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);%认为是瑞利阻尼，求系数
b=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);
c0=a*m+b*ik;
for i=1:ndzh
kxin=ik+(3/(ct*dt))*c0+(6/(ct*ct*dt*dt))*m; %kxin为新的刚度
dpxin=-m*unit*agtao(i)+m*(6/(ct*dt)*sdu1+3*jsdu1)+c0*(3*sdu1+ct*dt/2*jsdu1); %新的力增量
dxtao=kxin\dpxin;
dtjsdu=6*dxtao/(ct*(ct^2*dt^2))-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1;
jsdu=jsdu1+dtjsdu;
dtsdu=(dt/2)*(jsdu+jsdu1);
sdu=sdu1+dtsdu;
dtwyi=dt*sdu1+(1/3)*dt^2*jsdu1+(dt^2/6)*jsdu;
wyi=wyi1+dtwyi;
jsdu=-m\(m*unit*ag2(i)+c0*sdu+ik*wyi); % 调整加速度
wyi1=wyi;
sdu1=sdu;
jsdu1=jsdu;
wyimt=;
sdumt=;
jsdumt=;
end
mm=toc;
mm
t=0:dt:ndzh*dt;
figure(2)
subplot(2,2,1)
plot(t,wyimt(3,:)/1000,'r-')
subplot(2,2,2)
plot(t,jsdumt(3,:),'r-')
弹塑性
%%%%%%%%% MAIN PROGRAM %%%%%%%%%
%(made in July-October, 2000 by Xu Zhao-dong)
%(revised in January 6--16, 2002)%
% this program is three-poly stiffness degeneration model
% this program is the reinforced concrete structure
% structure parameters
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 elcen.dat
dzhbo=elcen(:,2);
ct=1.4;
dt=0.02;
ndzh=400;
tao=ct*dt;
xs=400/max(abs(dzhbo));
ag=xs*0.01*dzhbo;
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); % 求解结构的参数
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(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);
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 % 1:cn
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; % new acceleration
dtsdu=(dt/2)*(jsdu+jsdu1);
sdu=sdu1+dtsdu; % new velocity
dtwyi=dt*sdu1+(dt^2/3)*jsdu1+(dt^2/6)*jsdu;
wyi=wyi1+dtwyi; % new displacement
=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 % 1:ndzh
t=0:0.02:8;
subplot(2,2,1)
plot(t,wyimt(3,:),'k');
subplot(2,2,2)
plot(t,jsdumt(3,:),'k'); 不好意思，不知道怎么传图片，你运行下就可以看到了
可以了，红色为弹性的，下面的为弹塑性的
THX FOR YOUR SHARING 想看看，我是菜鸟
页:
[1]