dadiab

The ideal adiabatic model derivatives function 'dabiab'

We consider the numerical solution of the previously derived ideal adiabatic model equation set (refer to the Equation Summary and Method of Solution).

Note 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)

Our approach is to define two vectors, y and dy containing all the variables and derivatives. The first 7 variables are to be integrated numerically, and the remainder are evaluated algebraically in the derivative function dadiab shown below. This approach is called 'overloading' of the vectors, and the differential equation set is solved by using the Classical Fourth Order Runge-Kutta method as described in the Technical Note on ordinary differential equations. Notice how we have specified the indices of the vectors for improved readability.

function [y,dy] = dadiab(theta,y)
% Evaluate ideal adiabatic model derivatives
% Israel Urieli, 7/6/2002
% Arguments:  theta - current cycle angle [radians]
%             y(22) - vector of current variable values 
% Returned values: 
%             y(22) - updated vector of current variables
%             dy(16) vector of current derivatives
% Function invoked : volume.m
% global variables used from "define" functions
global vk % cooler void volume [m^3]
global vr % regen void volume [m^3]
global vh % heater void volume [m^3]
global rgas % gas constant [J/kg.K]
global cp % specific heat capacity at constant pressure [J/kg.K]
global cv % specific heat capacity at constant volume [J/kg.K]
global gama % ratio: cp/cv
global mgas % total mass of gas in engine [kg]
global tk tr th % cooler, regen, heater temperatures [K]
% Indices of the y, dy vectors:
 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)
%========================================================
% Volume and volume derivatives:
 [y(VC),y(VE),dy(VC),dy(VE)] = volume(theta);
% Pressure and pressure derivatives:
 vot = vk/tk + vr/tr + vh/th;
 y(P) = (mgas*rgas/(y(VC)/y(TC) + vot + y(VE)/y(TE)));
 top = -y(P)*(dy(VC)/y(TCK) + dy(VE)/y(THE));
 bottom = (y(VC)/(y(TCK)*gama) + vot + y(VE)/(y(THE)*gama));
 dy(P) = top/bottom;
% Mass accumulations and derivatives:
 y(MC) = y(P)*y(VC)/(rgas*y(TC));
 y(MK) = y(P)*vk/(rgas*tk);
 y(MR) = y(P)*vr/(rgas*tr);
 y(MH) = y(P)*vh/(rgas*th);
 y(ME) = y(P)*y(VE)/(rgas*y(TE));
 dy(MC) = (y(P)*dy(VC) + y(VC)*dy(P)/gama)/(rgas*y(TCK));
 dy(ME) = (y(P)*dy(VE) + y(VE)*dy(P)/gama)/(rgas*y(THE));
 dpop = dy(P)/y(P);
 dy(MK) = y(MK)*dpop;
 dy(MR) = y(MR)*dpop;
 dy(MH) = y(MH)*dpop;
% Mass flow between cells:
 y(GACK) = -dy(MC);
 y(GAKR) = y(GACK) - dy(MK);
 y(GAHE) = dy(ME);
 y(GARH) = y(GAHE) + dy(MH);
% Conditional temperatures between cells:
 y(TCK) = tk;
 if (y(GACK)>0)
       y(TCK) = y(TC);
 end
 y(THE) = y(TE);
 if (y(GAHE)>0)
       y(THE) = th;
 end
% 7 derivatives to be integrated by rk4:
% Working space temperatures:
 dy(TC) = y(TC)*(dpop + dy(VC)/y(VC) - dy(MC)/y(MC));
 dy(TE) = y(TE)*(dpop + dy(VE)/y(VE) - dy(ME)/y(ME));
% Energy:
 dy(QK) = vk*dy(P)*cv/rgas - cp*(y(TCK)*y(GACK) - tk*y(GAKR));
 dy(QR) = vr*dy(P)*cv/rgas - cp*(tk*y(GAKR) - th*y(GARH));
 dy(QH) = vh*dy(P)*cv/rgas - cp*(th*y(GARH) - y(THE)*y(GAHE));
 dy(WC) = y(P)*dy(VC);
 dy(WE) = y(P)*dy(VE);
% Net work done:
 dy(W) = dy(WC) + dy(WE);
 y(W) = y(WC) + y(WE);

The function volume is used to evaluate the working space volume variations and derivatives of the specific engine configuration as a function of the crank angle θ (theta). Recall from the engine function that we have included alpha engines as well as a free-piston beta drive engine in this resource, and it is intended that the function volume shown below be augmented as required.

function [vc,ve,dvc,dve] = volume(theta)
% determine working space volume variations and derivatives
% Israel Urieli, 7/6/2002
% Modified 2/14/2010 to include rockerV (rockvol)
% Modified 6/14/2016 to include beta free piston (sinevol)
% Argument:  theta - current cycle angle [radians]
% Returned values: 
%   vc, ve - compression, expansion space volumes [m^3]
%   dvc, dve - compression, expansion space volume derivatives 
global engine_type % s)inusoidal, y)oke r)ockerV (all alpha engines)

if (strncmp(engine_type,'s',1))
        [vc,ve,dvc,dve] = sinevol(theta);
elseif (strncmp(engine_type,'y',1))
        [vc,ve,dvc,dve] = yokevol(theta);
elseif (strncmp(engine_type,'r',1))
        [vc,ve,dvc,dve] = rockvol(theta);
elseif (strncmp(engine_type,'b',1))
        [vc,ve,dvc,dve] = sinevol(theta);
end
%========================================================
function [vc,ve,dvc,dve] = sinevol(theta)
% sinusoidal drive volume variations and derivatives
% Israel Urieli, 7/6/2002
% Argument:  theta - current cycle angle [radians]
% Returned values: 
%   vc, ve - compression, expansion space volumes [m^3]
%   dvc, dve - compression, expansion space volume derivatives 
global vclc vcle % compression,expansion clearence vols [m^3]
global vswc vswe % compression, expansion swept volumes [m^3]
global alpha % phase angle advance of expansion space [radians]
        
 vc = vclc + 0.5*vswc*(1 + cos(theta+pi));
 ve = vcle + 0.5*vswe*(1 + cos(theta + alpha+pi));
 dvc = -0.5*vswc*sin(theta+pi);
 dve = -0.5*vswe*sin(theta + alpha+pi);
%========================================================
function [vc,ve,dvc,dve] = yokevol(theta)
% Ross yoke drive volume variations and derivatives
% Israel Urieli, 7/6/2002
% Argument:  theta - current cycle angle [radians]
% Returned values: 
%   vc, ve - compression, expansion space volumes [m^3]
%   dvc, dve - compression, expansion space volume derivatives 
global vclc vcle % compression,expansion clearence vols [m^3]
global vswc vswe % compression, expansion swept volumes [m^3]
global alpha % phase angle advance of expansion space [radians]
global b1 % Ross yoke length (1/2 yoke base) [m]
global b2 % Ross yoke height [m]
global crank % crank radius [m]
global dcomp dexp % diameter of compression/expansion pistons [m]
global acomp aexp % area of compression/expansion pistons [m^2]
global ymin % minimum yoke vertical displacement [m]

 sinth = sin(theta);
 costh = cos(theta);
 bth = (b1^2 - (crank*costh)^2)^0.5;
 ye = crank*(sinth + (b2/b1)*costh) + bth;
 yc = crank*(sinth - (b2/b1)*costh) + bth;
 ve = vcle + aexp*(ye - ymin);
 vc = vclc + acomp*(yc - ymin);
 dvc = acomp*crank*(costh + (b2/b1)*sinth + crank*sinth*costh/bth);
 dve = aexp*crank*(costh - (b2/b1)*sinth + crank*sinth*costh/bth); 
%========================================================
function [vc,ve,dvc,dve] = rockvol(theta)
% Ross Rocker-V drive volume variations and derivatives
% Israel Urieli, 7/6/2002 & Martine Long 2/25/2005
% Argument:  theta - current cycle angle [radians]
% Returned values: 
%   vc, ve - compression, expansion space volumes [m^3]
%   dvc, dve - compression, expansion space volume derivatives 
global vclc vcle % compression,expansion clearence vols [m^3]
global crank % crank radius [m]
global acomp aexp % area of compression/expansion pistons [m^2]
global conrodc conrode % length of comp/exp piston connecting rods [m]
global ycmax yemax % maximum comp/exp piston vertical displacement [m]
        
 sinth = sin(theta);
 costh = cos(theta);
 beth = (conrode^2 - (crank*costh)^2)^0.5;
 bcth = (conrodc^2 - (crank*sinth)^2)^0.5;
 ye = beth - crank*sinth;
 yc = bcth + crank*costh;
 ve = vcle + aexp*(yemax - ye);
 vc = vclc + acomp*(ycmax - yc);
 dvc = acomp*crank*sinth*(crank*costh/bcth + 1);
 dve = -aexp*crank*costh*(crank*sinth/beth - 1);

______________________________________________________________________________________

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