The ideal adiabatic model function 'adiabatic' 

The purpose of the adiabatic function set is to determine the numerical solution of the ideal adiabatic model equation set (refer to the Equation Summary and Method of Solution). Recall in the equation set that apart from the constants, there are 22 variables and 16 derivatives to be solved over a complete cycle (θ = [0, 2π]) as follows:
Tc, Te, Qk, Qr, Qh, Wc, We - seven derivatives to be integrated numerically
W, p, Vc, Ve, mc, mk, mr, mh, me - nine analytical variables and derivatives
Tck, The, mck', mkr', mrh', mhe' - six conditional and mass flow variables (derivatives undefined)

Function adiabatic invokes the function set adiab to fill in the solution matrices var and dvar over a complete cycle, and then invokes function plotadiab (shown below) to do the required plots of the simulation results.

function [var,dvar] = adiabatic
% ideal adiabatic simulation and temperature/energy vs theta plots
% Israel Urieli, 7/20/2002
% 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 vk % cooler void volume [m^3]
global vr % regen void volume [m^3]
global vh % heater void volume [m^3]
% do ideal adiabatic analysis:
[var,dvar] = adiab;
% 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(' Thermal efficiency : %.1f[%%]\n', eff*100);
fprintf('========================================================\n')
% Various plots of the ideal adiabatic simulation results
plotadiab(var,dvar);


Function plotadiab is invoked both by function adiabatic (above) and by function simple.

function plotadiab(var,dvar)
% various plots of ideal adiabatic simulation results
% Israel Urieli, 7/21/2002 (corrected temp plots 12/3/2003) 
% Arguments: 
%   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 tk tr th % cooler, regenerator, heater temperatures [K]
global vk % cooler void volume [m^3]
global vr % regen void volume [m^3]
global vh % heater void volume [m^3]
choice = 'x';
while(~strncmp(choice,'q',1))
      fprintf('Choose plot type:\n');
      fprintf('   p - for a PV diagram\n');
      fprintf('   t - for a temperature vs crank angle plot\n');
      fprintf('   e - for an energy vs crank angle plot\n');
      fprintf('   q - to quit\n');
      choice = input('p)vdiagram, t)emperature, e)nergy, q)uit: ','s');
      if (strncmp(choice,'p',1))
            figure
            vol = (var(VC,:) + vk + vr + vh + var(VE,:))*1e6; % cubic centimeters
            pres = (var(P,:))*1e-5; % bar
            plot(vol,pres, 'k')
            grid on
            xlabel('Volume (cc)')
            ylabel('Pressure (bar [1bar = 100kPa])')
            title('P-v diagram')
      elseif (strncmp(choice,'t',1))
            figure
            x = 0:10:360;
            Tcomp = var(TC,:);
            Texp = var(TE,:);
            plot(x,Tcomp,'b-',x,Texp,'r-');
            hold on
            x = [0,360];
            y = [tk,tk];
            plot(x,y,'b-')
            y = [tr,tr];
            plot(x,y,'g-')
            y = [th,th];
            plot(x,y,'r-')
            hold off
            grid on
            xlabel('Crank angle (degrees)');
            ylabel('Temperature (K)');
            title('Temperature vs crank angle');
      elseif (strncmp(choice,'e',1))
            figure
            x = 0:10:360;
            Qkol = var(QK,:); % [J]
            Qreg = var(QR,:); % [J]
            Qhot = var(QH,:); % [J]
            Work = var(W,:); % [J]
            Wcom = var(WC,:); % [J]
            Wexp = var(WE,:); % [J]
            plot(x,Qkol,'b-',x,Qreg,'g-',x,Qhot,'r-',x,Work,'k-.',x,Wcom,'b--',x,Wexp,'r--');
            grid on
            xlabel('Crank angle (degrees)');
            ylabel('Energy [Joules]');
            title('Energy vs crank angle');
      end
end
fprintf('quitting ideal adiabatic plots...\n');



______________________________________________________________________________________


Stirling Cycle Machine Analysis by Israel Urieli is licensed under a Creative Commons Attribution-Noncommercial-Share Alike 3.0 United States License