cantbeam_ss_modred.m

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

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = b(index-1:index);

		cmode = cdisp(numdof,index-1:index);

		dmode = [0];

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

   		[magdisptip_mode,phsdisptip_mode]=bode(sysdisptip_mode,frad)  ;

   		magdisptip_modedb = 20*log10(magdisptip_mode);

	   	semilogx(f,magdisptip_modedb(1,:),'k-')
	
	end

	dc_gain_freq = freqlo*ones(size(freqnew));

	semilogx(dc_gain_freq(1:num_modes_used),20*log10(dc_gain(1:num_modes_used)),'ko:')

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

	hold off

%	now use lsim to calculate step response to a unit force

	ttotal = 0.1;

	t = linspace(0,ttotal,200);

	u = ones(size(t));

	[disptip,ts] = lsim(sysdisptip,u,t);

	plot(ts,disptip,'k-')
   	title(['tip disp for mid-length step force, all ',num2str(num_modes_used),' modes included'])
	xlabel('time, sec')
	ylabel('displacement, mm')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

else	

%	plot unsorted modal truncation

	semilogx(f,magdisptipdb(1,:),'k-',f,magdisptipdb_nosort(1,:),'k.-')
   	title(['unsorted modal truncation:  cantilever tip displacement for mid-length force, first ',num2str(num_modes_used),' modes included'])
	legend('all modes','unsorted partial modes',3)

	dcgain_error_percent_nosort = 100*(magdisptip_nosort(1) - magdisptip(1))/magdisptip(1)

   	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_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 = cdisp_nosort(numdof,index-1:index);

		dmode = [0];

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

	   	[magdisptip_mode,phsdisptip_mode]=bode(sysdisptip_mode,frad)  ;

		magdisptip_modedb = 20*log10(magdisptip_mode);

	   	semilogx(f,magdisptip_modedb(1,:),'k-')
	
	end

	dc_gain_freq_nosort = freqlo*ones(size(freqnew_nosort));

	semilogx(dc_gain_freq_nosort(1:num_modes_used),20*log10(dc_gain_nosort(1:num_modes_used)),'ko:')

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

	hold off

%	plot sorted modal truncation

	semilogx(f,magdisptipdb(1,:),'k-',f,magdisptipdb_sort(1,:),'k.-')
   	title(['sorted modal truncation:  cantilever tip displacement for mid-length force, first ',num2str(num_modes_used),' modes included'])
	legend('all modes','sorted partial modes',3)

	dcgain_error_percent_sort = 100*(magdisptip_sort(1) - magdisptip(1))/magdisptip(1)

   	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_used;

	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(numdof,index-1:index);

		dmode = [0];

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

	   	[magdisptip_mode,phsdisptip_mode]=bode(sysdisptip_mode,frad)  ;

		magdisptip_modedb = 20*log10(magdisptip_mode);

   		semilogx(f,magdisptip_modedb(1,:),'k-')
	
	end

	dc_gain_freq_sort = freqlo*ones(size(freqnew_nosort));

	semilogx(dc_gain_freq_sort(1:num_modes_used),20*log10(dc_gain_sort(1:num_modes_used)),'ko:')

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

	hold off

%	now use lsim to calculate step response to a unit force

	ttotal = 0.1;

	t = linspace(0,ttotal,200);

	u = ones(size(t));

	[disptip,ts] = lsim(sysdisptip,u,t);

	[disptip_nosort,ts_nosort] = lsim(sysdisptip_nosort,u,t);

	[disptip_sort,ts_sort] = lsim(sysdisptip_sort,u,t);

	plot(ts,disptip,'k-',ts_nosort,disptip_nosort,'k+-',ts_sort,disptip_sort,'k.-')
   	title(['tip disp for mid-length step force, first ',num2str(num_modes_used),' modes included'])
	legend('all modes','unsorted partial modes','sorted partial modes')
	xlabel('time, sec')
	ylabel('displacement, mm')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

%	use modred to reduce, select whether to use sorted or unsorted modes for the reduction

	modred_sort = input('modred: enter "1" to use sorted modes for reduced runs, "enter" to use unsorted ... ');

	if  isempty(modred_sort)
		modred_sort = 0
	end

	if  modred_sort == 1		%  use sorted mode order

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

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

	else						%  use original mode order

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

		freqnew = freqvec((1:num_modes_total));

	end

%	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 size of system matrix

	asize = 2*num_modes_total;

%  	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 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				

%  	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

%  	define output

   	d = [0];  %

%	define state space system for reduction, ordered defined by modred_sort

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

%	define reduced matrices using matched dc gain method "mdc"

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

	sysdisptip_mdc = modred(sysdisptip_red,states_elim,'mdc');

	[adisptip_mdc,bdisptip_mdc,cdisptip_mdc,ddisptip_mdc] = ssdata(sysdisptip_mdc);

%	define reduced matrices by eliminating high frequency states, 'del'

	sysdisptip_elim = modred(sysdisptip_red,states_elim,'del');

	[adisptip_elim,bdisptip_elim,cdisptip_elim,ddisptip_elim] = ssdata(sysdisptip_elim);

%	use "bode" command to generate magnitude/phase vectors for reduced systems

	[magdisptip_mdc,phsdisptip_mdc]=bode(sysdisptip_mdc,frad)  ;

  	[magdisptip_elim,phsdisptip_elim]=bode(sysdisptip_elim,frad)  ;

%	convert magnitude to db

	magdisptip_mdcdb = 20*log10(magdisptip_mdc);

   	magdisptip_elimdb = 20*log10(magdisptip_elim);

%	plot modred using 'elim'

	semilogx(f,magdisptipdb(1,:),'k-',f,magdisptip_elimdb(1,:),'k.-')

	if  modred_sort == 1
   		title(['reduced elimination:  tip disp for mid-length step force, first ',num2str(num_modes_used),' sorted modes included'])
		del_dcgain_error_percent_sort = 100*(magdisptip_elimdb(1) - magdisptip(1))/magdisptip(1)
	else
   		title(['reduced elimination:  tip disp for mid-length step force, first ',num2str(num_modes_used),' unsorted modes included'])
		del_dcgain_error_percent_nosort = 100*(magdisptip_elimdb(1) - magdisptip(1))/magdisptip(1)
	end
	
	legend('all modes','reduced elim',3)

   	xlabel('Frequency, hz')
	ylabel('Magnitude, db mm')
   	grid off

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

	hold on

%	now plot the overlay of the tip displacement magnitude with each mode contribution

	max_modes_plot = num_modes_used;

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = bdisptip_elim(index-1:index);

		cmode = cdisptip_elim(1,index-1:index);

		dmode = [0];

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

	   	[magdisptip_mode,phsdisptip_mode]=bode(sysdisptip_mode,frad)  ;

		magdisptip_modedb = 20*log10(magdisptip_mode);

   		semilogx(f,magdisptip_modedb(1,:),'k-')
	
	end

	dc_gain_freq_sort = freqlo*ones(size(freqnew_nosort));

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

	hold off

%	modred using 'mdc'

	semilogx(f,magdisptipdb(1,:),'k-',f,magdisptip_mdcdb(1,:),'k.-')

	if  modred_sort == 1
   		title(['reduced matched dc gain:  tip disp for mid-length step force, first ',num2str(num_modes_used),' sorted modes included'])
		mdc_dcgain_error_percent_sort = 100*(magdisptip_mdcdb(1) - magdisptip(1))/magdisptip(1)
	else
   		title(['reduced matched dc gain:  tip disp for mid-length step force, first ',num2str(num_modes_used),' unsorted modes included'])
		mdc_dcgain_error_percent_nosort = 100*(magdisptip_mdcdb(1) - magdisptip(1))/magdisptip(1)
	end

	legend('all modes','reduced mdc',3)


   	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_used;

	for  pcnt = 1:max_modes_plot

		index = 2*pcnt;

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

		bmode = bdisptip_mdc(index-1:index);

		cmode = cdisptip_mdc(1,index-1:index);

		dmode = [0];

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

	   	[magdisptip_mode,phsdisptip_mode]=bode(sysdisptip_mode,frad)  ;

		magdisptip_modedb = 20*log10(magdisptip_mode);

	   	semilogx(f,magdisptip_modedb(1,:),'k-')
	
	end

	dc_gain_freq_nosort = freqlo*ones(size(freqnew_nosort));

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

	hold off

%	now use lsim to calculate step response to a unit force

	[disptip,ts] = lsim(sysdisptip,u,t);

	[disptip_elim,ts_elim] = lsim(sysdisptip_elim,u,t);

	[disptip_mdc,ts_mdc] = lsim(sysdisptip_mdc,u,t);

	plot(ts,disptip,'k-',ts_mdc,disptip_mdc,'k.-',ts_elim,disptip_elim,'k+-')

	if  modred_sort == 1
   		title(['modred cantilever tip disp for mid-length step force, first ',num2str(num_modes_used),' sorted modes included'])
	else
   		title(['modred cantilever tip disp for mid-length step force, first ',num2str(num_modes_used),' unsorted modes included'])
	end

	legend('all modes','reduced - mdc','reduced - elim')
	xlabel('time, sec')
	ylabel('displacement, mm')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	end