cantbeam_ss_shkr_modred.m

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

	echo off
%	cantbeam_ss_shkr_modred.m  Creates a Matlab state space model 
%	using the results from ANSYS model cantbeam_ss_spring_shkr.inp.  
%	Ranks modes and then uses several reduction techniques to define 
%	smaller model.

   	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_shkr;		
	elseif  model == 2
		load cantbeam4red_shkr;		
	elseif  model == 3
		load cantbeam6red_shkr;		
	elseif  model == 4
		load cantbeam8red_shkr;		
	elseif  model == 5
		load cantbeam10red_shkr;		
	elseif  model == 6
		load cantbeam12red_shkr;		
	elseif  model == 7
		load cantbeam16red_shkr;		
	elseif  model == 8
		load cantbeam32red_shkr;		
	elseif  model == 9
		load cantbeam64red_shkr;		
	end

	kspring = 1000000;			% mN/mm from ANSYS run

	shaker_mass = 0.050;		% kg from ANSYS run

	mn2gm_conversion = 0.101968;	%  conversion factor from mn to gram-f, 1/9.807

%	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 shaker and tip nodes

	shaker_node_row = 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

%	define frequency range for frequency response

   	freqlo = 100;

   	freqhi = 100000;

   	flo=log10(freqlo) ;
   	fhi=log10(freqhi) ;
 
   	f=logspace(flo,fhi,200) ;
   	frad=f*2*pi ;
 
	dc_gain = abs([xn(shaker_node_row,1)*xn(tip_node_row,1)/(frad(1)^2)  ...
			  	   (xn(shaker_node_row,2:num_modes_total) ...
			  	  .*xn(tip_node_row,2:num_modes_total))./omega2(2:num_modes_total)]);

	[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
	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([freqlo; freqvec(2:num_modes_total)],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(shaker_node_row) = 9807*shaker_mass*1.0;		%  input force at shaker, 1g

%	now setup the principal force vector for the three cases, all modes, nosort, sort

%	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 tip force state space system with the "ss" command

	sysforce = ss(a,b,mn2gm_conversion*kspring*(cdisp(tip_node_row,:)-cdisp(shaker_node_row,:)),d);

%	define reduced system using nosort modes

	sysforce_nosort = ss(a_nosort,b_nosort,mn2gm_conversion*kspring*(cdisp_nosort(tip_node_row,:)-cdisp_nosort(shaker_node_row,:)),d);

%	define reduced system using sorted modes

	sysforce_sort = ss(a_sort,b_sort,mn2gm_conversion*kspring*(cdisp_sort(tip_node_row,:)-cdisp_sort(shaker_node_row,:)),d);

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

	[magforce,phsforce] = bode(sysforce,frad);

	[magforce_nosort,phsforce_nosort] = bode(sysforce_nosort,frad);

	[magforce_sort,phsforce_sort] = bode(sysforce_sort,frad);

%	start plotting

	if  num_modes_used == num_modes_total

%	plot all modes included response

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

	hold on

	max_modes_plot = num_modes_used;

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

		amode = a_nosort(index-1:index,index-1:index);

		bmode = b_nosort(index-1:index);

		cmode_shaker = cdisp_nosort(1,index-1:index);

		cmode_tip = cdisp_nosort(numdof,index-1:index);

		dmode = [0];

		sysforce_mode = ss(amode,bmode,mn2gm_conversion*kspring*(cmode_tip - cmode_shaker),dmode);

   		[magforce_mode,phsforce_mode]=bode(sysforce_mode,frad)  ;

	   	loglog(f,magforce_mode(1,:),'k-')
	
	end

	disp('execution paused to display figure, "enter" to continue'); pause
	
	hold off

%	now use lsim to calculate force due to a 0.002 sec half-sine 100g shock pulse