Can anybody help me to implement caughey's equivalent linearization method for bilinear structure in MATLAB? I am not able to get the ductility(x/x0) vs frequency ratio(w/ws) plot.

These are some images from "Passive control of bilinear hysteretic structures by tuned mass damper

for narrow band seismic motions" Z. Zhang, T. Balendra.

[darova]

Do you have any attempts?

[VAIBHAV PUJARI]

clc
clear all

result=zeros(1,300);
x_ax=(0:0.01:3);
aa=0; 
ws=sqrt(2.85e9/1.8e7); ms=1.8e7; k0=2.85e9; c0=0.015*2*sqrt(ms*k0); n=0.3; nu=0.02; x00=0.02; A=0.1*9.81;

ws=sqrt(2.85e9/1.8e7);

for w=x_ax*ws
% for w=ws
    aa=aa+1;
    x0=[1 1];

    [sol, fval]=fminunc(@(x)response(x,w, ms,k0,c0,n,nu,x00,A),x0);
    result(1,aa)=fval;
    ductility=result/x00;  
end

plot(x_ax,ductility,'r-*')


function DMF=response(x,w, ms,k0,c0,n,nu,x00,A)

wd=x(1,1); zid=x(1,2);
% ws=sqrt(2.85e9/1.8e7); ms=1.8e7; k0=2.85e9; c0=0.015*2*sqrt(ms*k0); n=0.3; nu=0.02; x00=0.02; A=0.1*9.81;
% ms=1.8e7; k0=2.85e6; c0=0.015*2*sqrt(ms*k0); n=0.3; nu=0.02; x00=0.02e-3; A=0.1;
md=nu*ms; kd=md*wd^2; cd=zid*2*md*wd;



    ke=k0; ce=c0; we=sqrt(ke/ms); zie=(ce/(2*ms*we));

    Z=md*(wd^2+2*1i*zid*wd*w)/(wd^2-w^2+2*1i*zid*wd*w);
    H=(ms+Z)/(ms*(we^2-w^2++2*1i*zie*we*w)-w^2*Z);     % should be multiplied by A*w^2

    RE=real(1/H);
    IM=imag(1/H);
    DMF=(1/sqrt(RE^2+IM^2))*ms*A/ke;

    tita=acos(1-2*x00/(DMF));    %%% Mistake: for inelastic, value of elastic DMF is being used 
                                            % see equation of ke, tita is
                                            % used
    if DMF<=x00
         disp('elastic')
        DMF
        
    else
        disp('inelastic')
        ke=(k0/pi)*((1-n)*tita+n*pi-0.5*(1-n)*sin(2*tita));
        we=sqrt(ke/ms);
        ce=c0+((k0/(pi*we))*(sin(tita))^2);
        zie=(ce/(2*ms*we));
        Z=md*(wd^2+2*1i*zid*wd*w)/(wd^2-w^2+2*1i*zid*wd*w);
        H=(ms+Z)/(ms*(we^2-w^2+2*1i*zie*we*w)-w^2*Z);   % should be multiplied by A*w^2

        RE=real(1/H); IM=imag(1/H);
        DMF=(1/sqrt(RE^2+IM^2))*ms*A/ke;
    end

end

[VAIBHAV PUJARI]

i have written this code. The values of various parameters are given in that paper, equations of H and Z are derived by me which match with the equations given in paper. please help..

I am stucked at this point from last many days and not able to understand the mistake

[darova]

I don't understand. HOw to check it? Where are original formulas?

[VAIBHAV PUJARI]

Please check. I have attached all required formulae.

[darova]

It's too complicated for me. Sorry, can't help

[VAIBHAV PUJARI]

ok.. Thanks.