cantbeam_ss_modred.m

vibration simulation using ansys and matlab 一书中的程序

	echo off
%	cantbeam_ss_modred.m  Creates a Matlab state space model 
%	using the eigenvalue and eigenvector results from previous 
%	ANSYS runs.  Modes are ranked for importance and several reduction 
%	techniques are used.

   	clear all;
	
	hold off;

	clf;
				
%	load the .mat file cantbeamXXred, containing evr - the modal matrix, freqvec -
%	the frequency vector and node_numbers - the vector of node numbers for the modal
%	matrix

	model = menu('choose which finite element model to use ... ', ....
				'2 beam elements', ...
				'4 beam elements', ...
				'6 beam elements', ...
				'8 beam elements', ...
				'10 beam elements', ...
				'12 beam elements', ...
				'16 beam elements', ...
				'32 beam elements', ...
				'64 beam elements');

	if  model == 1
		load cantbeam2red;		
	elseif  model == 2
		load cantbeam4red;		
	elseif  model == 3
		load cantbeam6red;		
	elseif  model == 4
		load cantbeam8red;		
	elseif  model == 5
		load cantbeam10red;		
	elseif  model == 6
		load cantbeam12red;		
	elseif  model == 7
		load cantbeam16red;		
	elseif  model == 8
		load cantbeam32red;		
	elseif  model == 9
		load cantbeam64red;		
	end

%	define the number of degrees of freedom and number of modes from size of modal matrix

	[numdof,num_modes_total] = size(evr);

%	define rows for middle and  tip nodes

	mid_node_row = 0.5*(numdof-1)+1;

	tip_node_row = numdof;

	xn = evr;
 
%   calculate the dc amplitude of the displacement of each mode by 
%   multiplying the forcing function row of the eigenvector by the output row

   	omega2 = (2*pi*freqvec)'.^2;     %  convert to radians and square

	dc_gain = abs(xn(mid_node_row,:).*xn(tip_node_row,:))./omega2;

	[dc_gain_sort,index_sort] = sort(dc_gain);

	dc_gain_sort = fliplr(dc_gain_sort);

	index_sort = fliplr(index_sort)

	dc_gain_nosort = dc_gain;

    index_orig = 1:num_modes_total;
 
   	semilogy(index_orig,freqvec,'ko-');
   	title('frequency versus mode number')
   	xlabel('mode number')
   	ylabel('frequency, hz')
   	grid on
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,dc_gain_nosort,'ko-')
   	title('dc value of each mode contribution versus mode number')
   	xlabel('mode number')
   	ylabel('dc value')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	loglog(freqvec,dc_gain_nosort,'ko-')
   	title('dc value of each mode contribution versus frequency')
   	xlabel('frequency, hz')
   	ylabel('dc value')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,dc_gain_sort,'ko-')
   	title('sorted dc value of each mode versus number of modes included')
   	xlabel('modes included')
   	ylabel('sorted dc value')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	num_modes_used = input(['enter how many modes to include, ',num2str(num_modes_total),' default, max ...  ']);
		
    if (isempty(num_modes_used))
        num_modes_used = num_modes_total;
    end

   	zeta = input('enter value for damping, .02 is 2% of critical (default) ... ');
	
    if (isempty(zeta))
    	zeta = .02;
    end

%	all modes included model, use original order

	xnnew = xn(:,(1:num_modes_total));

	freqnew = freqvec((1:num_modes_total));

%	reduced, no sorting, just use the first num_modes_used modes in xnnew_nosort

	xnnew_nosort = xn(:,1:num_modes_used);

	freqnew_nosort = freqvec(1:num_modes_used);

%	reduced, sorting, use the first num_modes_used sorted modes in xnnew_sort

	xnnew_sort = xn(:,index_sort(1:num_modes_used));

	freqnew_sort = freqvec(index_sort(1:num_modes_used));

%	define variables for all modes included system matrix, a

    w = freqnew*2*pi;		%  	frequencies in rad/sec

	w2 = w.^2;				

    zw = 2*zeta*w;

%	define variables for reduced, nosorted system matrix, a_nosort

    w_nosort = freqnew_nosort*2*pi;		%  	frequencies in rad/sec

	w2_nosort = w_nosort.^2;				

    zw_nosort = 2*zeta*w_nosort;

%	define variables for reduced, sorted system matrix, a_sort

    w_sort = freqnew_sort*2*pi;		%  	frequencies in rad/sec

	w2_sort = w_sort.^2;				

    zw_sort = 2*zeta*w_sort;

%	define size of system matrix

	asize = 2*num_modes_total;

   	asize_red = 2*num_modes_used;

   	disp('   ');
   	disp('   ');
   	disp(['size of system matrix a is ',num2str(asize)]);
   	disp(['size of reduced system matrix a is ',num2str(asize_red)]);

%  	setup all modes included "a" matrix, system matrix

   	a = zeros(asize);

   	for  col = 2:2:asize

      	row = col-1;

      	a(row,col) = 1;

   	end

   	for  col = 1:2:asize

      	row = col+1;

      	a(row,col) = -w2((col+1)/2);

   	end

   	for  col = 2:2:asize

      	row = col;

      	a(row,col) = -zw(col/2);

   	end

%  	setup reduced, nosorted "a_nosort" matrix, system matrix

   	a_nosort = zeros(asize_red);

   	for  col = 2:2:asize_red

      	row = col-1;

      	a_nosort(row,col) = 1;

   	end

   	for  col = 1:2:asize_red

      	row = col+1;

      	a_nosort(row,col) = -w2_nosort((col+1)/2);

   	end

   	for  col = 2:2:asize_red

      	row = col;

      	a_nosort(row,col) = -zw_nosort(col/2);

   	end

%  	setup reduced, sorted "a_sort" matrix, system matrix

   	a_sort = zeros(asize_red);

   	for  col = 2:2:asize_red

      	row = col-1;

      	a_sort(row,col) = 1;

   	end

   	for  col = 1:2:asize_red

      	row = col+1;

      	a_sort(row,col) = -w2_sort((col+1)/2);

   	end

   	for  col = 2:2:asize_red

      	row = col;

      	a_sort(row,col) = -zw_sort(col/2);

   	end

%	setup input matrix b, state space forcing function in principal coordinates

%  	f_physical is the vector of physical force
%  	zeros at each output DOF and input force at the input DOF

	f_physical = zeros(numdof,1);			%  	start out with zeros

	f_physical(mid_node_row) = 1.0;			%   input force at node 6, midpoint node

%	f_principal is the vector of forces in principal coordinates

   	f_principal = xnnew'*f_physical;

%	b is the vector of forces in principal coordinates, state space form

	b = zeros(2*num_modes_total,1);

	for  cnt = 1:num_modes_total
	
		b(2*cnt) = f_principal(cnt);
		
	end				

%	f_principal_nosort is the vector of forces in principal coordinates

   	f_principal_nosort = xnnew_nosort'*f_physical;

%	b_nosort is the vector of forces in principal coordinates, state space form

	b_nosort = zeros(2*num_modes_used,1);

	for  cnt = 1:num_modes_used
	
		b_nosort(2*cnt) = f_principal_nosort(cnt);
		
	end				

%	f_principal_sort is the vector of forces in principal coordinates

   	f_principal_sort = xnnew_sort'*f_physical;

%	b_sort is the vector of forces in principal coordinates, state space form

	b_sort = zeros(2*num_modes_used,1);

	for  cnt = 1:num_modes_used
	
		b_sort(2*cnt) = f_principal_sort(cnt);
		
	end				

%  	setup cdisp and cvel, padded xn matrices to give the displacement and velocity
%	vectors in physical coordinates
%	cdisp and cvel each have numdof rows and alternating columns consisting of columns
%   of xnnew and zeros to give total columns equal to the number of states

%	all modes included cdisp and cvel

   	for  col = 1:2:2*length(freqnew)

      	for  row = 1:numdof

         	cdisp(row,col) = xnnew(row,ceil(col/2));

         	cvel(row,col) = 0;
   
      	end

   	end

   	for  col = 2:2:2*length(freqnew)

      	for  row = 1:numdof

         	cdisp(row,col) = 0;

	        cvel(row,col) = xnnew(row,col/2);

      	end

   	end

%	reduced, nosorted cdisp and cvel

   	for  col = 1:2:2*length(freqnew_nosort)

      	for  row = 1:numdof

         	cdisp_nosort(row,col) = xnnew_nosort(row,ceil(col/2));

         	cvel_nosort(row,col) = 0;
   
      	end

   	end

   	for  col = 2:2:2*length(freqnew_nosort)

      	for  row = 1:numdof

         	cdisp_nosort(row,col) = 0;

	        cvel_nosort(row,col) = xnnew_nosort(row,col/2);

      	end

   	end

%	reduced, sorted cdisp and cvel

   	for  col = 1:2:2*length(freqnew_sort)

      	for  row = 1:numdof

         	cdisp_sort(row,col) = xnnew_sort(row,ceil(col/2));

         	cvel_sort(row,col) = 0;
   
      	end

   	end

   	for  col = 2:2:2*length(freqnew_sort)

      	for  row = 1:numdof

         	cdisp_sort(row,col) = 0;

	        cvel_sort(row,col) = xnnew_sort(row,col/2);

      	end

   	end

%  	define output

   	d = [0];  %

%	define frequency vector for frequency responses

   	freqlo = 10;

   	freqhi = 100000;

   	flo=log10(freqlo) ;
   	fhi=log10(freqhi) ;
 
   	f=logspace(flo,fhi,200) ;
   	frad=f*2*pi ;
 
%  	take transfer functions, outputting the midpoint and tip node rows of the displacement
%	vector cdisp

%	define displacement state space system with the "ss" command

	sysdisptip = ss(a,b,cdisp(tip_node_row,:),d);

%	defined reduced systems using num_modes_used nosort modes

	sysdisptip_nosort = ss(a_nosort,b_nosort,cdisp_nosort(tip_node_row,:),d);

%	define reduced systems using num_modes_used sorted modes

	sysdisptip_sort = ss(a_sort,b_sort,cdisp_sort(tip_node_row,:),d);

%	use "bode" command to generate magnitude/phase vectors

   	[magdisptip,phsdisptip]=bode(sysdisptip,frad)  ;
 
   	[magdisptip_nosort,phsdisptip_nosort]=bode(sysdisptip_nosort,frad)  ;
 
   	[magdisptip_sort,phsdisptip_sort]=bode(sysdisptip_sort,frad)  ;
 
%	convert magnitude to db

   	magdisptipdb = 20*log10(magdisptip);

   	magdisptipdb_nosort = 20*log10(magdisptip_nosort);

   	magdisptipdb_sort = 20*log10(magdisptip_sort);

%	start plotting

	if  num_modes_used == num_modes_total

%	plot all modes included response

	semilogx(f,magdisptipdb(1,:),'k.-')
	title(['cantilever tip displacement for mid-length force, all ',num2str(num_modes_used),' modes included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, db mm')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	hold on

	max_modes_plot = num_modes_total;