Dynamic Response of SDOF - Matlab Code

구조동역학

▣ Problems SDOF system : m=44.357kN・s2/m, k=1751.18kN/m, ξ=0.05 ω=2π rad/s, T=1.0s

Input loading; p(t)=44.48 sin(πt/0.6)kN Initial Condition; v(0)=v’(0)=p(0)=0

(1) Determine the theoretical solution sol) general solution of damped system subjected to harmonic loading is                  cos          

  si n    

  

               si n         cos                 



          

      cos       si n           

              si n         cos                 



accordance with initial condition      

                         



        si n           cos               cos       si n         

                  

 

     

              cos           si n                     



                                



                               



                                                                                                         







∴    



 

 



 

                  cos                                   





si n   

      

              si n        cos                 



1

구조동역학

▪ plot result

(2) Calculate the response by using the central difference method (dt=0.1s) ▪ matlab code %---------------------------------------% Central Diffenence Method %----------------------------------------

% input data m=44.357;

% mass (unit=kN.s2/m)

k=1751.18;

% stiffness (unit=kN/m)

h=0.05;

% dampin ratio (h=c/ccr=c/2mw)

w=2*pi;

% natural angular frequency of structure (unit=rad/s)

T=1.0;

% period (unit=sec)

c=2*h*m*w;

% damping coefficient

% define input loading dt=0.1;

% for stability, dt <= 0.318T

t=[0:dt:10]';

% time

P=44.48*sin(pi*t/0.6); n=size(P,1);

% input loading (unit=kN) % lengh of load vector(n=101)

% define initial values F(1)=0; F(2)=0; d(2)=0; v(2)=0; a(2)=1/m*(P(2)-c*v(2)-k*d(2)); d(1)=d(2)-dt*v(2)+dt^2/2*a(2);

% actually i=0 i.e d0,v0,a0

% actually i=-1

2

구조동역학

% coefficient of d(i+1) A=m/dt^2-c/(2*dt); B=k-2*m/dt^2; C=m/dt^2+c/(2*dt);

% iteration for i=3:n+2 TT(i)=t(i-2); F(i)=P(i-2); d(i)=2*dt*v(i-1)+d(i-2); d(i+1)=1/C*(F(i)-A*d(i-1)-B*d(i)); v(i)=(d(i+1)-d(i-1))/(2*dt); a(i)=(d(i+1)-2*d(i)+d(i-1))/dt^2; v(i+1)=v(i);

% virtual

a(i+1)=a(i);

% virtual

TT(i+1)=TT(i);

% virtual

F(i+1)=F(i);

% virtual

end

result=[TT' F' d' v' a']; xlswrite( 'cdm result.xls',result) result.xls' ,result)

figure(1) plot(t,P) xlabel('time(sec)' xlabel( 'time(sec)' ) ylabel('Load(kN)' ylabel( 'Load(kN)' ) title('Input title( 'Input loading' ) grid on

figure(2) plot(TT,d) xlabel('time(sec)' xlabel( 'time(sec)' ) ylabel('displacement(m)' ylabel( 'displacement(m)' ) title('Displacement title( 'Displacement Response' ) grid on

figure(3) plot(TT,v) xlabel('time(sec)' xlabel( 'time(sec)' ) ylabel('velocity(m/s)' ylabel( 'velocity(m/s)' ) title('Velocity title( 'Velocity Response' ) grid on

figure(4) plot(TT,a) xlabel('time(sec)' xlabel( 'time(sec)' ) ylabel('acceleration(m/s2)' ylabel( 'acceleration(m/s2)' ) title('Acceleration title( 'Acceleration Response' ) grid on

3

구조동역학

(3) Calculate the response by using the constant average acceleration method (dt=0.1s) ▪ matlab code

%---------------------------------------% Constant Average Acceleraion Method %----------------------------------------

% input data m=44.357;

% mass (unit=kN.s2/m)

k=1751.18;

% stiffness (unit=kN/m)

h=0.05;

% dampin ratio (h=c/ccr= (h=c/ccr=c/2mw) c/2mw)

w=2*pi;

% natural angular frequency of structure (unit=rad/s)

T=1.0;

% period (unit=sec)

c=2*h*m*w;

% damping coefficient

% define input loading dt=0.1;

% for stability, dt <= 0.318T

t=[0:dt:10]'; P=44.48*sin(pi*t/0.6); n=size(P,1);

% time % input loading (unit=kN) % lengh of load vector(n=101)

% define initial values (i=0 i.e d0,v0,a0) d(1)=0; v(1)=0; a(1)=1/m*(P(1)-c*v(1)-k*d(1));

% iteration for i=1:n-1

d(i+1)=1/(k+2*c/dt+4*m/dt^2)*(P(i+1)+m*(4/dt^2*d(i)+4/dt*v(i)+a(i))+c*(2/dt*d(i)+ v(i))); v(i+1)=2/dt*(d(i+1)-d(i))-v(i); a(i+1)=4/dt^2*(d(i+1)-d(i))-4/dt*v(i)-a(i); end

4

구조동역학

(4) Calculate the response by using the linear acceleration method (dt=0.1s) ▪ matlab code %---------------------------------------% Linear Acceleraion Method %----------------------------------------

% input data m=44.357;

% mass (unit=kN.s2/m)

k=1751.18;

% stiffness (unit=kN/m)

h=0.05;

% dampin ratio (h=c/ccr= (h=c/ccr=c/2mw) c/2mw)

w=2*pi;

% natural angular frequency of structure (unit=rad/s)

T=1.0;

% period (unit=sec)

c=2*h*m*w;

% damping coefficient

% define input loading dt=0.1;

% for stability, dt <= 0.318T

t=[0:dt:10]'; P=44.48*sin(pi*t/0.6); n=size(P,1);

% time % input loading (unit=kN) % lengh of load vector(n=101)

% define initial values (i=0 i.e d0,v0,a0) d(1)=0; v(1)=0; a(1)=1/m*(P(1)-c*v(1)-k*d(1));

% iteration for i=1:n-1

d(i+1)=1/(k+6*m/dt^2+3*c/dt)*(P(i+1)+(6*m/dt^2+3*c/dt)*d(i)+(6*m/dt+2*c)*v(i)+(2* m+c*dt/2)*a(i)); v(i+1)=3/dt*(d(i+1)-d(i))-2*v(i)-dt/2*a(i); a(i+1)=6/dt^2*(d(i+1)-d(i))-6/dt*v(i)-2*a(i); end

5

구조동역학

(5) Plot Result

6