v2_004.m

美国辛辛那提大学的振动教材

% v2_004.m   
%
% This is a script file to solve a mdof system,
% given the mass, damping and stiffness terms
% in dimensionless units,  for the time response of the 
% system.
%            
% The method uses the Euler approach to solving the ODE.
% To achieve reasonable accuracy, a small time step size
% is required.  Accuracy is limited to (delta t ^^2).
%
% Consider the second order differential equation for a multiple
% degree of freedom vibrating system
%
% [A] * { y' } + [B] * { y } = { f }
%
% where the prime ' denotes the derivative operator.  
% We can rewrite this equation as follows:
%
% { y' } = inv([A]) * [ { f } - [B] * { y } ]
 
%********************************************************************
% Author: Randall J. Allemang
% Date:	18-Apr-94
% Structural Dynamics Research Lab
% University of Cincinnati
% Cincinnati, Ohio  45221-0072
% TEL:  [phone]
% FAX:  [phone]
% E-MAIL: [e-mail]
%*********************************************************************
%
%  Set up State-Space Problem
%
clear,clg
plt=input('Store plots to file (Yes=1): (0)');if isempty(plt),plt=0;end;
pi=3.14159265;
NDOF=2;
mass=[5,0;0,10];
stiff=500*[4,-2;-2,6];    
damp=0.01*stiff;
null=[0,0;0,0];
%    Form 2N x 2N state space equation.  
A=[null,mass;mass,damp];
B=[-mass,null;null,stiff];
[x,d]=eig(-inv(A)*B);
%
% Sort Modal Frequencies  
%
orig_lambda=diag(d);
[Y,I]=sort(imag(orig_lambda));
lambda=orig_lambda(I);
psi=x(:,I);
%     Modal Frequencies
clear d
xxx=input('Hit any key to continue');
%********************************************************************
%	Define time limits and increment
tinit = 0;
tfinal = 2;
t=linspace(tinit,tfinal,1001);
dt=t(2)-t(1)
%	Define initial conditions
x1_init = 1;
v1_init = 0;
x2_init = 1;
v2_init = 0;
%       Clear force function
for ii=1:1001
f1(ii)=0;
f2(ii)=0;
Yf(1,ii)=0;
Yf(2,ii)=0;
Yf(3,ii)=0;
Yf(4,ii)=0;
end
%	Define transient excitation (sine pulse)
time_dur=0.4;
numb_int=round(time_dur/dt)
for ii=1:numb_int
f1(ii)=100*sin(2.*pi./(2.*time_dur).*t(ii));
end
plot(t,f1), title('Force time history'), pause
clg
%	Start Euler ODE approximation
x(1,1)=x1_init;
x(2,1)=x2_init;
v(1,1)=v1_init;
v(2,1)=v2_init;
Yx=[v(:,1);x(:,1)];
Yf(3,:)=f1;
Yv=inv(A)*(Yf(:,1)-(B*Yx));
a(1,1)=Yv(1);
a(2,1)=Yv(2);
for ii=2:1001
x(:,ii)=x(:,ii-1)+v(:,ii-1).*dt;
v(:,ii)=v(:,ii-1)+a(:,ii-1).*dt;
clear Yx
Yx=[v(:,ii);x(:,ii)];
Yv=inv(A)*(Yf(:,ii)-(B*Yx));
a(1,ii)=Yv(1);
a(2,ii)=Yv(2);
clear Yv
end
%	Plot results
subplot(211),plot(t,f1),grid
subplot(212),plot(t,f2),grid
title('Force time history'), pause
if plt==1, meta v2_004a, end;
clg
subplot(211),plot(t,x(1,:)),grid
subplot(212),plot(t,x(2,:)),grid
title('2DOF time history (displacement)'), pause
if plt==1, meta v2_004b, end;
clg
subplot(211),plot(t,v(1,:)),grid
subplot(212),plot(t,v(2,:)),grid
title('2DOF time history (velocity)'), pause
if plt==1, meta v2_004c, end;
clg
subplot(211),plot(t,a(1,:)),grid
subplot(212),plot(t,a(2,:)),grid
title('2DOF time history (acceleration)'), pause
if plt==1, meta v2_004d, end;