nonlinear.m

钢、混凝土结构时频反应谱程序

function [yu,d,v,a,ep]=nonlinear(para,sa,np,per,xg,t,dt,eb)
% 能量分析则分别有ei、ek、ed、ep、ee；
st=length(xg);
m=1;%先假定为单位质量；
k1=m*(2*pi/per)^2;%初始弹性刚度，计算时频反应谱是一个系列；
%此为恢复力模型输入参数；
eta=0.05;
afa=0.02;%%钢混一般取0.02，理想弹塑性为0,可进行调整；
k0=k1/2;k3=k0*afa;%假定初始刚度为k0的2～5倍；
%假定px/py=bta，一般为0.6～0.7；由此确定k2；
bta=0.8;%%%%%py=px/bta;
k2=k1*(1-bta)/(k1/k0-bta);
%对于周期较大的工况，对应屈服位移也应较大，因为结构柔度大；
ebuxi=eb;
py=ebuxi*sa*m;%一般工程结构屈服强度系数ebuxi为0.25~0.4左右；
px=py*bta;
ux=px/k1;
uy=py/k0;
f1=0;
kk1=k1;%迭代前的刚度为k1，此处将初始弹性刚度赋予它；
txg=zeros(1,2);
x1=zeros(1,3);
x2=zeros(1,3);   
x1(3)=-xg(1);
nst=0;
for nt=1:st   
	if (nt==1)%初始时刻位移和速度为零，加速度为负初始时刻地震加速度；
        x(1)=0;
        v(1)=0;
        a(1)=-xg(1);
        ff(1)=0;
        ep=0;%能量分析值；
	else
		txg(1)=xg(nt-1);
		txg(2)=xg(nt);
        [x2]=newmark_single(dt,kk1,m,eta,txg,x1);
        [kk2,x2,f2,nst]=kb_trilinear(per,nt,m,eta,kk1,k1,k2,k3,ux,uy,dt,txg,x1,x2,f1,nst);
        %能量分析值，采用近似梯形公式，最好是考虑拐点分割为精确；
        if (f1*f2>=0)
            ep=ep+(f2+f1)*(x2(1)-x1(1));  
        else
            ep=ep-f1^2/kk1+f2^2/kk2;
        end
        kk(nt)=kk2;
	    kk1=kk2;
	    f1=f2;
	    x1=x2;
		d(nt)=x2(1);
		v(nt)=x2(2);
		a(nt)=x2(3);
    end
    ff(nt)=f1;
end
yu=abs(d)/uy;
ee=f2^2/kk2/2;
ep=ep-ee;




% if (per==0.2)
%     for ni=1:500
% 	figure(31)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
%     for ni=501:1000
% 	figure(32)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
%     for ni=1001:st
% 	figure(33)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
% end
% if (per==0.2)
%     for ni=1:300
% 	figure(41)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
%     for ni=301:400
% 	figure(42)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
%     for ni=401:500
% 	figure(43)
% 	scatter(d(ni),ff(ni));
% 	hold on;
% 	end
% end