concrete_trilinear.m

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

function [kk2,x2,f2,pd,ucom,up,fp,un,fn]=concrete_trilinear(pd,per,nt,m,eta,kk1,ke,kc,ky,uc,uy,dt,xg,x1,x2,f1,ucom,up,fp,un,fn)
%退化三线型恢复力模型，适用于混凝土结构非线性地震分析；
%钢混开裂前后分别为ke、kc；钢混屈服后分别为ky； 
kk2=kk1;
f2=f1+kk1*(x2(1)-x1(1));
fc=ke*uc;fy=fc+kc*(uy-uc);
%采用狄龙示意图编程；
if pd==0 
    if x2(1)>uc
        pd=2; kk2=kc;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tc,xc]=guaidian1(dt,ke,m,eta,xg,x1,uc);
        txg(1)=xg(1)+(xg(2)-xg(1))*tc/dt;
		txg(2)=xg(2);
		dt2=dt-tc;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xc);  
        if(x2(1)>uy)
            x2(1)=uy;
        elseif(x2(1)<uc)
            x2(1)=uc;
        end
        if(x2(2)<0)
            x2(2)=0;% simplify;
        end
        f2=fc+kc*(x2(1)-uc);
        if (x2(1)>ucom)%%%the crack point;
            ucom=x2(1);
        end
    elseif x2(1)<-uc
        pd=-2; kk2=kc;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tc,xc]=guaidian1(dt,ke,m,eta,xg,x1,-uc);
        txg(1)=xg(1)+(xg(2)-xg(1))*tc/dt;
		txg(2)=xg(2);
		dt2=dt-tc;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xc);  
        if(x2(1)<-uy)
            x2(1)=-uy;
        elseif(x2(1)>-uc)
            x2(1)=-uc;
        end
        if(x2(2)>0)
            x2(2)=0;% simplify;
        end
        f2=-fc+kc*(x2(1)+uc);
        if (abs(x2(1))>ucom)%%%the crack point;
            ucom=abs(x2(1));
        end
    end    
elseif pd==2 
    if(x2(1)==uy&x2(2)>0)
        pd=3;kk2=ky;
    elseif(x2(1)>uy )
        pd=3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[ty,xy]=guaidian1(dt,kk1,m,eta,xg,x1,uy);
        txg(1)=xg(1)+(xg(2)-xg(1))*ty/dt;
		txg(2)=xg(2);
		dt2=dt-ty;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xy);  
        if(x2(1)<uy)
            x2(1)=uy;
%           刚度和状态值不变；
        end
        if(x2(2)<0)
            x2(2)=0;% simplify;
        end
        f2=fy+kk2*(x2(1)-uy);      
    elseif (x2(2)<0) 
        pd=1; %%%kk2后面处理；
        %第二类拐点处理；采用精确拐点非迭代法；
        if(x1(2)==0)
            td=0;xd=x1;
        else
			[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1);
        end
        if xd(1)>uy
            xd(1)=uy;%近似处理；骨架曲线原本粗糙，无必要太精确；
        end
        fd=f1+kk1*(xd(1)-x1(1));
        %设立最值点，确定屈服前卸载刚度；
        if xd(1)>ucom
            ucom=xd(1);
        end
        kk2=fd/xd(1);
        txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
        txg(2)=xg(2);
        dt2=dt-td;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xd);  
        if x2(1)>xd(1)
            x2(1)=xd(1)
        end
        f2=fd+kk2*(x2(1)-xd(1));
    end 
elseif pd==-2
    if(x2(1)==-uy&x2(2)<0)
        pd=-3;kk2=ky;
    elseif(x2(1)<-uy)
        pd=-3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[ty,xy]=guaidian1(dt,kk1,m,eta,xg,x1,-uy);
        txg(1)=xg(1)+(xg(2)-xg(1))*ty/dt;
		txg(2)=xg(2);
		dt2=dt-ty;
        [x2]=newmark_single(dt2,ky,m,eta,txg,xy);  
        if(x2(1)>-uy)
            x2(1)=-uy;
%           刚度和状态值不变；
        end
        if(x2(2)>0)
            x2(2)=0;% simplify;
        end        
        f2=-fy+kk2*(x2(1)+uy);
	elseif (x2(2)>0) 
		pd=1; 
        if(nt==280&per==1.2)
            bug=1%%%%好奇怪，就是读不出per=1.2，matlab有bug;
        end
		%第二类拐点处理；采用精确拐点非迭代法；
        if(x1(2)==0)
            td=0;xd=x1;
        else
			[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1);
        end
		if xd(1)<-uy
			xd(1)=-uy;%近似处理；骨架曲线原本粗糙，无必要太精确；
		end
		fd=f1+kk1*(xd(1)-x1(1));
		%设立最值点，确定屈服前卸载刚度；
        if abs(xd(1))>ucom
            ucom=abs(xd(1));
        end
		kk2=fd/xd(1);
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);     
        if x2(1)<xd(1)
            x2(1)=xd(1)
        end
		f2=fd+kk2*(x2(1)-xd(1));
    end   
%%%%%
elseif pd==1 
    if x2(1)>ucom
        pd=2;kk2=kc;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tcom,xcom]=guaidian1(dt,kk1,m,eta,xg,x1,ucom);
        txg(1)=xg(1)+(xg(2)-xg(1))*tcom/dt;
		txg(2)=xg(2);
		dt2=dt-tcom;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xcom);  
        if(x2(1)>uy)
            x2(1)=uy;
            if(x2(2)>=0)
                pd=3;kk2=ky;
            else
                pd=4;kk2=ke;
            end
        elseif(x2(1)<ucom)
            x2(1)=ucom;
        end
        if(x2(2)<0)
            x2(2)=0;% simplify;
        end
        f2=f1+kk1*(xcom(1)-x1(1))+kk2*(x2(1)-xcom(1));
        if(x2(1)>ucom)
            ucom=x2(1);
        end
    elseif x2(1)<-ucom
        pd=-2;kk2=kc;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tcom,xcom]=guaidian1(dt,kk1,m,eta,xg,x1,-ucom);
        txg(1)=xg(1)+(xg(2)-xg(1))*tcom/dt;
		txg(2)=xg(2);
		dt2=dt-tcom;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xcom);  
        if(x2(1)<-uy)
            x2(1)=-uy;
            if(x2(2)<=0)
                pd=-3;kk2=ky;
            else
                pd=-4;kk2=ke;
            end
        elseif(x2(1)>-ucom)
            x2(1)=-ucom;
        end
        if(x2(2)>0)
            x2(2)=0;% simplify;
        end
        f2=f1+kk1*(xcom(1)-x1(1))+kk2*(x2(1)-xcom(1));
        if(abs(x2(1))>ucom)
            ucom=abs(x2(1));
        end
    end
%%%%%%%%
elseif pd==3
    if x2(2)<0
        pd=4;kk2=ke;
        %第二类拐点处理；采用精确拐点非迭代法；           
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
        fd=f1+kk1*(xd(1)-x1(1));
        if (xd(1)>up)
            up=xd(1);fp=fd;
        end
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd); 
        if (x2(1)>up)
            x2(1)=up;
        end
		f2=fd+kk2*(x2(1)-xd(1)); 
    end
    if (x2(1)>up)%%理论上肯定要大于;
        up=x2(1);fp=f2;
    end
%%%%%%%%%
elseif pd==-3
    if x2(2)>0
        pd=-4;kk2=ke;
        %第二类拐点处理；采用精确拐点非迭代法；
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
		fd=f1+kk1*(xd(1)-x1(1));
        if (xd(1)<un)
            un=xd(1);fn=fd;
        end
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);  
        if (x2(1)<un)
            x2(1)=un;
        end
		f2=fd+kk2*(x2(1)-xd(1));
    end
    if (x2(1)<un)%%理论上肯定要小于;
        un=x2(1);fn=f2;
    end
%%%%%%%%%
elseif pd==4
    if f2<0
        pd=-6;
        %第三类拐点处理；采用精确拐点非迭代法；
        uz=x1(1)-f1/kk1;%求零力处位移；
        kk2=fn/(un-uz);
		[tz,xz]=guaidian1(dt,kk1,m,eta,xg,x1,uz);
        txg(1)=xg(1)+(xg(2)-xg(1))*tz/dt;
		txg(2)=xg(2);
		dt2=dt-tz;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xz);  
        if(x2(1)>uz)
            x2(1)=uz;
        elseif(x2(1)<un)
            x2(1)=un;
        end
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
        f2=kk2*(x2(1)-xz(1));
        if (x2(1)==uz)
            f2=0;
        end
    elseif x2(1)>up
        pd=3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tp,xp]=guaidian1(dt,kk1,m,eta,xg,x1,up);
        txg(1)=xg(1)+(xg(2)-xg(1))*tp/dt;
		txg(2)=xg(2);
		dt2=dt-tp;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xp);  
        if(x2(1)<up)
            x2(1)=up;
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
        f2=f1+kk1*(xp(1)-x1(1))+kk2*(x2(1)-xp(1));
        up=x2(1);fp=f2;
    end
%%%%%%%%%%%%%%
elseif pd==-4
    if f2>0
        pd=6;
        %第三类拐点处理；采用精确拐点非迭代法；
        uz=x1(1)-f1/kk1;%求零力处位移；
        kk2=fp/(up-uz);
		[tz,xz]=guaidian1(dt,kk1,m,eta,xg,x1,uz);
        txg(1)=xg(1)+(xg(2)-xg(1))*tz/dt;
		txg(2)=xg(2);
		dt2=dt-tz;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xz);  
        if(x2(1)<uz)
            x2(1)=uz;
        elseif(x2(1)>up)
            x2(1)=up;
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
        f2=kk2*(x2(1)-xz(1));
        if (x2(1)==uz)
            f2=0;
        end
    elseif x2(1)<un
        pd=-3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tn,xn]=guaidian1(dt,kk1,m,eta,xg,x1,un);
        txg(1)=xg(1)+(xg(2)-xg(1))*tn/dt;
		txg(2)=xg(2);
		dt2=dt-tn;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xn);  
        if(x2(1)>un)
            x2(1)=un;
        end
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
        f2=f1+kk1*(xn(1)-x1(1))+kk2*(x2(1)-xn(1));
        un=x2(1);fn=f2;
    end
%%%%%%%%%%%%%%
elseif pd==-6
    if (x2(1)<un)
        pd=-3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tn,xn]=guaidian1(dt,kk1,m,eta,xg,x1,un);
        txg(1)=xg(1)+(xg(2)-xg(1))*tn/dt;
		txg(2)=xg(2);
		dt2=dt-tn;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xn);  
        if(x2(1)>un)
            x2(1)=un;
        end
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
        f2=f1+kk1*(xn(1)-x1(1))+kk2*(x2(1)-xn(1));
        un=x2(1); fn=f2;
    elseif (x2(2)>0)
        pd=-5;kk2=ke;
        %第二类拐点处理；采用精确拐点非迭代法；
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
        if(xd(1)<un)
            xd(1)=un;%近似处理；骨架曲线原本粗糙，无必要太精确；
        end
		fd=f1+kk1*(xd(1)-x1(1));
        uz=xd(1)-fd/kk2;%求零力处位移；
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);   
        if(x2(1)<un)
            x2(1)=un;%近似处理；骨架曲线原本粗糙，无必要太精确；
        elseif(x2(1)>uz)
            x2(1)=uz;
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
		f2=fd+kk2*(x2(1)-xd(1));
        if (x2(1)==uz)
            f2=0;
        end
    end
%%%%%%%%%%%%%%
elseif pd==6
    if (x2(1)>up)
        pd=3;kk2=ky;
        %第一类拐点处理；采用精确拐点非迭代法；
		[tp,xp]=guaidian1(dt,kk1,m,eta,xg,x1,up);
        txg(1)=xg(1)+(xg(2)-xg(1))*tp/dt;
		txg(2)=xg(2);
		dt2=dt-tp;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xp);  
        if(x2(1)<up)
            x2(1)=up;
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
        f2=f1+kk1*(xp(1)-x1(1))+kk2*(x2(1)-xp(1));
        up=x2(1); fp=f2;
    elseif (x2(2)<0)
        pd=5;kk2=ke;
        %第二类拐点处理；采用精确拐点非迭代法；
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
        if(xd(1)>up)
            xd(1)=up;%近似处理；骨架曲线原本粗糙，无必要太精确；
        end
		fd=f1+kk1*(xd(1)-x1(1));
        uz=xd(1)-fd/kk2;%求零力处位移；
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);   
        if(x2(1)>up)
            x2(1)=up;%近似处理；骨架曲线原本粗糙，无必要太精确；
        elseif(x2(1)<uz)
            x2(1)=uz;
        end  
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
		f2=fd+kk2*(x2(1)-xd(1));
        if (x2(1)==uz)
            f2=0;
        end
    end
%%%%%%%%%%%%%%
elseif pd==-5    
    if (f2>0)
        pd=6;
        %第三类拐点处理；采用精确拐点非迭代法；
        uz=x1(1)-f1/kk1;%求零力处位移；
        kk2=fp/(up-uz);
		[tz,xz]=guaidian1(dt,kk1,m,eta,xg,x1,uz);
        txg(1)=xg(1)+(xg(2)-xg(1))*tz/dt;
		txg(2)=xg(2);
		dt2=dt-tz;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xz);  
        if(x2(1)<uz)
            x2(1)=uz;
        elseif(x2(1)>up)
            x2(1)=up;
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
        f2=kk2*(x2(1)-xz(1));
    elseif (x2(2)<0)
        pd=-6;
        %第二类拐点处理；采用精确拐点非迭代法；
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
        uz=x1(1)-f1/kk1;
        if(xd(1)>uz)
            xd(1)=uz;%近似处理；骨架曲线原本粗糙，无必要太精确；
        end
		fd=f1+kk1*(xd(1)-x1(1));
        kk2=(fn-fd)/(un-xd(1));
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);   
        if(x2(1)<un)
            x2(1)=un;%近似处理；骨架曲线原本粗糙，无必要太精确；
        elseif(x2(1)>xd(1))
            x2(1)=xd(1);
        end
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
		f2=fd+kk2*(x2(1)-xd(1));
    end
%%%%%%%%%%%%%%
elseif pd==5    
    if (f2<0)
        pd=-6;
        %第三类拐点处理；采用精确拐点非迭代法；
        uz=x1(1)-f1/kk1;%求零力处位移；
        kk2=fn/(un-uz);
		[tz,xz]=guaidian1(dt,kk1,m,eta,xg,x1,uz);
        txg(1)=xg(1)+(xg(2)-xg(1))*tz/dt;
		txg(2)=xg(2);
		dt2=dt-tz;
        [x2]=newmark_single(dt2,kk2,m,eta,txg,xz);  
        if(x2(1)>uz)
            x2(1)=uz;
        elseif(x2(1)<un)
            x2(1)=un;
        end
        if(x2(2)>0)
            x2(2)=0;%%%simplify;
        end
        f2=kk2*(x2(1)-xz(1));
    elseif (x2(2)>0)
        pd=6;
        %第二类拐点处理；采用精确拐点非迭代法；
		[td,xd]=guaidian2(dt,kk1,m,eta,xg,x1); 
        uz=x1(1)-f1/kk1;
        if(xd(1)<uz)
            xd(1)=uz;%近似处理；骨架曲线原本粗糙，无必要太精确；
        end
		fd=f1+kk1*(xd(1)-x1(1));
        kk2=(fp-fd)/(up-xd(1));
		txg(1)=xg(1)+(xg(2)-xg(1))*td/dt;
		txg(2)=xg(2);
		dt2=dt-td;
		[x2]=newmark_single(dt2,kk2,m,eta,txg,xd);   
        if(x2(1)>up)
            x2(1)=up;%近似处理；骨架曲线原本粗糙，无必要太精确；
        elseif(x2(1)<xd(1))
            x2(1)=xd(1);
        end
        if(x2(2)<0)
            x2(2)=0;%%%simplify;
        end
		f2=fd+kk2*(x2(1)-xd(1));
    end
end