act8pz.m

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

	num_modes_used = 1 + osc_states_used/2;		% number of modes for overlaid plots

%	use modred to order oscillatory states from balreal to define reduced order
%	oscillatory system using both "del" and "mdc"

	rsys_delo = modred(sysob,osc_states_used+1:2*num_modes_total-2,'del');

	rsys_mdco = modred(sysob,osc_states_used+1:2*num_modes_total-2,'mdc');

%	rebuild system by appending balanced realization of oscillatory modes to rigid body mode

	[a_delo_bal,b_delo_bal,c_delo_bal,d_delo_bal] = ssdata(rsys_delo);

	a_del_bal = [   a(1:2,1:2)      zeros(2,osc_states_used)
			     zeros(osc_states_used,2)      a_delo_bal    ];
		
	b_del_bal = [b(1:2,:)
			     b_delo_bal];

	c_del_bal = [c_disp(1:2,1:2) c_delo_bal];

	rsys_del = ss(a_del_bal,b_del_bal,c_del_bal,d);

	[a_mdco_bal,b_mdco_bal,c_mdco_bal,d_mdco_bal] = ssdata(rsys_mdco);

	a_mdc_bal = [   a(1:2,1:2)      zeros(2,osc_states_used)
			     zeros(osc_states_used,2)      a_mdco_bal    ];
		
	b_mdc_bal = [b(1:2,:)
			     b_mdco_bal];

	c_mdc_bal = [c_disp(1:2,1:2) c_mdco_bal];

	rsys_mdc = ss(a_mdc_bal,b_mdc_bal,c_mdc_bal,d);

%	frequency response for unsorted system

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

%	plot original system response, output of bode command has dimensions of "i" x "j" x "k"
%	where "i" is output row, "j" is input column and "k" is the vector of frequencies

	magh0coil = mag(2,1,:);
	magh1coil = mag(1,1,:);
	magh0pz   = mag(2,2,:);
	magh1pz   = mag(1,2,:);

	loglog(f,magh0coil(1,:),'k.-',f,magh1coil(1,:),'k-')
	title(['gap displacement, all ',num2str(num_modes_total),' modes included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','head 1, coil input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0pz(1,:),'k.-',f,magh1pz(1,:),'k-')
	title(['gap displacement, all ',num2str(num_modes_total),' modes included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0 piezo input','head 1 piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0coil(1,:),'k.-',f,magh1coil(1,:),'k.-',f,magh0pz(1,:),'k-',f,magh1pz(1,:),'k-')
	title(['gap displacement, all ',num2str(num_modes_total),' modes included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','head 1, coil input','head 0 piezo input','head 1 piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

%	frequency response for balanced reduced modred "del"

	[magr_del,phsr_del] = bode(rsys_del,frad);

	magr_delh0coil = magr_del(2,1,:);
	magr_delh1coil = magr_del(1,1,:);
	magr_delh0pz   = magr_del(2,2,:);
	magr_delh1pz   = magr_del(1,2,:);

	loglog(f,magr_delh0coil(1,:),'k-',f,magr_delh1coil(1,:),'k.-',f,magr_delh0pz(1,:),'k.-',f,magr_delh1pz(1,:),'k-')
	title(['gap displacement, modred "del", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','head 1, coil input','head 0 piezo input','head 1 piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0coil(1,:),'k-',f,magr_delh0coil(1,:),'k.-')
	title(['gap displacement, modred "del", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','"del" reduced head 0, coil input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh1coil(1,:),'k-',f,magr_delh1coil(1,:),'k.-')
	title(['gap displacement, modred "del", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 1, coil input','"del" reduced head 1, coil input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0pz(1,:),'k-',f,magr_delh0pz(1,:),'k.-')
	title(['gap displacement, modred "del", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, piezo input','"del" reduced head 0, piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh1pz(1,:),'k-',f,magr_delh1pz(1,:),'k.-')
	title(['gap displacement, modred "del", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 1, piezo input','"del" reduced head 1, piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

%	frequency response for balanced reduced modred "mdc"

	[magr_mdc,phsr_mdc] = bode(rsys_mdc,frad);

	magr_mdch0coil = magr_mdc(2,1,:);
	magr_mdch1coil = magr_mdc(1,1,:);
	magr_mdch0pz   = magr_mdc(2,2,:);
	magr_mdch1pz   = magr_mdc(1,2,:);

	loglog(f,magr_mdch0coil(1,:),'k-',f,magr_mdch1coil(1,:),'k.-',f,magr_mdch0pz(1,:),'k.-',f,magr_mdch1pz(1,:),'k-')
	title(['gap displacement, modred "mdc", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','head 1, coil input','head 0 piezo input','head 1 piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0coil(1,:),'k-',f,magr_mdch0coil(1,:),'k.-')
	title(['gap displacement, modred "mdc", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, coil input','"mdc" reduced head 0, coil input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh1coil(1,:),'k-',f,magr_mdch1coil(1,:),'k.-')
	title(['gap displacement, modred "mdc", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 1, coil input','"mdc" reduced head 1, coil input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh0pz(1,:),'k-',f,magr_mdch0pz(1,:),'k.-')
	title(['gap displacement, modred "mdc", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 0, piezo input','"mdc" reduced head 0, piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

	loglog(f,magh1pz(1,:),'k-',f,magr_mdch1pz(1,:),'k.-')
	title(['gap displacement, modred "mdc", ',num2str(osc_states_used),' oscillatory states included'])
   	xlabel('Frequency, hz')
	ylabel('Magnitude, mm')
	axis([500 25000 1e-9 2e-4])
	legend('head 1, piezo input','"mdc" reduced head 1, piezo input',3)
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

%	calculate impulse responses

	ttotal = 0.0025;

	t = linspace(0,ttotal,400)';

	[disp_syso,t_syso] = impulse(syso,t);

	[disp_rsys_delo,t_rsys_delo] = impulse(rsys_delo,t);

	[disp_rsys_mdco,t_rsys_mdco] = impulse(rsys_mdco,t);

	disph0coil = disp_syso(:,2,1);
	disph1coil = disp_syso(:,1,1);
	disph0pz   = disp_syso(:,2,2);
	disph1pz   = disp_syso(:,1,2);

	dispr_delh0coil = disp_rsys_delo(:,2,1);
	dispr_delh1coil = disp_rsys_delo(:,1,1);
	dispr_delh0pz   = disp_rsys_delo(:,2,2);
	dispr_delh1pz   = disp_rsys_delo(:,1,2);

	dispr_mdch0coil = disp_rsys_mdco(:,2,1);
	dispr_mdch1coil = disp_rsys_mdco(:,1,1);
	dispr_mdch0pz   = disp_rsys_mdco(:,2,2);
	dispr_mdch1pz   = disp_rsys_mdco(:,1,2);

%	build matrix of results

	dispo = [disph0coil disph1coil disph0pz disph1pz ...
			 dispr_delh0coil dispr_delh1coil dispr_delh0pz dispr_delh1pz ...
			 dispr_mdch0coil dispr_mdch1coil dispr_mdch0pz dispr_mdch1pz];

	h0coil_del_del =  dispo(:,1) - dispo(:,5);

	h1coil_del_del =  dispo(:,2) - dispo(:,6);

	h0piezo_del_del = dispo(:,3) - dispo(:,7);

	h1piezo_del_del = dispo(:,4) - dispo(:,8);

	h0coil_mdc_del =  dispo(:,1) - dispo(:,9);

	h1coil_mdc_del =  dispo(:,2) - dispo(:,10);

	h0piezo_mdc_del = dispo(:,3) - dispo(:,11);

	h1piezo_mdc_del = dispo(:,4) - dispo(:,12);

	index_h0coil_del = sqrt(sum(h0coil_del_del.*h0coil_del_del))/sqrt(sum(dispo(:,1).*dispo(:,1)));

	index_h1coil_del = sqrt(sum(h1coil_del_del.*h1coil_del_del))/sqrt(sum(dispo(:,2).*dispo(:,2)));

	index_h0piezo_del = sqrt(sum(h0piezo_del_del.*h0piezo_del_del))/sqrt(sum(dispo(:,3).*dispo(:,3)));

	index_h1piezo_del = sqrt(sum(h1piezo_del_del.*h1piezo_del_del))/sqrt(sum(dispo(:,4).*dispo(:,4)));

	index_h0coil_mdc = sqrt(sum(h0coil_mdc_del.*h0coil_mdc_del))/sqrt(sum(dispo(:,1).*dispo(:,1)));

	index_h1coil_mdc = sqrt(sum(h1coil_mdc_del.*h1coil_mdc_del))/sqrt(sum(dispo(:,2).*dispo(:,2)));

	index_h0piezo_mdc = sqrt(sum(h0piezo_mdc_del.*h0piezo_mdc_del))/sqrt(sum(dispo(:,3).*dispo(:,3)));

	index_h1piezo_mdc = sqrt(sum(h1piezo_mdc_del.*h1piezo_mdc_del))/sqrt(sum(dispo(:,4).*dispo(:,4)));

	[index_h0coil_del index_h1coil_del index_h0piezo_del index_h1piezo_del ...
		index_h0coil_mdc index_h1coil_mdc index_h0piezo_mdc index_h1piezo_mdc]

	plot(t_syso,disph0coil,'k.-',t_rsys_delo,dispr_delh0coil,'k-',t_rsys_mdco,dispr_mdch0coil,'k--')
	title(['head 0, displacement vs time, coil impulse input, ',num2str(osc_states_used),' oscillatory states included'])
	xlabel('time, sec')
	ylabel('displacement, mm')
	legend('all modes','modred del','modred mdc',4)
	grid off

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

	plot(t_syso,disph1coil,'k.-',t_rsys_delo,dispr_delh1coil,'k-',t_rsys_mdco,dispr_mdch1coil,'k--')
	title(['head 1, displacement vs time, coil impulse input, ',num2str(osc_states_used),' oscillatory states included'])
	xlabel('time, sec')
	ylabel('displacement, mm')
	legend('all modes','modred del','modred mdc',4)
	grid off

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

	plot(t_syso,disph0pz,'k.-',t_rsys_delo,dispr_delh0pz,'k-',t_rsys_mdco,dispr_mdch0pz,'k--')
	title(['head 0, displacement vs time, piezo impulse input, ',num2str(osc_states_used),' oscillatory states included'])
	xlabel('time, sec')
	ylabel('displacement, mm')
	legend('all modes','modred del','modred mdc',4)
	grid off

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

	plot(t_syso,disph1pz,'k.-',t_rsys_delo,dispr_delh1pz,'k-',t_rsys_mdco,dispr_mdch1pz,'k--')
	title(['head 1, displacement vs time, piezo impulse input, ',num2str(osc_states_used),' oscillatory states included'])
	xlabel('time, sec')
	ylabel('displacement, mm')
	legend('all modes','modred del','modred mdc',4)
	grid off

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

%		  states h0cd	h1cd   h0pd   h1pd   h0cm   h1cm   h0pm   h1pm

	error = [10	0.1081 0.1075 0.4162 0.3963 0.1081 0.1075 0.4165 0.3964

			 12	0.1079 0.1072 0.3154 0.3058 0.1079 0.1073 0.3157 0.3061

			 16	0.1075 0.1070 0.1393 0.1421 0.1074 0.1070 0.1393 0.1419

			 20	0.0395 0.0425 0.1391 0.1410 0.0397 0.0425 0.1391 0.1411

			 24	0.0363 0.0374 0.0839 0.0873 0.0463 0.0473 0.0841 0.0875

			 28	0.0161 0.0178 0.0469 0.0495 0.0160 0.0191 0.0791 0.0794

			 32	0.0140 0.0142 0.0145 0.0160 0.0142 0.0143 0.0146 0.0163	];

	nmode = error(:,1)/2;
	
	error_h0coil_del = error(:,2);
	
	error_h1coil_del = error(:,3);
	
	error_h0piezo_del = error(:,4);
	
	error_h1piezo_del = error(:,5);
	
	error_h0coil_mdc = error(:,6);
	
	error_h1coil_mdc = error(:,7);
	
	error_h0piezo_mdc = error(:,8);
	
	error_h1piezo_mdc = error(:,9);
	
	plot(nmode,error_h0coil_del,'k.-',nmode,error_h0coil_mdc,'k-')
	title('head 0, coil input normalized reduction index')
	xlabel('number of modes included')
	ylabel('normalized reduction index')
	legend('modred del','modred mdc')
	axis([0 20 0 0.5])
	grid off
	
	disp('execution paused to display figure, "enter" to continue'); pause

	plot(nmode,error_h1coil_del,'k.-',nmode,error_h1coil_mdc,'k-')
	title('head 1, coil input normalized reduction index')
	xlabel('number of modes included')
	ylabel('normalized reduction index')
	legend('modred del','modred mdc')
	axis([0 20 0 0.5])
	grid off
	
	disp('execution paused to display figure, "enter" to continue'); pause

	plot(nmode,error_h0piezo_del,'k.-',nmode,error_h0piezo_mdc,'k-')
	title('head 0, piezo input normalized reduction index')
	xlabel('number of modes included')
	ylabel('normalized reduction index')
	legend('modred del','modred mdc')
	axis([0 20 0 0.5])
	grid off
	
	disp('execution paused to display figure, "enter" to continue'); pause

	plot(nmode,error_h1piezo_del,'k.-',nmode,error_h1piezo_mdc,'k-')
	title('head 1, piezo input normalized reduction index')
	xlabel('number of modes included')
	ylabel('normalized reduction index')
	legend('modred del','modred mdc')
	axis([0 20 0 0.5])
	grid off
	
	disp('execution paused to display figure, "enter" to continue'); pause