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