v2_007.m

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

% v2_007.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 using numerical integration.
%            
% The method uses the ODE45 function in MATLAB.
%
% Consider the second order differential equation for a multiple
% degree of freedom vibrating system
%
% [M] * { x''} +  [C] * { x' } + [K] * { x } = { f }
%
% where the prime ' denotes the derivative operator.  
% We can rewrite this equation as follows:
%
% { x''} = inv([M]) * ({f} - [C]*{x'} - [K]*{x})
 
%**********************************************************************
% Author: Randall J. Allemang
% Date:	22-FEB-2000
% Structural Dynamics Research Lab
% University of Cincinnati
% Cincinnati, Ohio  45221-0072
% TEL:  [phone]
% FAX:  [phone]
% E-MAIL: [e-mail]
%*********************************************************************
%
clear
% close all
plt=input('Store plots to file (Yes=1): (0)');if isempty(plt),plt=0;end;
pi=3.14159265;
NDOF=2;
M=[10,0,0;0,14,0;0,0,12];
K=1000*[5,-3,0;-3,5.5,-2.5;0,-2.5,2.5];    
C=0.001*K;
Z=[0,0,0;0,0,0;0,0,0];
%    Form 2N x 2N state space equation.  
A=[Z,M;M,C];
B=[-M,Z;Z,K];
[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
freq=imag(lambda)/(2*pi);
period=1.0./(freq');
lambda,freq,period
clear d x
xxx=input('Hit any key to continue');
%********************************************************************
%	Define time limits and increment
tinit = 0;
dt=input('Time Increment (delta t): (.002)');if isempty(dt),dt=0.002;end;
Nt=input('Number of Time Points (Nt): (8001)');if isempty(Nt),Nt=8001;end;
tfinal = (Nt-1)*dt;
t=linspace(tinit,tfinal,Nt);
%	Define initial conditions

x1_init = 0;
v1_init = 0;
x2_init = 0;
v2_init = 0;
x3_init = 0;
v3_init = 0;

%  Define force function
f1(Nt)=0;
f2(Nt)=0;
f3(Nt)=0;
%	Define steady state excitation (cosine)
f1=100*cos(18.*t);
f=[f1;f2;f3];
[t,x]=ode45('mkc',[0 20],[0 0 0 0 0 0]);

%	Plot results
figure
subplot(211),plot(t,f(1,:)),grid,title('Force 1')
subplot(212),plot(t,x(1,:)),grid,title('Displacement 1')
if plt==1, print -deps v2_007a, end;

figure
subplot(211),plot(t,f(2,:)),grid,title('Force 2')
subplot(212),plot(t,x(2,:)),grid,title('Displacement 2')
if plt==1, print -deps v2_007b, end;

figure
subplot(211),plot(t,f(3,:)),grid,title('Force 3')
subplot(212),plot(t,x(3,:)),grid,title('Displacement 3')
if plt==1, print -deps v2_006c, end;

function dx=mkc(t,x)

M=[10,0,0;0,14,0;0,0,12];
K=1000*[5,-3,0;-3,5.5,-2.5;0,-2.5,2.5];    
C=0.001*K;