The ideal adiabatic model function 'adiab'
The purpose of
function set adiab
is to fill in the arrays var
(variable values) and dvar
(derivative vaues) every 10 degrees over a cycle. The function set
includes the derivatives function dadiab,
the Classical fourth order Runge-Kutta function rk4,
and the function filmatrix
(both shown below).
Recall from the Equation
Summary and Method of Solution
that apart from the constants, there are 22 variables and 16
derivatives to be solved 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).
Notice how we have specified the indices
of the arrays for improved readability.
function [var,dvar] = adiab % ideal adiabatic model simulation % Israel Urieli, 7/6/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) global tk th % cooler, heater temperatures [K] % Row indices of the var, dvar matrices, and the y,dy variable 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) % Size of var(ROWV,COL), y(ROWV), dvar(ROWD,COL), dy(ROWD) 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) %=========================================================== fprintf('============Ideal Adiabatic Analysis====================\n') fprintf('Cooler Tk = %.1f[K], Heater Th = %.1f[K]\n', tk, th); epsilon = 1; % Allowable error in temperature (K) max_iteration = 20; % Maximum number of iterations to convergence ninc = 360; % number if integration increments (every degree) step = ninc/36; % for saving values in var, dvar matrices dtheta = 2.0*pi/ninc; % integration increment (radians) % Initial conditions: y(THE) = th; y(TCK) = tk; y(TE) = th; y(TC) = tk; iter = 0; terror = 10*epsilon; % Initial error to enter the loop % Iteration loop to cyclic convergence while ((terror >= epsilon)&(iter < max_iteration)) % cyclic initial conditions tc0 = y(TC); te0 = y(TE); theta = 0; y(QK) = 0; y(QR) = 0; y(QH) = 0; y(WC) = 0; y(WE) = 0; y(W) = 0; fprintf('iteration %d: Tc = %.1f[K], Te = %.1f[K]\n',iter,y(TC),y(TE)) for(i = 1:1:ninc) [theta,y,dy] = rk4('dadiab',7,theta,dtheta,y); end terror = abs(tc0 - y(TC)) + abs(te0 - y(TE)); iter = iter + 1; end if (iter >= max_iteration) fprintf('No convergence within %d iteration\n',max_iteration) end % Initial var and dvar matrix var = zeros(22,37); dvar = zeros(16,37); % a final cycle, to fill the var, dvar matrices theta=0; y(QK)=0; y(QR)=0; y(QH)=0; y(WC)=0; y(WE)=0; y(W)=0; [var,dvar] = filmatrix(1,y,dy,var,dvar); for(i = 2:1:COL) for(j = 1:1:step) [theta,y,dy] = rk4('dadiab',7,theta,dtheta,y); end [var,dvar] = filmatrix(i,y,dy,var,dvar); end |
The
function rk4
(below) was
developed in the Technical Note on ordinary
differential equations.
function [x, y, dy] = rk4(deriv,n,x,dx,y) %Classical fourth order Runge-Kutta method %Integrates n first order differential equations %dy(x,y) over interval x to x+dx %Israel Urieli - Jan 21, 2002 x0 = x; y0 = y; [y,dy1] = feval(deriv,x0,y); for i = 1:n y(i) = y0(i) + 0.5*dx*dy1(i); end xm = x0 + 0.5*dx; [y,dy2] = feval(deriv,xm,y); for i = 1:n y(i) = y0(i) + 0.5*dx*dy2(i); end [y,dy3] = feval(deriv,xm,y); for i = 1:n y(i) = y0(i) + dx*dy3(i); end x = x0 + dx; [y,dy] = feval(deriv,x,y); for i = 1:n dy(i) = (dy1(i) + 2*(dy2(i) + dy3(i)) + dy(i))/6; y(i) = y0(i) + dx*dy(i); end |
function [var,dvar]=filmatrix(j,y,dy,var,dvar); % Fill in the j-th column of the var, dvar matrices with values of y, dy % Israel Urieli, 7/20/2002 % Arguments: j - column index (1 - 37, every 10 degrees of cycle angle) % y(ROWV) - vector of current variable values % dy(ROWD) vector of current derivatives % var(ROWV,37) - matrix of current variables vs cycle angle % dvar(ROWD,37) - matrix of current derivatives vs cycle angle % Returned values: % var(ROWV,37) - matrix of updated variables vs cycle angle % dvar(ROWD,37) - matrix of updated derivatives vs cycle angle ROWV = 22; % number of rows in the var matrix ROWD = 16; % number of rows in the dvar matrix for (i = 1:1:ROWV) var(i,j) = y(i); end for (i = 1:1:ROWD) dvar(i,j) = dy(i); end |
______________________________________________________________________________________
Stirling Cycle Machine Analysis by Israel
Urieli is licensed under a Creative
Commons Attribution-Noncommercial-Share Alike 3.0 United States
License