function 
 
 
   plotmass 
 
 
   % Kyle Wilson 10-2-02 
 
 
   % Stirling Cycle Machine Analysis ME 589 
 
 
   % Particle Trajectory Map 
 
 
   % Equations from Organ's "'Natural' coordinates for analysis of the practical 
 
 
   % Stirling cycle" and Oegik  Soegihardjo's 1993 project on the same topic 
 
 
   % Modified by Israel Urieli (11/27/2010) to obtain correct phase advance 
 
 
   % angle alpha subsequent to error determined by Zack Alexy (March 2010) 
 
 
   % Global values from define program 
 
 
   global 
 
 
   vclc vcle 
 
 
   % compression,expansion clearance vols [m^3] 
 
 
   global 
 
 
   vswc vswe 
 
 
   % compression, expansion swept volumes [m^3] 
 
 
   global 
 
 
   alpha 
 
 
   % phase angle advance of expansion space [radians] 
 
 
   global 
 
 
   vk 
 
 
   % cooler void volume [m^3] 
 
 
   global 
 
 
   vh 
 
 
   % heater void volume [m^3] 
 
 
   global 
 
 
   vr 
 
 
   % regen void volume [m^3] 
 
 
   global 
 
 
   pmean 
 
 
   % mean (charge) pressure [Pa] 
 
 
   global 
 
 
   tk tr th 
 
 
   % cooler, regenerator, heater temperatures [K] 
 
 
   NT = th/tk; 
 
 
   % Temperature ratio 
 
 
   Vref = vswe; 
 
 
   % Reference volume (m^3) 
 
 
   %% Fixed reduced volumes (dimensionless) 
 
 
   vswe_r = (vswe/Vref)/NT; 
 
 
   % Reduced expansion swept volume 
 
 
   vcle_r = (vcle/Vref)/NT; 
 
 
   % Reduced expansion clearance volume 
 
 
   vh_r = (vh/Vref)/NT; 
 
 
   % Reduced heater void volume 
 
 
   vr_r = (vr/Vref)*log(NT)/(NT-1); 
 
 
   % Reduced regenerator void volume 
 
 
   vk_r = (vk/Vref); 
 
 
   % Reduced cooler void volume 
 
 
   vswc_r = (vswc/Vref); 
 
 
   % Reduced compression swept volume 
 
 
   vclc_r = (vclc/Vref); 
 
 
   % Reduced compression clearance volume 
 
 
   %% Phase domain 
 
 
   angi = 0; 
 
 
   angf = 2*pi; 
 
 
   dang = 0.1; 
 
 
   ang = [angi:dang:angf]; 
 
 
   n = size(ang); 
 
 
   %% Reduced volume variations 
 
 
   for 
 
 
   i = 1:n(2) 
 
 
   deg(i) = ang(i)*180/pi; 
 
 
   Ve(i) = (vswe/2)*(1- cos(ang(i))); 
 
 
   % Expansion volume vs phase 
 
 
   Vc(i) = (vswc/2)*(1+ cos(ang(i) - alpha)); 
 
 
   % Compression volume vs phase 
 
 
   ve(i) = (Ve(i)/Vref)/NT; 
 
 
   % Reduced expansion vs phase 
 
 
   vc(i) = Vc(i)/Vref; 
 
 
   % Reduced compression vs phase 
 
 
   vt(i) = vswe_r + vcle_r + vh_r + vr_r + vk_r + vclc_r + vc(i); 
 
 
   % Total volume vs phase 
 
 
   end 
 
 
   figure 
 
 
   step = 30; 
 
 
   for 
 
 
   m = 1:step-1 
 
 
   for 
 
 
   i = 1:n(2) 
 
 
   v(i) = ve(i) + (m/step)*(vt(i)-ve(i)); 
 
   % Reduced volume segments 
 
 
 
   end 
 
 
   hold on 
 
 
   plot(v,deg, 
 
 
   'k:' 
 
 
   ) 
 
 
   end 
 
 
   hold on 
 
 
   plot(ve,deg, 
 
 
   'k' 
 
 
   ) 
 
 
   plot(vt,deg, 
 
 
   'k' 
 
 
   ) 
 
 
   %% Vertical lines 
 
 
   L1 = vswe_r; 
 
 
   % Boundary of reduced expansion swept volume 
 
 
   L2 = L1 + vcle_r; 
 
 
   % Boundary of reduced expansion clearance volume 
 
 
   L3 = L2 + vh_r; 
 
 
   % Boundary of reduced heater void volume 
 
 
   L4 = L3 + vr_r; 
 
 
   % Boundary of reduced regenerator void volume 
 
 
   L5 = L4 + vk_r; 
 
 
   % Boundary of reduced cooler void volume 
 
 
   L6 = min(vt); 
 
 
   % Boundary of reduced expansion swept volume 
 
 
   point1 = [L1;L1]; 
 
 
   % Preparing for plot 
 
 
   point2 = [L2;L2]; 
 
 
   point3 = [L3;L3]; 
 
 
   point4 = [L4;L4]; 
 
 
   point5 = [L5;L5]; 
 
 
   point6 = [L6;L6]; 
 
 
   point = [0;deg(n(2))]; 
 
 
   plot(point1, point, 
 
 
   'r--' 
 
 
   , point2, point, 
 
 
   'r--' 
 
 
   , point3, point, 
 
 
   'g--' 
 
 
   ) 
 
 
   plot(point4, point, 
 
 
   'g--' 
 
 
   , point5, point, 
 
 
   'b--' 
 
 
   , point6, point, 
 
 
   'b--' 
 
 
   ) 
 
 
   axis([0 max(vt) 0 deg(n(2))]) 
 
 
   xlabel( 
 
 
   'Reduced volume' 
 
 
   ) 
 
 
   ylabel( 
 
 
   'Crank Angle (deg)' 
 
 
   ) 
 
 
   title( 
 
 
   'Particle mass plot' 
 
 
   ) 
 
 
   hold off