ks_bilinear.m

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

function [kk2,x2,f2]=ks_bilinear(per,nt,m,eta,kk1,k1,k2,uy,dt,xg,x1,x2,f1)
%双线型恢复力模型，适用于钢结构非线性地震分析；
%刚度开裂前为k1，开裂后为k2；拐点处理后xt和f2会改变；
py=k1*uy;
f0=(k1-k2)*uy;
kk2=kk1;
f2=f1+kk1*(x2(1)-x1(1));
%%%%控制精度，重新编写程序，更清晰和实用；拐点处理更明了。
if(kk1==k1)
    xup=(f0+x1(1)*k1-f1)/(k1-k2);
    xdn=xup-2*uy;
    fup=f0+k2*xup;
    fdn=fup-2*py;
    if(f2>fup)%第一类拐点处理；采用精确拐点非迭代法；
        kk2=k2;
        [tq,xq]=guaidian1(dt,k1,m,eta,xg,x1,xup);
        fq=fup;
%         fq=f1+k1*(xq(1)-x1(1));%机器误差，出现问题；
        txg(1)=xg(1)+(xg(2)-xg(1))*tq/dt;
        txg(2)=xg(2);
        dt1=dt-tq;
        [x2]=newmark_single(dt1,k2,m,eta,txg,xq);  
        f2=fq+k2*(x2(1)-xq(1));
        if(x2(2)<0)%采用精确拐点非迭代法求出第二类零速度拐点；
            kk2=k1;
            [tm,xm]=guaidian2(dt1,k2,m,eta,txg,xq);     
%             fm=fq+k2*(xm(1)-xq(1));
            xup=xm(1);
            xdn=xup-2*uy;
            fup=f0+k2*xup;
            fdn=fup-2*py;
            fm=fup;%防止机器误差；
            dt2=dt1-tm;
            [x2]=newmark_single(dt2,k1,m,eta,txg,xm);              
            if(x2(1)>xup)%此处近似处理是否恰当？
                x2(1)=xup;
            elseif(x2(1)<xdn)
                x2(1)=xdn;
            end 
            f2=fm+k1*(x2(1)-xm(1));
%             x2=xm;f2=fm;x2(2)=-1e-9;
        end
    elseif(f2<fdn)
        kk2=k2;
        [tq,xq]=guaidian1(dt,k1,m,eta,xg,x1,xdn);
%         fq=f1+k1*(xq(1)-x1(1));
        fq=fdn;
        txg(1)=xg(1)+(xg(2)-xg(1))*tq/dt;
        txg(2)=xg(2);
        dt1=dt-tq;
        [x2]=newmark_single(dt1,k2,m,eta,txg,xq);  
        f2=fq+k2*(x2(1)-xq(1));
        if(x2(2)>0)%采用精确拐点非迭代法求出第二类零速度拐点；
            kk2=k1;
            [tm,xm]=guaidian2(dt1,k2,m,eta,txg,xq);     
%             fm=fq+k2*(xm(1)-xq(1));
            xdn=xm(1);
            xup=xdn+2*uy;
            fup=f0+k2*xup;
            fdn=fup-2*py;
            fm=fdn;
            dt2=dt1-tm;
            [x2]=newmark_single(dt2,k1,m,eta,txg,xm);              
            if(x2(1)>xup)%此处近似处理是否恰当？
                x2(1)=xup;
            elseif(x2(1)<xdn)
                x2(1)=xdn;
            end 
            f2=fm+k1*(x2(1)-xm(1));
%             x2=xm;f2=fm;x2(2)=1e-9;
        end
    else
    end
elseif(kk1==k2)
    if(x1(2)>0)
         if(x2(2)<0)%采用精确拐点非迭代法求出第二类零速度拐点；
            kk2=k1;
            [tm,xm]=guaidian2(dt,k2,m,eta,xg,x1);     
%             fm=f1+k2*(xm(1)-x1(1));
            xup=xm(1);
            xdn=xup-2*uy;
            fup=f0+k2*xup;
            fdn=fup-2*py;
            fm=fup;
            txg(1)=xg(1)+(xg(2)-xg(1))*tm/dt;
            txg(2)=xg(2);
            dt1=dt-tm;
            [x2]=newmark_single(dt1,k1,m,eta,txg,xm);  
            if(x2(1)>xup)%此处近似处理是否恰当？
                x2(1)=xup;
            elseif(x2(1)<xdn)
                x2(1)=xdn;
            end 
            f2=fm+k1*(x2(1)-xm(1));
        else
        end
    elseif(x1(2)<0)
        if(x2(2)>0)%采用精确拐点非迭代法求出第二类零速度拐点；
            kk2=k1;
            [tm,xm]=guaidian2(dt,k2,m,eta,xg,x1);     
%             fm=f1+k2*(xm(1)-x1(1));
            xdn=xm(1);
            xup=xdn+2*uy;
            fup=f0+k2*xup;
            fdn=fup-2*py;
            fm=fdn;
            txg(1)=xg(1)+(xg(2)-xg(1))*tm/dt;
            txg(2)=xg(2);
            dt1=dt-tm;
            [x2]=newmark_single(dt1,k1,m,eta,txg,xm);  
            if(x2(1)>xup)%此处近似处理是否恰当？
                x2(1)=xup;
            elseif(x2(1)<xdn)
                x2(1)=xdn;
            end  
            f2=fm+k1*(x2(1)-xm(1));
        else
        end
    else
    end
else
end