function 
 
 
   [var,dvar] = simple 
 
 
   % simple analysis - including heat transfer and pressure drop effects 
 
 
   % Israel Urieli, 7/22/2002 (modified 12/3/2003 for temp plots) 
 
 
   % Modified 2/7/2010 to include qrloss in hotsim & kolsim 
 
 
   % Modified 2/15/2010 logical reorganization 
 
 
   % Returned values: 
 
 
   %   var(22,37) array of variable values every 10 degrees (0 - 360) 
 
 
   %   dvar(16,37) array of derivatives every 10 degrees (0 - 360) 
 
 
   % Row indices of the var, dvar arrays: 
 
 
   TC = 1; 
 
 
   % Compression space temperature (K) 
 
 
   TE = 2; 
 
 
   % Expansion space temperature (K) 
 
 
   QK = 3; 
 
 
   % Heat transferred to the cooler (J) 
 
 
   QR = 4; 
 
 
   % Heat transferred to the regenerator (J) 
 
 
   QH = 5; 
 
 
   % Heat transferred to the heater (J) 
 
 
   WC = 6; 
 
 
   % Work done by the compression space (J) 
 
 
   WE = 7; 
 
 
   % Work done by the expansion space (J) 
 
 
   W  = 8; 
 
 
   % Total work done (WC + WE) (J) 
 
 
   P  = 9; 
 
 
   % Pressure (Pa) 
 
 
   VC = 10; 
 
 
   % Compression space volume (m^3) 
 
 
   VE = 11; 
 
 
   % Expansion space volume (m^3) 
 
 
   MC = 12; 
 
 
   % Mass of gas in the compression space (kg) 
 
 
   MK = 13; 
 
 
   % Mass of gas in the cooler (kg) 
 
 
   MR = 14; 
 
 
   % Mass of gas in the regenerator (kg) 
 
 
   MH = 15; 
 
 
   % Mass of gas in the heater (kg) 
 
 
   ME = 16; 
 
 
   % Mass of gas in the expansion space (kg) 
 
 
   TCK = 17; 
 
 
   % Conditional temperature compression space / cooler (K) 
 
 
   THE = 18; 
 
 
   % Conditional temeprature heater / expansion space (K) 
 
 
   GACK = 19; 
 
 
   % Conditional mass flow compression space / cooler (kg/rad) 
 
 
   GAKR = 20; 
 
 
   % Conditional mass flow cooler / regenerator (kg/rad) 
 
 
   GARH = 21; 
 
 
   % Conditional mass flow regenerator / heater (kg/rad) 
 
 
   GAHE = 22; 
 
 
   % Conditional mass flow heater / expansion space (kg/rad) 
 
 
   % Size of var(ROWV,COL), dvar(ROWD,COL) 
 
 
   ROWV = 22; 
 
 
   % number of rows in the var matrix 
 
 
   ROWD = 16; 
 
 
   % number of rows in the dvar matrix 
 
 
   COL = 37; 
 
 
   % number of columns in the matrices (every 10 degrees) 
 
 
   %============================================================== 
 
 
   global 
 
 
   freq 
 
 
   % cycle frequency [herz] 
 
 
   global 
 
 
   tk tr th 
 
 
   % cooler, regenerator, heater temperatures [K] 
 
 
   global 
 
 
   cqwr 
 
 
   % regenerator housing thermal conductance [W/K] 
 
 
   twk = tk; 
 
 
   % Cooler wall temp - equal to initial cooler gas temp 
 
 
   twh = th; 
 
 
   % Heater wall temp - equal to initial heater gas temp 
 
 
   epsilon = 1; 
 
 
   % allowable temperature error bound for cyclic convergence 
 
 
   terror = 10*epsilon; 
 
 
   % Initial temperature error (to enter loop) 
 
 
   while 
 
 
   (terror>epsilon) 
 
 
   [var,dvar] = adiab; 
 
 
   qrloss = regsim(var); 
 
 
   % included 2/7/2010 
 
 
   tgh = hotsim(var,twh,qrloss); 
 
 
   % new heater gas temperature 
 
 
   tgk = kolsim(var,twk,qrloss); 
 
 
   % new cooler gas temperature 
 
 
   terror = abs(th - tgh) + abs(tk - tgk); 
 
 
   th = tgh; 
 
 
   tk = tgk; 
 
 
   tr = (th-tk)/log(th/tk); 
 
 
   end 
 
 
   fprintf( 
 
 
   '===== converged heater and cooler mean temperatures =====\n' 
 
 
   ); 
 
 
   fprintf( 
 
 
   'heater wall/gas temperatures: Twh = %.1f[K], Th = %.1f[K]\n' 
 
 
   ,twh,th); 
 
 
   fprintf( 
 
 
   'cooler wall/gas temperatures: Twk = %.1f[K], Tk = %.1f[K]\n' 
 
 
   ,twk,tk); 
 
 
   % Print out ideal adiabatic analysis results 
 
 
   eff = var(W,COL)/var(QH,COL); 
 
 
   % engine thermal efficency 
 
 
   Qkpower = var(QK,COL)*freq; 
 
 
   % Heat transferred to the cooler (W) 
 
 
   Qrpower = var(QR,COL)*freq; 
 
 
   % Heat transferred to the regenerator (W) 
 
 
   Qhpower = var(QH,COL)*freq; 
 
 
   % Heat transferred to the heater (W) 
 
 
   Wpower = var(W,COL)*freq; 
 
 
   % Total power output (W) 
 
 
   fprintf( 
 
 
   '========== ideal adiabatic analysis results ==========\n' 
 
 
   ); 
 
 
   fprintf( 
 
 
   ' Heat transferred to the cooler: %.2f[W]\n' 
 
 
   , Qkpower); 
 
 
   fprintf( 
 
 
   ' Net heat transferred to the regenerator: %.2f[W]\n' 
 
 
   , Qrpower); 
 
 
   fprintf( 
 
 
   ' Heat transferred to the heater: %.2f[W]\n' 
 
 
   , Qhpower); 
 
 
   fprintf( 
 
 
   ' Total power output: %.2f[W]\n', 
 
 
   Wpower); 
 
 
   fprintf( 
 
 
   ' Ideal Adiabatic Thermal efficiency: %.1f[%%]\n' 
 
 
   , eff*100); 
 
 
   fprintf( 
 
 
   '============= Regenerator analysis results============\n' 
 
 
   ); 
 
 
   fprintf( 
 
 
   ' Regenerator net enthalpy loss: %.1f[W]\n' 
 
 
   , qrloss*freq); 
 
 
   qwrl = cqwr*(twh - twk)/freq; 
 
 
   fprintf( 
 
 
   ' Regenerator wall heat leakage: %.1f[W]\n' 
 
 
   , qwrl*freq); 
 
 
   % Temperature plot of the simple simulation results 
 
 
   figure 
 
 
   x = 0:10:360; 
 
 
   Tcomp = var(TC,:); 
 
 
   Texp = var(TE,:); 
 
 
   plot(x,Tcomp, 
 
 
   'b-' 
 
 
   ,x,Texp, 
 
 
   'r-' 
 
 
   ); 
 
 
   hold on 
 
 
   x = [0,360]; 
 
 
   y = [twk,twk]; 
 
 
   plot(x,y, 
 
 
   'b--' 
 
 
   ) 
 
 
   y = [tk,tk]; 
 
 
   plot(x,y, 
 
 
   'b-' 
 
 
   ) 
 
 
   y = [tr,tr]; 
 
 
   plot(x,y, 
 
 
   'g-' 
 
 
   ) 
 
 
   y = [th,th]; 
 
 
   plot(x,y, 
 
 
   'r-' 
 
 
   ) 
 
 
   y = [twh,twh]; 
 
 
   plot(x,y, 
 
 
   'r--' 
 
 
   ) 
 
 
   hold off 
 
 
   grid on 
 
 
   xlabel( 
 
 
   'Crank angle (degrees)' 
 
 
   ); 
 
 
   ylabel( 
 
 
   'Temperature (K)' 
 
 
   ); 
 
 
   title( 
 
 
   'Simple Simulation - Wall and Gas Temps vs crank angle' 
 
 
   ); 
 
 
   % Various plots of the ideal adiabatic simulation results 
 
 
   plotadiab(var,dvar); 
 
 
   fprintf( 
 
 
   '========= pressure drop simple analysis =============\n' 
 
 
   ); 
 
 
   dwork = worksim(var,dvar); 
 
 
   fprintf( 
 
 
   ' Pressure drop available work loss: %.1f[W]\n' 
 
 
   , dwork*freq) 
 
 
   actWpower = Wpower - dwork*freq; 
 
 
   actQhpower = Qhpower + qrloss*freq + qwrl*freq; 
 
 
   acteff = actWpower/actQhpower; 
 
 
   fprintf( 
 
 
   ' Actual power from simple analysis: %.1f[W]\n' 
 
 
   , actWpower); 
 
 
   fprintf( 
 
 
   ' Actual heat power in from simple analysis: %.1f[W]\n' 
 
 
   , actQhpower); 
 
 
   fprintf( 
 
 
   ' Actual efficiency from simple analysis: %.1f[%%]\n' 
 
 
   , acteff*100);