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