v2_006.m

来自「美国辛辛那提大学的振动教材」· M 代码 · 共 171 行

% v2_006.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 either the Euler approach or a Runge Kutta approach
% to solving the ODE.
% To achieve reasonable accuracy, a small time step size
% is required.  Accuracy for the Euler approach is limited to (delta t ^^2).
% Accuracy for the Runge Kutta approach is limited to (delta t ^^4).
% While the Euler method is simple to understand and explain, the method
% is not stable and may diverge.
%
% The Runge Kutta method is developed without the use 
% of any matlab functions (ODE23 or ODE45) based upon the presentation
% in Thomson's vibration book.
%            
% Reference:  Theory of Vibration with Applications
%             4th Edition, 1993
%             William T. Thomson
%
% For better accuracy, utilize the ODE scripts provided 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.1569*2;
% v1_init = 0.2081*2*18;
% x2_init = -0.0795*2;
% v2_init = 0.1188*2*18;
% x3_init = 0.1537*2;
% v3_init = -0.2064*2*18;

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];
method=input('Choose numerical integration method (Euler=0, Runge Kutta=1): (0)');
if isempty(method),method=0;end;

% Start Euler Method
if method == 0
	x(1,1)=x1_init;
   x(2,1)=x2_init;
   x(3,1)=x3_init;
	v(1,1)=v1_init;
   v(2,1)=v2_init;
   v(3,1)=v3_init;
	a(:,1)=inv(M)*(f(:,1)-C*v(:,1)-K*x(:,1));
   for ii=2:Nt
   x(:,ii)=x(:,ii-1)+v(:,ii-1).*dt;
	v(:,ii)=v(:,ii-1)+a(:,ii-1).*dt;
	a(:,ii)=inv(M)*(f(:,ii)-C*v(:,ii)-K*x(:,ii));
	end

%	Start Runge Kutta Method 
elseif method == 1
   x(1,1)=x1_init;
   x(2,1)=x2_init;
   x(3,1)=x3_init;
	v(1,1)=v1_init;
   v(2,1)=v2_init;
   v(3,1)=v3_init;
   a(:,1)=inv(M)*(f(:,1)-C*v(:,1)-K*x(:,1));
	for ii=2:Nt
	T1=t(ii-1);
	X1=x(:,ii-1);
	V1=v(:,ii-1);
	A1=inv(M)*(f(:,ii-1)-C*V1-K*X1);
	T2=t(ii-1)+dt/2;
	X2=x(:,ii-1)+V1*dt/2;
	V2=v(:,ii-1)+A1*dt/2;
	A2=inv(M)*(((f(:,ii-1)+f(:,ii))./2)-K*X2-C*V2);
	T3=t(ii-1)+dt/2;
	X3=x(:,ii-1)+V2*dt/2;
	V3=v(:,ii-1)+A2*dt/2;
	A3=inv(M)*(((f(:,ii-1)+f(:,ii))./2)-K*X3-C*V3);
	T4=t(ii);
	X4=x(:,ii-1)+V3*dt;
	V4=v(:,ii-1)+A3*dt;
	A4=inv(M)*(f(:,ii)-K*X3-C*V3);
	x(:,ii)=x(:,ii-1) + (V1+2*V2+2*V3+V4)*dt/6;
	v(:,ii)=v(:,ii-1) + (A1+2*A2+2*A3+A4)*dt/6;
	a(:,ii)=inv(M)*(f(:,ii)-K*x(:,ii)-C*v(:,ii));
	end
end   

%	Plot results
figure
subplot(221),plot(t,f(1,:)),grid,title('Force 1')
subplot(222),plot(t,x(1,:)),grid,title('Displacement 1')
subplot(223),plot(t,v(1,:)),grid,title('Velocity 1')
subplot(224),plot(t,a(1,:)),grid,title('Acceleration 1')
if plt==1, print -deps v2_006a, end;
figure
subplot(221),plot(t,f(2,:)),grid,title('Force 2')
subplot(222),plot(t,x(2,:)),grid,title('Displacement 2')
subplot(223),plot(t,v(2,:)),grid,title('Velocity 2')
subplot(224),plot(t,a(2,:)),grid,title('Acceleration 2')
if plt==1, print -deps v2_006b, end;
figure
subplot(221),plot(t,f(3,:)),grid,title('Force 3')
subplot(222),plot(t,x(3,:)),grid,title('Displacement 3')
subplot(223),plot(t,v(3,:)),grid,title('Velocity 3')
subplot(224),plot(t,a(3,:)),grid,title('Acceleration 3')
if plt==1, print -deps v2_006c, end;