cantbeam_ss_shkr_modred.m

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

	ttotal = 0.03;

	shock_amplitude = 100;

	pulse_width = input('enter half-sine shock pulse width, sec, default is 0.002 ...  ');

	if  isempty(pulse_width)
		pulse_width = 0.002;
	end

	t = linspace(0,ttotal,1000);

  	dt = t(2) - t(1);

	for  cnt = 1:length(t)

 		if  t(cnt) < pulse_width  
 	
        	u(cnt) = shock_amplitude*sin(2*pi*(1/(2*pulse_width))*t(cnt));

	 	else

    		u(cnt) = 0;

     	end

	end

	plot(t,u,'k-')
	title('acceleration of shaker mass')
	xlabel('time, sec')
	ylabel('acceleration, g')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	[force,ts] = lsim(sysforce,u,t);

	plot(ts,force,'k-')
   	title(['cantilever tip force for ',num2str(shock_amplitude),'g, ',num2str(pulse_width) ...
   			,' sec input, all ',num2str(num_modes_used),' modes included'])
	xlabel('time, sec')
	ylabel('Force, gm')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	peak_force = max(abs(force))

	else
	
%	unsorted modal truncation			

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

	dcgain_error_percent_nosort = 100*(magforce_nosort(1) - magforce(1))/magforce(1)

   	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

%	sorted modal truncation

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

	dcgain_error_percent_sort = 100*(magforce_sort(1) - magforce(1))/magforce(1)

   	xlabel('Frequency, hz')
	ylabel('Force, gm')
   	grid off

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

	hold on

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

	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_shaker = cdisp_sort(1,index-1:index);

		cmode_tip = cdisp_sort(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

	ttotal = 0.03;

	shock_amplitude = 100;

	pulse_width = input('enter half-sine shock pulse width, sec, default is 0.002 ...  ');

	if  isempty(pulse_width)
		pulse_width = 0.002;
	end

	t = linspace(0,ttotal,1000);

  	dt = t(2) - t(1);

	for  cnt = 1:length(t)

     	if  t(cnt) < pulse_width  
 	
	    	u(cnt) = shock_amplitude*sin(2*pi*(1/(2*pulse_width))*t(cnt));

     	else

	    	u(cnt) = 0;

 		end

  	end

	plot(t,u,'k-')
	title('acceleration of shaker mass')
	xlabel('time, sec')
	ylabel('acceleration, g')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause
	
	[force,ts] = lsim(sysforce,u,t);

	[force_nosort,ts_nosort] = lsim(sysforce_nosort,u,t);

	[force_sort,ts_sort] = lsim(sysforce_sort,u,t);

	plot(ts,force,'k-',ts_nosort,force_nosort,'k+:',ts_sort,force_sort,'k.-')
   	title(['cantilever tip force for ',num2str(shock_amplitude),'g, ',num2str(pulse_width) ...
   			,' sec input, ',num2str(num_modes_used),' modes included'])
	legend('all modes','unsorted partial modes','sorted partial modes',4)
	xlabel('time, sec')
	ylabel('Force, gm')
	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	max_force = max(abs(force));

	max_force_nosort = max(abs(force_nosort));
	max_force_sort = max(abs(force_sort));

	peak_error_nosort_percent = 100*(max_force_nosort - max_force)/max_force
	peak_error_sort_percent = 100*(max_force_sort - max_force)/max_force

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

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

%  	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

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

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

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

	sysforce_mdc = modred(sysforce_red,states_elim,'mdc');

	[aforce_mdc,bforce_mdc,cforce_mdc,dforce_mdc] = ssdata(sysforce_mdc);

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

	sysforce_elim = modred(sysforce_red,states_elim,'del');

	[aforce_elim,bforce_elim,cforce_elim,dforce_elim] = ssdata(sysforce_elim);

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

	[magforce_mdc,phsforce_mdc]=bode(sysforce_mdc,frad)  ;

  	[magforce_elim,phsforce_elim]=bode(sysforce_elim,frad)  ;

%	convert magnitude to db

	magforce_mdcdb = 20*log10(magforce_mdc);

   	magforce_elimdb = 20*log10(magforce_elim);

%	start plotting

%	modred using 'elim'

	loglog(f,magforce(1,:),'k-',f,magforce_elim(1,:),'k.-')

	if  modred_sort == 1
   		title(['reduced elimination:  cantilever tip force for mid-length force, first ',num2str(num_modes_used),' sorted modes included'])
		dcgain_error_percent_elim_sort = 100*(magforce_elim(1) - magforce(1))/magforce(1)
	else
   		title(['reduced elimination:  cantilever tip force for mid-length force, first ',num2str(num_modes_used),' unsorted modes included'])
		dcgain_error_percent_elim_nosort = 100*(magforce_elim(1) - magforce(1))/magforce(1)
	end
	
	legend('all modes','reduced elimination',3)

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

		bmode = bforce_elim(index-1:index);

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

		dmode = [0];

		sysforce_mode = ss(amode,bmode,cmode,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

%	modred using 'mdc'

	loglog(f,magforce(1,:),'k-',f,magforce_mdc(1,:),'k.-')

	if  modred_sort == 1
   		title(['reduced matched dc gain:  cantilever tip force for mid-length force, first ',num2str(num_modes_used),' sorted modes included'])
		dcgain_error_percent_mdc_sort = 100*(magforce_mdc(1) - magforce(1))/magforce(1)
	else
   		title(['reduced matched dc gain:  cantilever tip force for mid-length force, first ',num2str(num_modes_used),' unsorted modes included'])
		dcgain_error_percent_mdc_nosort = 100*(magforce_mdc(1) - magforce(1))/magforce(1)
	end

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


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

		bmode = bforce_mdc(index-1:index);

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

		dmode = [0];

		sysforce_mode = ss(amode,bmode,cmode,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

	[force_mdc,ts_mdc] = lsim(sysforce_mdc,u,t);

	[force_elim,ts_elim] = lsim(sysforce_elim,u,t);

	plot(ts,force,'k-',ts_mdc,force_mdc,'k.-',ts_elim,force_elim,'k+-')

	if  modred_sort == 1
	   	title(['modred cantilever tip force for ',num2str(shock_amplitude),'g, ',num2str(pulse_width) ...
   				,' sec input, ',num2str(num_modes_used),' sorted modes included'])
	else
	   	title(['modred cantilever tip force for ',num2str(shock_amplitude),'g, ',num2str(pulse_width) ...
   				,' sec input, ',num2str(num_modes_used),' unsorted modes included'])
	end

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

	max_force_mdc = max(abs(force_mdc));
	max_force_elim = max(abs(force_elim));

	peak_error_mdc_percent = 100*(max_force_mdc - max_force)/max_force
	peak_error_elim_percent = 100*(max_force_elim - max_force)/max_force

	end