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