Document Actions
v1_005.m
% v1_005.m
%
% This a script file to compare the euler solution
% with the exact solution for a transient excitation
% case. Solution is for a single DOF system with a
% sine pulse excitation.
%**********************************************************************
% Author: Randall J. Allemang
% Date: 25-Aug-97
% Structural Dynamics Research Lab
% University of Cincinnati
% Cincinnati, Ohio 45221-0072
% TEL: 513-556-2725
% FAX: 513-556-3390
% E-MAIL: randy.allemang@uc.edu
%*********************************************************************
%
clg, clear
plt=input('Store plots to file (Yes=1): (0)');
if isempty(plt),plt=0;end;
mass=input('Mass value: (20)');
if isempty(mass), mass=20;end;
stiff=input('Stiffness value: (9000)');
if isempty(stiff),stiff=9000;end;
damp=input('Damping value: (130)');
if isempty(damp),damp=130;end;
%
lambda=roots([mass,damp,stiff]);
sigma=real(lambda(1));
om1=abs(imag(lambda(1)));
lambda,sigma,om1
pause
om2=1.0/0.8*(2.0*pi);
om3=om2-om1;
om4=om2+om1;
om1,om2,om3,om4
pause
% Define time limits and increment
tinit = 0;
tfinal = 2;
t=linspace(tinit,tfinal,501);
dt=t(2)-t(1)
% Define initial conditions
xinit = 0;
vinit = 0;
% Clear force function
for ii=1:501
f(ii)=0;
end
% Define transient excitation (sine pulse)
time_dur=0.1;
numb_int=round(time_dur/dt)
for ii=1:numb_int
f(ii)=15*sin(2.*pi./(2.*time_dur).*t(ii));
end
plot(t,f), title('SDOF time history (force)'), pause
clg
% Start Computation
for ii=1:501
tt=t(ii);
if tt <= time_dur
A1=15.*exp(sigma.*tt).*sin(om1.*tt)./(om1.*mass);
B1=A1*(om3+exp(-sigma*tt)*(-sigma.*sin(om3*tt)-om3*cos(om3*tt)));
C1=B1/(2*sigma*sigma + 2*om3*om3);
D1=A1*(om4+exp(-sigma*tt)*(-sigma.*sin(om4*tt)-om4*cos(om4*tt)));
E1=D1/(2*sigma*sigma + 2*om4*om4);
x1(ii)=C1+E1;
A2=-15.*exp(sigma.*tt).*cos(om1.*tt)./(om1.*mass);
B2=A2*(sigma+exp(-sigma*tt)*(om3.*sin(om3*tt)-sigma*cos(om3*tt)));
C2=B2/(2*sigma*sigma + 2*om3*om3);
D2=A2*(-sigma-exp(-sigma*tt)*(om4.*sin(om4*tt)-sigma*cos(om4*tt)));
E2=D2/(2*sigma*sigma + 2*om4*om4);
x2(ii)=C2+E2;
x(ii)=x1(ii)+x2(ii);
else
A1=15.*exp(sigma.*tt).*sin(om1.*tt)./(om1.*mass);
B1=A1*(om3+exp(-sigma*0.1)*(-sigma.*sin(om3*0.1)-om3*cos(om3*0.1)));
C1=B1/(2*sigma*sigma + 2*om3*om3);
D1=A1*(om4+exp(-sigma*0.1)*(-sigma.*sin(om4*0.1)-om4*cos(om4*0.1)));
E1=D1/(2*sigma*sigma + 2*om4*om4);
x1(ii)=C1+E1;
A2=-15.*exp(sigma.*tt).*cos(om1.*tt)./(om1.*mass);
B2=A2*(sigma+exp(-sigma*0.1)*(om3.*sin(om3*0.1)-sigma*cos(om3*0.1)));
C2=B2/(2*sigma*sigma + 2*om3*om3);
D2=A2*(-sigma-exp(-sigma*0.1)*(om4.*sin(om4*0.1)-sigma*cos(om4*0.1)));
E2=D2/(2*sigma*sigma + 2*om4*om4);
x2(ii)=C2+E2;
x(ii)=x1(ii)+x2(ii);
end
end
% Plot results
fig1=figure(1);
plot(t,f), title('SDOF time history (force)'), grid, pause
fig2=figure(2)
if plt==1, print -ddeskjet -f1 v1_005a, end;
plot(t,x), title('Exact Solution - SDOF time history (displacement)')
grid, pause
if plt==1, print -ddeskjet -f2 v1_005b, end;