act8.m

振动仿真 vibration simulation using matlab and ansys!

%	act8.m		

   	clear all;
	
	hold off;

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

%	the output for the ANSYS run is the following dof's

%  dof	node	dir		 where
%
%	 1	    22	ux - radial, top head gap
%    2   10022	ux - radial, bottom head gap
%    3   24061	ux - radial, coil
%    4   24066	ux - radial, coil
%    5   24082	ux - radial, coil
%    6   24087	ux - radial, coil
%	 7  	22	uy - circumferential, top head gap
%    8   10022	uy - circumferential, bottom head gap
%    9   24061	uy - circumferential, coil
%   10   24066	uy - circumferential, coil
%   11   24082	uy - circumferential, coil
%   12   24087	uy - circumferential, coil

	load actrl_eig;

	[numdof,num_modes_total] = size(evr);

	freqvec(1) = 0;					% set frequency of rigid body mode to zero
	
	xn = evr;
 
%   calculate the dc amplitude of the displacement of each mode by 
%   multiplying the composite forcing function by the output row

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

%	define frequency range for frequency response

   	freqlo = 501;

   	freqhi = 25000;

   	flo=log10(freqlo) ;
   	fhi=log10(freqhi) ;
 
   	f=logspace(flo,fhi,300) ;
   	frad=f*2*pi ;
 
%	define radial and circumferential forces applied at four coil force nodes
%	"x" is radial, "y" is circumferential, total force is unity

	n24061fx = 0.25*sin(9.1148*pi/180);
	n24061fy = 0.25*cos(9.1148*pi/180);

	n24066fx = 0.25*sin(15.1657*pi/180);
	n24066fy = 0.25*cos(15.1657*pi/180);
			  
	n24082fx = -0.25*sin(15.1657*pi/180);
	n24082fy = 0.25*cos(15.1657*pi/180);

	n24087fx = -0.25*sin(9.1148*pi/180);
	n24087fy = 0.25*cos(9.1148*pi/180);

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

	f_physical = [    	0
					  	0
					n24061fx
					n24066fx
					n24082fx
					n24087fx
						0
						0
					n24061fy
					n24066fy
					n24082fy
					n24087fy  ];

% 	define composite forcing function, force applied to each node times	eigenvector value
%	for that node

	force = f_physical'*xn;

%	choose which head to use for frequency responses

	head = input('enter "0" default for head 0 or "1" for head 1 ...  ');

	if  isempty(head)
		head = 0;
	end

%	prompt for uniform or variable zeta

	zeta_type = input('enter "1" to read in damping vector (zetain.m) or "enter" for uniform damping ...  ');

	if  (isempty(zeta_type))

		zeta_type = 0;
		
	   	zeta_uniform = input('enter value for uniform damping, .005 is 0.5% of critical (default) ... ');
	
    	if (isempty(zeta_uniform))
    		zeta_uniform = 0.005;
	    end

		zeta_unsort = zeta_uniform*ones(num_modes_total,1);

		gainstr = 'dc gain';

	else

		zetain;		% read in zeta_unsort damping vector from zetain.m file

		gainstr = 'peak gain';

	end

	if  length(zeta_unsort) ~= num_modes_total

		error(['error - zeta_unsort vector has ',num2str(length(zeta_unsort)),' entries instead of ', ...
				num2str(num_modes_total)]);

	end

%	calculate dc gains if uniform damping, peak gains if non-uniform

	if  zeta_type == 0			% dc gain

		gain_h0 = abs([force(1)*xn(8,1)/frad(1) ...
					force(2:num_modes_total).*xn(8,2:num_modes_total) ...
					./omega2(2:num_modes_total)]);
	
		gain_h1 = abs([force(1)*xn(7,1)/frad(1) ...
					force(2:num_modes_total).*xn(7,2:num_modes_total) ...
					./omega2(2:num_modes_total)]);

	elseif  zeta_type == 1		% peak gain

		gain_h0 = abs([force(1)*xn(8,1)/frad(1) ...
					force(2:num_modes_total).*xn(8,2:num_modes_total) ...
					./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);
	
		gain_h1 = abs([force(1)*xn(7,1)/frad(1) ...
					force(2:num_modes_total).*xn(7,2:num_modes_total) ...
					./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);

	end


%	sort gains, keeping track of original and new indices so can rearrange 
%	eigenvalues and eigenvectors

	[gain_h0_sort,index_h0_sort] = sort(gain_h0);

	[gain_h1_sort,index_h1_sort] = sort(gain_h1);

	gain_h0_sort = fliplr(gain_h0_sort);		% max to min

	gain_h1_sort = fliplr(gain_h1_sort);		% max to min

	index_h0_sort = fliplr(index_h0_sort)				% max to min indices

	index_h1_sort = fliplr(index_h1_sort)				% max to min indices

    index_orig = 1:num_modes_total;

	if  head == 0

		index_sort = index_h0_sort;

		headstr = 'head 0';

		index_out = 2;

	elseif  head == 1

		index_sort = index_h1_sort;

		headstr = 'head 1';

		index_out = 1;

	end

%	plot results
 
   	semilogy(index_orig(2:num_modes_total),freqvec(2:num_modes_total),'k-');
   	title('frequency versus mode number')
   	xlabel('mode number')
   	ylabel('frequency, hz')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0,'k-',index_orig,gain_h1,'k.-')
   	title([gainstr ' of each mode contribution versus mode number'])
   	xlabel('mode number')
   	ylabel('dc value')
	legend('head 0','head 1')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	loglog(freqvec(2:num_modes_total),gain_h0(2:num_modes_total),'k-', ...
   		   freqvec(2:num_modes_total),gain_h1(2:num_modes_total),'k.-')
	title([gainstr ' of each mode contribution versus frequency'])
 	xlabel('frequency, hz')
	ylabel('dc value')
	legend('head 0','head 1')
	axis([500 25000 -inf 1e-4])
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_sort,'k-',index_orig,gain_h1_sort,'k.-')
   	title([gainstr ' of each mode versus number of modes included'])
   	xlabel('modes included')
   	ylabel('sorted dc value')
	legend('head 0','head 1')
   	grid off

%	choose number of modes to use based on ranking of dc gain values

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

	num_states_used = 2*num_modes_used;

%	define eigenvalues and eigenvectors for unsorted and sorted modes

%  	all modes included model, use original order

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

	freqnew = freqvec((1:num_modes_total));

	zeta = zeta_unsort;

%	all modes included, sorted

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

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

	zeta_sort = zeta_unsort(index_sort(1:num_modes_total));

%	define variables for all modes included system matrix, a

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

	w2 = w.^2;				

    zw = 2*zeta_unsort.*w;

%	define variables for all modes included sorted system matrix, a_sort

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

	w2_sort = w_sort.^2;				

    zw_sort = 2*zeta_sort.*w_sort;

%	define size of system matrix

	asize = 2*num_modes_total;

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

%  	setup system matrix for all modes included model

   	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 system matrix for sorted all modes included model

   	a_sort = zeros(asize);

   	for  col = 2:2:asize

      	row = col-1;

      	a_sort(row,col) = 1;

   	end

   	for  col = 1:2:asize

      	row = col+1;

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

   	end

   	for  col = 2:2:asize

      	row = col;

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

   	end

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

%	now setup the principal force vector for the three cases, all modes, 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);