act8.m

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

	for  cnt = 1:num_modes_total
	
		b(2*cnt) = f_principal(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_total,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

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

         	cvel(row,col) = 0;
   
      	end

   	end

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

      	for  row = 1:numdof

         	c_disp(row,col) = 0;

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

      	end

   	end

%	all modes included 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 state space systems with the "ss" command, outputs are the two gap displacements

%	define unsorted all modes included system

	sys = ss(a,b,c_disp(7:8,:),d);

%	define sorted all modes included system

	sys_sort = ss(a_sort,b_sort,cdisp_sort(7:8,:),d);

%	define sorted reduced system

	a_sort_red = a_sort(1:num_states_used,1:num_states_used);
	
	b_sort_red = b_sort(1:num_states_used);

	cdisp_sort_red = cdisp_sort(7:8,1:num_states_used);
	
	sys_sort_red = ss(a_sort_red,b_sort_red,cdisp_sort_red,d);

%	define modred "mdc" reduced system, modred "del" option same as sorted reduced above

	states_del = (2*num_modes_used+1):2*num_modes_total;

	sys_mdc = modred(sys_sort,states_del,'mdc');

	sys_mdc_nosort = modred(sys,[17:100],'mdc');

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

	[mag,phs] = bode(sys,frad);

	[mag_sort_red,phs_sort_red] = bode(sys_sort_red,frad);

	[mag_mdc,phs_mdc]=bode(sys_mdc,frad)  ;

	[mag_mdc_nosort,phs_mdc_nosort]=bode(sys_mdc_nosort,frad)  ;

%	convert magnitude to db

	magdb = 20*log10(mag);

	mag_sort_reddb = 20*log10(mag_sort_red);

	mag_mdcdb = 20*log10(mag_mdc);

%	check on discretized system aliasing

	sample_freq = input('enter sample frequency, khz, default 20 khz ...  ');
		
	if  isempty(sample_freq)
		
		sample_freq = 20;

	end

	nyquist_freq = sample_freq/2;

	disp(['Nyquist frequency is ',num2str(nyquist_freq),' khz']);

	ts = 1/(1000*sample_freq);
 
   	freqdlo = 500;

   	freqdhi = 1000*nyquist_freq;	% only take frequency response to nyquist_freq

   	fdlo=log10(freqdlo) ;
   	fdhi=log10(freqdhi) ;
 
   	fd=logspace(fdlo,fdhi,400) ;
   	fdrad=fd*2*pi ;

	sysd = c2d(sys,ts);
	
	[magd,phsd] = bode(sysd,fdrad);
	
	magddb = 20*log10(magd);		

%	start plotting

%	plot all modes included response

	loglog(f,mag(index_out,:),'k.-')
	title([headstr ', gap displacement, all ',num2str(num_modes_total),' modes included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 -inf 1e-4])
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	hold on

	max_modes_plot = num_modes_total;

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = b(index-1:index);

		cmode = c_disp(7:8,index-1:index);

		dmode = [0];

		sys_mode = ss(amode,bmode,cmode,dmode);

   		[mag_mode,phs_mode]=bode(sys_mode,frad)  ;

		mag_modedb = 20*log10(mag_mode);

		loglog(f,mag_mode(index_out,:),'k-')
	
	end

	axis([500 25000 -inf 1e-4])

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

	loglog(f,mag(index_out,:),'k-',fd,magd(index_out,:),'k.-')
	title([headstr ', gap displacement, all ',num2str(num_modes_total), ...
			' modes included, Nyquist frequency ',num2str(nyquist_freq),' hz'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	legend('continuous','discrete')
	axis([500 25000 1e-8 1e-4])
   	grid off

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

	if  num_modes_used < num_modes_total		% calculate and plot reduced models
	
%	sorted modal truncation

	loglog(f,mag(index_out,:),'k-',f,mag_sort_red(index_out,:),'k.-')
   	title([headstr ', sorted modal truncation:  gap displacement, first ', ...
   						num2str(num_modes_used),' modes included'])
	legend('all modes','sorted partial modes',3)
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-8 1e-4])
   	grid off

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

	hold on

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = b_sort(index-1:index);

		cmode = cdisp_sort(7:8,index-1:index);

		dmode = [0];

		sys_mode = ss(amode,bmode,cmode,dmode);

	   	[mag_mode,phs_mode]=bode(sys_mode,frad)  ;

   		loglog(f,mag_mode(index_out,:),'k-')
	
	end

	axis([500 25000 -inf 1e-4])

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

	hold off

%	modred using 'mdc'

	loglog(f,mag(index_out,:),'k-',f,mag_mdc(index_out,:),'k.-')
	title([headstr ', reduced matched dc gain:  gap displacement, first ', ...
						num2str(num_modes_used),' sorted modes included'])
	legend('all modes','reduced mdc',3)
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-8 1e-4])
   	grid off

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

	hold on

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = b_sort(index-1:index);

		cmode = cdisp_sort(7:8,index-1:index);

		dmode = [0];

		sys_mode = ss(amode,bmode,cmode,dmode);

	   	[mag_mode,phs_mode]=bode(sys_mode,frad)  ;

   		loglog(f,mag_mode(index_out,:),'k-')
	
	end

	axis([500 25000 -inf 1e-4])

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

	hold off

%	modred using 'mdc' with unsorted modes

	loglog(f,mag(index_out,:),'k-',f,mag_mdc_nosort(index_out,:),'k.-')
	title([headstr ', reduced unsorted matched dc gain:  gap displacement, first ', ...
						num2str(num_modes_used),' sorted modes included'])
	legend('all modes','reduced mdc',3)
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-8 1e-4])
   	grid off

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

	hold on

	for  pcnt = 1:num_modes_used

		index = 2*pcnt;

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

		bmode = b(index-1:index);

		cmode = c_disp(7:8,index-1:index);

		dmode = [0];

		sys_mode = ss(amode,bmode,cmode,dmode);

	   	[mag_mode,phs_mode]=bode(sys_mode,frad)  ;

   		loglog(f,mag_mode(index_out,:),'k-')
	
	end

	axis([500 25000 -inf 1e-4])

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

	hold off

	end

%	save the workspace for use in balred.m 
 
	save act8_data