Index Circular Stirling Stirling program Source code Links

Stirling engine design quickbasic program
The engine of the future project
The following is a translation of the isothermal second order calculation program, written by William R. Martini in Fortran.
I found it on:
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19830022057_1983022057.pdf
and translated it into quickbasic.
This version is very close to the original Fortran program.
How do you download the program? In fact it is already on your computer.
You can select the part between the ===== lines and copy it in a simple editor and save it with a name with the extention .bas
Download a version of Quick Basic or QBasic and load the program. The input is in the DATA lines on the end, and can be changed there.
I wrote a version that makes more use of the input and output features that modern programming languages offer.
This page is meant for progammers.
For users, these two versions are much better. Click here for the optimized version of this program.

For the html-javascript version of the Isothermal second order calculation Program (Written by William R. Martini in Fortran) Click here
The program is in the source code.It works in the browser, and also offline.

Click here for HOME

Problems

In this program is a line were
CX=P4*ID*LE*NE
Where CX is the dead volume in the connecting duct,
P4=pi/4, ID is the inside diameter of the connecting duct,
LE is the length of the connecting duct and NE is the number of connecting ducts per cylinder.
I am sure this must be an error.
It has to be:
CX=P4*ID*ID*LE*NE
As long as ID is almost one, the difference will not be much.
In the example inputs it is .76
With a much bigger or smaller diameter it gives a serious error.

The program works quite well with hydrogen, but the outputs for helium and air are not convincing.
For helium there is the line
M3=-9.3E8
I think it should be
M3=-9.3E-8
In the comment line is written:
"The properties of oxygen are used."
That has to be helium.
For air the program is not giving convincing outputs, but I still don't know how to fix it.
There the properties of oxygen could be used.


In the program are the lines
'AC = PI * IC! * LD! * NC * N
I corrected it to:
AC = PI * IC! * LD! * NC
' AC and AH were multiplied by N, but in all further
' calculations they are multiplied by N again.
PRINT "heat transfer area for cooler="; AC; "cm^2"
'AH = PI * IH! * LI! * NH * N
I corrected it to:
AH = PI * IH! * LI! * NH

N = number of cylinders per engine.
AC heat transfer area for cooler, cm^2
AH heat transfer area of heater, cm^2

The bugger factor BF to convert the power output to nearly what General Motors says it should be is .4 .
That factor differs too much from 1.
It means that without that correction the outputs differ a lot from reality.

On page 328 is written:
Note in editing: The program is valid for four cylinder engines only.
I think it is written specialy for the Siemens type Stirling engine, and checked for the GM4L23 engine.

The inputs for the heater tube outside diameter, OH and the cooler tube outside diameter, OC are not used in the program. They could have been used for conduction calculation.
'====================================================================
'Isothermal second order calculation
'Program 150 -10 oct 1979-
'Written by William R. Martini
'Program written with the Primos operating system
'program must have access to both the input file and an output file
'See attached reference for list and description of nomenclature
'**********************************************************************
'A  n/rm
'A1% counter for finding right average pressure
'AA .435 correlation of power with pressure
'AC heat transfer area for cooler, cm^2
'AF Area of flow, cm^2
'AH heat transfer area of heater, cm^2
'AL phase angle alpha = 90 degrees
'AS (was AS) area to volume ratio for regenerator matrix =
'   179 cm^/cm^2 for Met/Net 0.05-0.20 
'B  table spacing constant
'BA .1532 = exponent of correlation of power with pressure
'BF bugger factor to convert power outputs to nearly what GM says they
'   should be
'BH basic heat input, watts (BHI)
'BP basic power, watts
'c1 = QN * (100 / FE - 1)
'C() cold volumes at 360/nd points/cycle
'CD  cold dead volume, cm^3
'CF  cooler windage, watts
'CM  2.54 cm/inch            ??? not used
'CN  minimum FC()
'CP  heat capacity of hydrogen at constant p=14.62 j/g K @ 700 K (assumed
'    not to vary importantly with temperature)
'CR  length of connecting rod, cm
'CRT logical Unit no. for input file
'CV  heat capacity of hydrogen at constant volume=10.49 j/g K @ 700 k
'CW  friction factor for Met Net and others
'CX  cold dead volume outside cooler tubes, cm^3
'CY  maximum FC()
'd   overall eff.% = FE * NP! / QN
'DC  diameter engine cylinder, cm
'DD  diameter of piston drive rod, cm
'DN  360/nd
'DP  pressure drop, MPa
'DR  diameter of regenerator, cm
'DT  temperature rise in cooling water, K
'DU  temperature change for cold space, K
'DV  temperature change for hot space, K
'DW  diameter of "wire" in regenerator, cm=.0017(2.54)=0.00432 cm
'FuelInp   fuel input= 100 * QN / FE
'E1
'E2
'E3
'E4
'EC  piston end clearance, cm
'F   crank angle, degrees
'FC1 (F3+F4)/2
'FC() fraction of gas mass in cold spaces at 360/no points/cycle
'FE   furnace efficiency, %
'FF   filler factor, fraction of regenerator volume filled with solid
'FH1  (F1+F2)/2
'FH() fraction of gas mass in hot spaces at 360/no points/cycle
'FQ   60 Hz/rpm
'FR   (fh+fc)/2
'FW   flow of cooling water, g/sec
'FX   cooling water flow GPM @ 2000 rpm per cilinder
'F1   fraction of cycle time gas is assumed to leave hot space at constant
'     rate
'F2   fraction of cycle time gas is assumed to enter hot space at constant
'     rate
'F3   fraction of cycle time that flow out of cold space is assumed to occur
'     at constant rate
'F4   fraction of cycle time that flow into cold space is assumed to occur
'     at constant rate
'G    gap in hot cap, cm=0.56 cm
'GC   mass velocity through cooler, g/sec cm^2
'GD   mass velocity in connecting duct, g/sec cm^2
'GH   mass velocity in heater, g/sec cm^2
'GR   mass velocity in regenerator, g/sec cm^2
'H()  hot volumes at 360/nd points/cycle
'HC   heat transfer coefficient at cooler, w/cm^2 K
'HD   hot dead volume, cm^3
'HH   heat transfer coefficient in heater, w/cm^2 K
'HN   minimum FH()
'HP   1.341e-3 hp/watt
'HW = DP * WH * 2 * FH1 * A  ??? not in list
'HX   maximum FH()
'i    iteration counter
'IC!   ID! of cooler tube, cm

'ID!   inside diameter of connecting duct, cm
'IH!   InsDiam of heater tubes, cm
'IP!   indicated power, watts
'J    iteration counter
'KA!   coefficient for gas thermal conductivity calculation
'KB!   coefficient for gas thermal conductivity calculation
'KG!   gas thermal conductivity, watts/cm K
'KM!   metal thermal conductivity, w/cm K
'K3!   constant in reheat loss equation
'L1!   fraction of total gas charge leaking per MPa P per second
'L!()  gas inventory X gas constant, j/k (changes due to leak
'LB!   length of hot cap, cm
'LC!   length of cooler tube, cm
'LD!   heat transfer length of cooler tube, cm
'LE!   length of connecting duct, cm
'LH!   heater tube length, cm
'LI!   heater tube heat transfer length, cm
'LP   logical unit no. for output file
'LR!   length of regenerator, cm
'LX!   fraction of gas charge leaking per time increment per _P
'LY!   accumulation of MR!'s
'M!    number of moles of gas in working fluid, g mol
'ME!   mechanical efficiency, %
'MF!   mechanical friction loss
'MR!   gas inventory times gas constant, j/K
'MU!   gas viscosity, g/cm sec
'MW!   molecular weight, g/g mol
'MX!   mass of regenerator matrix
'M1! }
'M2! } }coefficient in viscosity equation
'M3! }
'N    number of cylinders per engine
'NC   number of cooler tubes per cylinder
'ND!   degree increment in time step (normally 30 degrees)
'NE   number of connecting ducts per cylinder
'NH   number of heat tubes per cylinder

'NP!   net power, watts
'NR   number of regenerators per cylinder
'NT!   number of tranfer units in regenerator, NTUB
'NU!   engine frequency, Hz
'n$   "Name"     ??? (can't find it in the program)
'OC   od of cooler tubes, cm
'OD   outside diameter of connecting duct, cm
'OG%   operating gas, 1=hydrogen, 2=helium, 3=air
'OH   heater tube OD, cm
'P()  pressure first with MR!=1, later at average pressure
'PG   average gas pressure, MPa
'PI   3.14159
'PM   mean pressure, of all p's
'PN   minimum pressure. MPs

'PP   0.006894 MPa/psia
'PR   Prandtl number to the 2/3 power = (Pr)^(2/3)
'PS   average pressure, psia
'PX   maximum pressure, MPa
'P4   pie/4=.785398
'QC   heat absorbed by cooler, watts
'QN   net heat required, watts
'QP   pumping loss for all N cylinders
'QS   shuttle loss, watts
'R    gas constant, 8.314 j/g mol K
'RA   0.0174533 radians/degree
'RC   crank radius, cm
'RD   regenerator dead volume, cm^3
'RE   Reynolds number, heater or cooler
'RH   reheat loss, watts
'RM   gas density for regenerator, g/cm^2
'RP   sum and average of power ratios
'RQ   sum and average of efficiency ratios
'RR   regenerator Reynolds number
'RT   Reynolds number, heater
'RW   regenerator windage, watts, for all cylinders in engine
'RZ   Reynolds number, cooler

'SC   wall thickness of hot cap, cm
'SE   wall thickness of expansion cylinder wall, cm
'SL   temp swing loss, watts=QTS
'SP   engine speed, RPM
'SR   wall thickness of regenerator housing, cm
'ST   Stanton number X(Pr)^(2/3)
'TC   effective cold space temperature, K
'TF   inside heater tube wall temperature, K
'TH   effective hot space temperature, K
'TM   inside heater tube wall temperature, K
'TR   regenerator temperature, K
'TS   matrix temp swing, K=DELTMX
'TW   inlet cooling water, K
'TX   cooler tube metal temperature, K
'TY   inlet cooling water temperature, F
'V()  total gas volume at 360/ND Points/cycle
'VC   velocity through gas cooler or connecting duct, cm/sec
'VH   velocity through gas heater, cm/sec
'VN   minimum total volume, cm^3
'VX   maximum total volume, cm^3
'v$   "value"
'WC   flow rate into or out of cold space, g/sec
'WH   flow rate into or out of hot space, g/sec
'WR   (WH+WC)/2=g/sec through regenerator= WRS
'W1   work for one cycle and one cylinder, joules
'X    temporary variable
'XX   correction factor to work diagram for large angle increments
'Y    temporary variable
'YY   temporary variable
'Z    temporary variable
'ZA%   0 for rapid iteration method, =1 for slower iteration method when
'     rapid method does not work
'ZB%   iteration counter
'ZH   specified static heat conduction loss, watts
'ZZ%   0 for specified static conduction, 1 for calculated static conduction

'init  ****************************************************************
'names with i,j,k,l,m and n are integers
DEFINT I-N
'sets up arrays (dimensions)
DIM H(13), C(13), P(13), FH(13), FC(13), V(14)
'sets up integers
'INTeger A1%, OG%, ZA%,ZB%,ZZ%,CRT,trm,N
'Sets up real numbers
'real IC!,ID!,IH!,IP!,KA!,KB!,KG!,KM!,K3!,L1!,LB!,LD!,LE!,LI!,LR!,LX!,LY!,M!,ME!,MF!,
'MR!,MU!,MW!,MX!,N1!,M2!,M3!,NP!,NU!,LC!,LH!,L!(14),NT!,ND!
DIM L!(14)
'Sets up logical unit numbers. 'CRT' is the logical unit number for
'the input file, and 'LP' is the logical unit for the output
'file
PI = 3.141593
P4 = PI / 4
cm = 2.54 'cm/inch   not used

'main
'Program reads in engine dimensions, operating conditions, and
'conversion constants from the input file. Also this is the return
'point after a case has been completed. If there are no more cases to
'run (i.e. an end of file occurs), the program calls exit.
300 'read(CRT,*,end=45) DC,LC!,LD!,IC!,OC,NC,PI
INPUT "new input"; I$
READ start$
WHILE start$ <> "bottom"

READ DC, LC!, LD!, IC!, OC, NC
'read(CRT,*) P4,DW,FX,ME!,FE,OG%,ZZ%
READ DW, FX, ME!, FE, OG%, ZZ%
'read(CRT,*) ZH,LH!,LI!,IH!,OH,NH,DD
READ ZH, LH!, LI!, IH!, OH, NH, DD
'read(CRT,*) RA,G,LB!,PS,KM!,SC,SE
READ RA, G, LB!, PS, KM!, SC, SE
'read(CRT,*)SR,LR!,DR,NR,FF,CR,RC
READ SR, LR!, DR, NR, FF, CR, RC
'read(CRT,*) N,AL,TF,TY,SP,AA,BA
READ N, AL, TF, TY, SP, AA, BA
'read(CRT,*) ID!,LE!,NE,BF,PP,CM,FQ
READ ID!, LE!, NE, BF, PP, cm, FQ
'read(CRT,*) R,HP,EC,L1!,AS
READ R, HP, EC, L1!, AS

PRINT DC; LC!; LD!; IC!; OC; NC; PI
PRINT P4; DW; FX; ME!; FE; OG%; ZZ%
'read(CRT,*) ZH,LH!,LI!,IH!,OH,NH,DD
PRINT ZH; LH!; LI!; IH!; OH; NH; DD
'read(CRT,*) RA,G,LB!,PS,KM!,SC,SE
PRINT RA; G; LB!; PS; KM!; SC; SE
'read(CRT,*)SR,LR!,DR,NR,FF,CR,RC
PRINT SR; LR!; DR; NR; FF; CR; RC
'read(CRT,*) N,AL,TF,TY,SP,AA,BA
PRINT N; AL; TF; TY; SP; AA; BA
'read(CRT,*) ID!,LE!,NE,BF,PP,CM,FQ
PRINT ID!; LE!; NE; BF; PP; cm; FQ
'read(CRT,*) R,HP,EC,L1!,AS
PRINT R; HP; EC; L1!'; AS
INPUT I$
'write(LP,4) DC,DR,IC!,OC,DW,DD,IH!,OH,G,LB!,LR!,CR,RC,LC!,LD!,LH!
PRINT "current dimensions are:"
PRINT "DC="; DC, "DR="; DR, "IC!="; IC!, "OC="; OC
PRINT "DW="; DW, "DD="; DD, "IH!="; IH!, "OH="; OH
PRINT "G="; G, "LB!="; LB!, "LR!="; LR!, "CR="; CR
PRINT "RC="; RC, "LC!="; LC!, "LD!="; LD!, "LH!="; LH!

'write(LP,5) LI!,NC,NR,N,NH,FF,AL,CX,NE,FE,EC,SC,SE,SR,ZZ%,ZH,KM!,ID!,
'LE!,NE,BF
PRINT "LI!="; LI!, "NC="; NC, "NR="; NR, "N="; N
PRINT "NH="; NH, "FF="; FF, "AL="; AL,
CX = P4 * ID! ^ 2 * LE! * NE
PRINT "CX="; CX
PRINT "ME!="; ME!, "FE="; FE, "EC="; EC, "SC="; SC
PRINT "SE="; SE, "SR="; SR, "ZZ%="; ZZ%, "ZH="; ZH
PRINT "KM!="; KM!, "ID!="; ID!, "LE!="; LE!, "NE="; NE
PRINT "BF="; BF
PRINT "------------------------------------------------"
PRINT "PI="; PI, "P4="; P4, "FX="; FX
INPUT I$

GOSUB loops
READ start$

WEND

IF start$ = "bottom" THEN
  PRINT "End of data."
  INPUT "End of program"; I$
  END
ELSEIF start$ <> ".null." THEN
  PRINT "First input must be .null."
  INPUT "Enter to continue"; I$
END IF
'End of main program
' call exit
END

loops:
'The degree increment is set at 30 degrees
ND! = 30
'A correction factor is calculated which increases the accuracy in
'calculating the work integrals with 30 degree increments.
XX = 1 + 5.321E-05 * ND! ^ 1.9797
'Temperature change for cold space (DU) and temperature change for hot
'space (DV) are set.
DU = 64
DV = 64
'The first thing the program does is to compute a list of engine
'volumes.

'Conversion to Kelvin degrees from input fahrenheit degrees.
TM = (TF + 460) / 1.8
PRINT "Heater tube inside temp.="; TF; "fahrenheit="; TM; "Kelvin="; TM - 273.15; "Celcius"
TW = (TY + 460) / 1.8
PRINT "Cooling water temperature="; TY; "fahrenheit="; TW; "Kelvin="; TW - 273.15; "Celcius"
'conversion to hertz and to npa.
NU! = SP / 60
PRINT "engine speed="; SP; " RPM="; NU!; "RPS"
PG = .006894 * PS
PRINT "average pressure="; PS; " psia="; PG; "MPa"

'Determines gas property values from 'OG%' (
'if 'OG%'= 1, the property values for hydrogen are used.
'If 'OG%' = 2, the property values for oxygen are used.
'If 'OG%' = 3. the property values for air are used.)
'Property values for additional gases may be added if desired'
'IF OG% = 1 THEN GOTO 20
'IF OG% = 2 THEN GOTO 21
IF OG% = 3 THEN 'air property values
  PRINT "WORKING GAS= AIR"
  KA! = -12.6824
  KB! = .782
  CP = 1.0752
  CV = .7883
  M1! = 1.8194E-04
  M2! = 5.36E-07
  M3! = 1.22E-06
  MW! = 29
  PR = .9071
ELSEIF OG% = 1 THEN 'hydrogen property values
  PRINT "WORKING GAS= HYDROGEN"
  KA! = -11.0004
  KB! = .813
  CP = 14.62
  CV = 10.49
  M1! = 8.873E-05
  M2! = .0000002
  M3! = 1.18E-07
  MW! = 2.02
  PR = .8408
ELSEIF OG% = 2 THEN 'oxygen property values
  PRINT "WORKING GAS= OXYGEN"
  KA! = -10.1309
  KB! = .6335
  CP = 5.2
  CV = 3.12
  M1! = 1.6614E-04
  M2! = 4.63E-07
  M3! = -9.3E+08
  MW! = 4
  PR = .8018
ELSE
  PRINT "  _______ #### ERROR, INPUT FAILURE #### ___________"
  PRINT "OG% = wrong, must be 1, 2 or 3 for hydrogen, oxygen, or air."
  INPUT "ENTER for end program"; I$
END IF
'Conversion of cooling water flow to grams/second. Initially cooler
'tube metal temperature is made the same as the inlet cooling water
'temperature. The total heat transfer areas for all the engine's
'coolers and all the engine's heaters are calculated.

FW = 63.12 * FX
' GPM @ 2000 rpm per cylinder
PRINT "flow of cooling water="; FX; "GPM at 2000 rpm per cylinder="; FW; "g/sec"
TX = TW
'AC = PI * IC! * LD! * NC * N
AC = PI * IC! * LD! * NC
'This must be an error.
' AC and AH were multiplied by N, but in all further
' calculations they are multiplied by N again.
PRINT "heat transfer area for cooler="; AC; "cm^2"
'AH = PI * IH! * LI! * NH * N
AH = PI * IH! * LI! * NH
PRINT "heat transfer area of heater="; AH; "cm^2"

'calculates engine dead volumes and initializes pressures and volumes.
'Initializes for determination of average pressure and maximum and
'minimum volumes.

HD = P4 * IH! * IH! * LH! * NH + EC * DC ^ 2 * P4
PRINT
PRINT P4 * DC ^ 2 * EC; "= piston end clearance, dead volume, cm^3"
PRINT P4 * IH! ^ 2 * LH! * NH; "= heater tube dead volume, cm^3"
PRINT "+ -----------"
PRINT HD; "= total hot dead volume, cm^3"
PRINT
'CX = P4 * ID! * LE! * NE
'          ID!   inside diameter of connecting duct, cm
'                LE!   length of connecting duct, cm
'                      NE   number of connecting ducts per cylinder

CX = P4 * ID! * LE! * NE '  ^2 was forgotten in original
'CX = P4 * ID! ^ 2 * LE! * NE'  ^2 was forgotten in original

RD = (1 - FF) * P4 * DR ^ 2 * LR! * NR + PI * DC * G * LB!
PRINT (1 - FF) * P4 * DR ^ 2 * LR! * NR; "= regenerator inside dead volume, cm^3"

PRINT PI * DC * G * LB!; "= gap in hot cap dead volume, cm^3"
                        'pi* DiameterEngineCylinder*GapInHotCap*LengthHotCap

PRINT "+ -----------"
PRINT RD; "= total regenerator dead volume, cm^3"
PRINT
CD = CX + P4 * IC! ^ 2 * LC! * NC + EC * P4 * (DC ^ 2 - DD ^ 2)
PRINT CX; "= cold dead volume outside cooler tubes, in connecting duct cm^3"
PRINT P4 * IC! ^ 2 * LC! * NC; "= cooler tube dead volume cm^3"
PRINT EC * P4 * (DC ^ 2 - DD ^ 2); "= end clearance cold dead volume cm^3"
'EC  piston end clearance, cm
'DC  diameter engine cylinder, cm
'DD  diameter of piston drive rod, cm
PRINT "+ -----------"
PRINT CD; "= total cold dead volume, cm^3 (incl. connecting duct)"

PN = 0
VX = 0
VN = 1E+30
'initially sets the effective hot space temperature to the hot metal
'temperature and the effective cold space temperature to the cooling
'metal temperature for the first time around. Calculates the log mean
'temperature for the regenerator. Calculates the leakage coefficient
'for 30 degree increments.
TH = TM
TC = TW
TR = (TM - TW) / LOG(TM / TW)
LX! = L1! * ND! / (360 * NU!)
'Since the thermoconductivity enters the calculation only at the
'regenerator temperature it can be calculated before the main
'iteration loop
KG! = EXP(KA! + KB! * LOG(TR))
'Start of do loop 23 to calculate engine volumes.

FOR I = 1 TO 13
'Calculates the hot volume and cold volume for each angle increment for
'crank operated pistons. Since a double acting machine has a piston
'drive rod (DD) and a single acting machine does not, 'DD' is used as
'an indicator of whether the cold volume of the engine is above the
'piston or below it.
  X = 30 * (I - 1) * RA
  J = I
  IF DD <> 0 THEN
    Y = (30 * (I - 1) + AL) * RA
  ELSE
    Y = (30 * (I - 1) - AL) * RA
  END IF
  H(J) = P4 * DC ^ 2 * (RC - SQR(CR ^ 2 - (RC * SIN(X)) ^ 2) + RC * COS(X) + CR) + HD

  IF DD <> 0 THEN
    C(J) = P4 * (DC ^ 2 - DD ^ 2) * (SQR(CR ^ 2 - (RC * SIN(Y)) ^ 2) - RC * COS(Y) - CR + RC) + CD
  ELSE
    C(J) = P4 * DC ^ 2 * (RC - SQR(CR ^ 2 - (RC * SIN(Y)) ^ 2) + RC * COS(Y) + CR) + CD
  END IF
  'calculates the total gas volume and finds the maximum volume.
  V(J) = H(J) + RD + C(J)
  IF V(J) > VX THEN VX = V(J)
  'finds the minimum volume.
  IF V(J) < VN THEN VN = V(J)
  'calculates the initial gas inventory.
  IF J = 3 THEN L!(1) = PG * (H(J) / TH + RD / TR + C(J) / TC)
  'end of loop to calculate engine volumes
NEXT I'continue
' 'ZA%' is set at zero so that the fastest way of arriving at the proper
'effective hot space and cold space temperature will be tried first.
'also a counter, 'ZB%', is set at zero.
ZA% = 0
ZB% = 0
GOSUB 200
RETURN

'initialization
200 A = 0
WHILE (X > .0001) OR (Z > .01)

 PM = 0
LY! = 0
'start of loop (to calculate pressures for one engine cycle).
FOR I = 1 TO 13
  'calculate pressure
  P(I) = L!(I) / (H(I) / TH + RD / TR + C(I) / TC)
  'Calculate gas inventory for next increment due to leakage
  L!(I + 1) = L!(I) * (1 - LX! * (P(I) - PG))
  'Accumulate values, mean pressure and mean gas inventory.
  IF I <> 1 THEN
    PM = PM + P(I)
    LY! = LY! + L!(I)
  END IF
  'end of loop (to calculate pressures for one engine cycle)
NEXT I
'indexes cycle counter, calculates mean pressure, readjusts gas
'inventory taking into account gas leakage.
A = A + 1
PM = PM / 12
IF A >= 3 THEN
  L!(1) = L!(13)
ELSE
  L!(1) = L!(13) * PG / PM
END IF
'convergence criteria: pressure from beginning to the end of cycle
'must not change by more than one hundredth of a percent and the mean
'pressure must be within one percent of the desired gas pressure.
'Usually one or two cycles are required to meet this criteria.
X = ABS(P(1) - P(13))
PRINT X; "Pressure difference after loop"; A
Z = ABS(PM - PG)
PRINT Z; "Volume difference after loop"; A
WEND

PRINT "Pressure and Volume are O.K."
INPUT I$
'initializing
W1 = 0
PX = P(1)' was 0, this is better
PN = P(1)' was 10000, this is better
MR! = LY! * ND! / 360
'Start of loop (finds the maximum and minimum pressure).
FOR I = 1 TO 13
  IF P(I) > PX THEN PX = P(I)
  IF P(I) < PN THEN PN = P(I)
NEXT I'continue
PRINT PX; "= Maximum pressure"
PRINT PN; "= Minimum pressure"
INPUT I$
'Start of loop (finds the work per cycle by integrating the
'pressure volume loop).
FOR I = 1 TO 12
  W1 = W1 + (P(I) + P(I + 1)) * (V(I + 1) - V(I)) / 2
NEXT I'continue
'Basic power for the whole engine is calculated from the integrated
'power using the correction factor XX which compensates for the
'truncation error of using only a small number of points to integrate.
BP = NU! * XX * W1 * N
'Initializing
HX = 0
CY = 0
HN = 1
CN = 1
'Calculates an array giving the fraction of the total gas inventory in
'the hot space and in the cold space for each point during the cycle.
'INPUT "ok?"; i$
FOR I = 1 TO 13
  FH(I) = P(I) * H(I) / (MR! * TH)
  IF FH(I) > HX THEN HX = FH(I)
  IF FH(I) < HN THEN HN = FH(I)
  FC(I) = P(I) * C(I) / (MR! * TC)
  IF FC(I) > CY THEN CY = FC(I)
  IF FC(I) < CN THEN CN = FC(I)
NEXT I' continue
'If FH(i) and FC(i) are graphed as a function of the angle, it is seen
'that a good approximation of the graph is to have two periods per
'cycle of constant mass flow interspersed with periods of no flow at
'all. F1 to F4 are the fractions of the total cycle time when
'different flows are assumed to occur (see nomenclature).
'When 'FH1' and 'FC1' are calculated, the average cycle time, when flow
'is assumed to occur either into or out of the hot space and either
'into or out of the cold space, is calculated.
F1 = (HX - HN) / (6 * (FH(1) - FH(3)))
F2 = (HX - HN) / (6 * (FH(10) - FH(8)))
F3 = (CY - CN) / (6 * (FC(8) - FC(10)))
F4 = (CY - CN) / (6 * (FC(3) - FC(1)))
FH1 = (F1 + F2) / 2
FC1 = (F3 + F4) / 2
'Effective mass flow into or out of the hot space is calculated.
M! = MR! / R
WH = (HX - HN) * M! * MW! * NU! / FH1
'Effective mass flow into or out of the cold space is calculated.
WC = (CY - CN) * M! * MW! * NU! / FC1
'Fraction of the time the flow is assumed to pass through the
'regenerator and the flow rate of the regenerator is calculated as the
'average between the hot and cold flows.
FR = (FH1 + FC1) / 2
WR = (WH + WC) / 2
'Regenerator gas density.
RM = .1202 * MW! * PG / TR
'Calculates regenerator windage loss.
MU! = M1! + M2! * (TR - 293) + M3! * PG
GR = WR / (P4 * DR ^ 2 * NR)
RR = DW * GR / MU!
CW = 2.7312 * (1 + 10.397 / RR)
DP = CW * GR ^ 2 * LR! / (2E+07 * DW * RM)
A = N / RM
RW = DP * WR * 2 * FR * A
'Calculates heater windage loss. In this calculation the viscosity for
'the input temperature and subroutine 'rest' returns the friction
'factor for the input Reynolds number. The calculation takes into
'account friction losses, as well as 4.4 velocity heads for an
'entrance and an exit loss, one 180 degree bend, and two 90 degree
'bends.
MU! = M1! + M2! * (TM - 293) + M3! * PG
RM = .1202 * MW! * PG / TM
A = N / RM
GH = WH / (P4 * IH! ^ 2 * NH)
RE = IH! * GH / MU!
RT = RE

IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF

AF = P4 * IH! ^ 2 * NH
VH = WH / (RM * AF)
DP = 2 * CW * GH ^ 2 * LH! / (1E+07 * IH! * RM) + VH ^ 2 * 4.4 * RM / 2E+07
HW = DP * WH * 2 * FH1 * A
'This calculates the windage loss through the gas cooler and the
'connecting tube. The same comments for the gas heating windage loss
'apply here as well. The velocity heads charge to the gas cooler is
'1.5 for a simple entrance and exit loss. In the connecting head line,
'three velocity heads are changed to account for entrance and exit loss
'plus two 90 degree bends.
'NU! = M1! + M2! * (TX - 293) + M3! * PG ??? NU is input
MU! = M1! + M2! * (TX - 293) + M3! * PG
'MU!   gas viscosity, g/cm sec
RM = .1202 * MW! * PG / TX
A = N / RM
GC = WC / (P4 * IC! ^ 2 * NC)

RE = IC! * GC / MU!
RZ = RE
IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF
AF = P4 * IC! ^ 2 * NC
VC = WC / (RM * AF)
DP = 2 * CW * GC ^ 2 * LC! / (1E+07 * IC! * RM) + VC ^ 2 * 1.5 * RM / 2E+07
GD = WC / (P4 * ID! ^ 2 * NE)
RE = ID! * GD / MU!

IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF

AF = P4 * ID! ^ 2 * NE
VC = WC / (RM * AF)
DP = DP + 2 * CW * GD ^ 2 * LE! / (1E+07 * ID! * RM) + VC ^ 2 * 3 * RM / 2E+07
CF = DP * WC * 2 * FC1 * A
'Calculates indicated power.
IP! = BP - HW - RW - CF

'Calculates mechanical friction loss.
MF! = (1 - ME! / 100) * IP!
'Calculates net power.
NP! = IP! - MF!
'Calculates basic heat input.
BH = BP / (1 - TC / TH)
'Calculates reheat loss for net net .05-.20 which is used in the 4l23
'machine. This section is specific for this type of regenerator
'material.

IF RR < 42 THEN
  X = EXP(-.1826 - .0583 * LOG(RR))
ELSEIF RR < 140 THEN
  X = EXP(.5078 - .243 * LOG(RR))
ELSE
  X = EXP(1.78 - .5044 * LOG(RR))
END IF

NT! = X * LR! / DW
X = WR * CP * (TH - TW)
Y = RD * CV * (PX - PN) * NU! * MW! / (R * FR)
K3! = FR * (X - Y)
RH = K3! / (NT! + 2) * N * 2
'Calculates temperature swing loss
MX! = NR * P4 * DR ^ 2 * LR! * FF * 7.5
TS = K3! / (NU! * MX! * 1.05)
SL = K3! * TS * N / (2 * (TH - TX))
'Calculates pumping or appendix loss.
X = (PI * DC / KG!) ^ .6
Y = ((PX - PN) * MW! * NU! * CP * 2 / ((TH + TX) * R)) ^ 1.6
Z = G ^ 2.6
QP = N * X * 2 * LB! * (TH - TX) * Y * Z / 1.5
'Calculates shuttle heat loss.
QS = 2 * P4 * RC * RC * KG! * (TH - TC) * DC / (G * LB!) * N
'Calculates static heat loss. This can be either specified or
'calculated from the basic dimensions.
IF ZZ% = 1 THEN ZH = (TH - TC) * (KM! * ((DR ^ 2 * P4 * FF + PI * DR * SR) / LR! + PI * DC * (SC + SE) / LB!) + KG! * (DR ^ 2 * P4 * (1 - FF) / LR! + DC ^ 2 * P4 / LB!))

'Sums all losses to calculate net heat demand.

QN = BH + ZH + SL + RH - HW - RW / 2 + QS + QP
'Calculates cooler heat load.
QC = QN - NP!
'Temperature rise in cooling water.
DT = QC / (FW * 4.185)
'Effective cold metal temperature.
TX = TW + DT / 2
'Calculates heat transfer coefficient in the cold heat exchanger.
RE = RZ
J = 1
'GOTO subroutine rest
GOSUB 100
RETURN

47  DV = 64
    DU = 64
    ZA% = 1
    TC = TW
    TH = TM
    GOSUB 46
RETURN


'On the first time through 'TC' = 'TW' and the error in the cold space,
'E2, is made equal to the required heat transfer through the gas
'coolers, 'QC'. Then the next estimate for 'TC' is made by adding
' 'DU', 64 degrees, to 'TX', the average temperature of the gas
'cooler metal. The program then goes to 48, skipping over the rest of
'the ad+fusiment program for the cold space.
46 IF TC = TW THEN
  'Location of control if 'TC' equals 'TW'.
49 E2 = QC
  TC = TX + DU
  GOSUB 48
  RETURN
END IF

'if 'TC' is not equal to 'TW', as it will be for anything except
'for the first time through, the previous error is saved as 'E1'.
'than 'E2' is calculated as the difference between the heat that
'should be transferred and the heat that can be transferred by the
'capabilities of the heat exchanger.
E1 = E2
E2 = QC - HC * FC1 * AC * N * (TC - TX) * BF
'if this error is positive, then the correction number, 'DU', is
'added to the cold temperature, 'TC', and the program goes on to the
'hot space analysis.
IF E2 > 0 THEN
' Location of control if 'E2' is greater than 0.
  TC = TC + DU
  GOSUB 48
  RETURN
END IF

'If this error is negative and the previous error was positive,
'then the degree increment 'DU', is just divided by 4, for future
'corrections.
IF (E2 < 0) AND (E1 > 0) THEN DU = DU / 4
'The degree increment is subtracted from 'TC'. If 'TC' becomes
'greater than 'TM', the hot metal temperature, obviously there is
'insufficient cooler heat transfer area and the program stops for
'this case. This can occur for small cooler areas and specified heat
'leaks.
TC = TC - DU
PRINT "TC="; TC
PRINT "DU="; DU
PRINT "TM="; TM
PRINT "tc-tm="; TC - TM
INPUT I$
IF TC > TM THEN
 'Because of insufficient cooler area the program is terminated for
 'this case.
 'write(LP+1)
 '1 format(10('*')'Insufficient cooler area',10('*'))
  PRINT "Insufficient cooler area"
  RETURN 'GOTO 300
END IF
GOSUB 48
RETURN

'Calculates heat transfer coefficient for gas heater. Flag 'ZA%'
'indicates whether the fast method of convergence at 59 or the slow
'method at 52 should be used.
48 RE = RT
J = 2
'GOTO subroutine rest
GOSUB 100
RETURN


100 'Subroutine rest
'Calculates Stanton number from Reynolds number
IF RE >= 10000 THEN
  ST = EXP(-3.57024 - .229496 * LOG(RE))
ELSEIF RE > 7000 THEN
  ST = .0034
ELSEIF RE > 4000 THEN
  ST = EXP(-13.3071 + .861016 * LOG(RE))
ELSEIF RE > 3000 THEN
  ST = .0021
ELSE
  ST = EXP(.337046 - .812212 * LOG(RE))
END IF

IF J = 1 THEN
   HC = ST * CP * GC / PR

 'Two different methods of arriving at the proper effective hot space
 'and cold space temperature are interspersed. The fastest way,
 'which is usually tried first, involves calculating what the
 'temperature difference has to be between the metal temperature and
 'the effective gas temperature considering the heat transfer
 'capability of the heat exchanger and the correction factor.
 'However, if the heat exchanger is too small, the first iteration
 'method goes unstable and a second, more cautious, method must be
 'employed. The 'ZA%' is the flag which shows that the second
 'method is called in.
  IF ZA% = 1 THEN
    GOSUB 46
    RETURN
  END IF
 ' 'X' is used as a temporary variable for the previous cold
 'temperature. The cold temperature is calculated, assuming there is
 'no error between the heat that can be transferred and the heat that
 'should be transferred. Counter 'ZB%' is indexed. A test is now made
 'of the 'TC' value just calculated. If the effective cold gas
 'temperature is greater than the effective hot gas temperature or
 'less than the cooling water temperature this iteration method has
 'gone unstable and the second, more cautious, method is called in.
 'Also if the first iteration method has not come to an answer within
 '10 iterations, ('ZB%' greater than 10), the second iteration method
 'is called in. The initial change in the hot gas temperature, 'DV',
 'and in the cold gas temperature. 'DU', are both set at 64 degrees.
 'The flag 'ZA%' is set at 1 and 'TC' and 'TH' are set at the initial

 'values. Control pass to 46 where the second approach begins. If
 'the value of 'TC' does not indicate the second approach is needed
 'control passes to 48 to start calculation of the effective
 'temperature in the hot space.
  X = TC
  YY = HC * FC1 * AC * N * BF
  TC = QC / YY + TX
  E2 = QC - YY * (TC - TX)
  ZB% = ZB% + 1
  IF (TC > TH) OR (TC < TX) OR (ZB% > 10) THEN
    GOSUB 47
    RETURN
  END IF
  GOSUB 48
  RETURN
 ELSE
    HH = ST * CP * GH / PR           'ST of RT @@@
   'HH   heat transfer coefficient in heater, w/cm^2 K
   'ST   Stanton number X(Pr)^(2/3)
'RT   Reynolds number, heater
'CP  heat capacity of hydrogen at constant p=14.62 j/g K @ 700 K (assumed
'    not to vary importantly with temperature)
'GH   mass velocity in heater, g/sec cm^2
'PR   Prandtl number to the 2/3 power = (Pr)^(2/3)

  IF ZA% = 1 THEN

   'This is analogous to 46 to 48, except this is for the hot space.
    IF TH = TM THEN
     'Location of control if 'TH' equals 'TH'.
      E4 = QN
      TH = TM - DV
    ELSE

      E3 = E4
      E4 = QN - HH * FH1 * AH * N*(TM - TH) * BF
      IF E4 <= 0 THEN
        IF (E4 < 0) AND (E3 > 0) THEN DV = DV / 4
        TH = TH + DV
        IF TH < TW THEN
         'Because of insufficient heater area the program is terminated for
         'this case
          PRINT "Insufficient heater area"
          RETURN 'GOTO 300
        END IF
      END IF

    ' Location of control if 'TH' is not less than 'TM'.
      TH = TH - DV
    END IF

'   Convergence criteria for the slower method of iteration.
'   Convergence is complete when the air in the hot space and the air in
'   the cold space are both less than 1% of the heat transferred through

'   the heat exchangers.
    x1 = ABS(E4)
    x2 = ABS(E2)
    x3 = QN / 100
    x4 = QC / 100
    IF (x1 > x3) OR (x2 > x4) THEN
      GOSUB 200
      RETURN
    END IF
 '  Location of control for the second iteration method.
    GOSUB 47
    RETURN
  END IF

'This is analogous to the comment made after 44 on the cold space,
'except this is for the hot space.
Y = TH
YY = HH * FH1 * AH * N * BF
TH = TM - QN / YY
E4 = QN - YY * (TM - TH)
IF (TH > TM) OR (TH < TC) THEN
  GOSUB 47
  RETURN
END IF
 'Convergence criteria for the first iteration method. The iteration
'is complete when change in the effective hot space and cold space
'temperature is less than one degree kelvin per iteration.
x1 = ABS(TH - Y)
x2 = ABS(TC - X)
IF (x1 > 1) OR (x2 > 1) THEN
  GOSUB 200
  RETURN
END IF

'completes preparation for output
 A = -HW - RW / 2
B = 100 * IP! / QN
C1 = QN * (100 / FE - 1)
D = FE * NP! / QN
FuelInp = 100 * QN / FE
'Reinitializing
I = I + 1
ZA% = 0
ZB% = 0

'This is where the printing of the output starts. To compress output
'the operating conditions and engine dimensions are identified only by
'their Fortran symbol.

'Prints program heading
'write(LP,10)
'Prints current operating conditions
'write (LP+3) SP,PS,ND!,TF,L1!,TY,FX,OG%
PRINT "Current operating conditions are:"
PRINT "SP="; SP, "PS="; PS, "ND!="; ND!, "TF="; TF
PRINT "L1!="; L1!, "TY="; TY, "FX="; FX, "OG%="; OG%
PRINT
'Prints current dimensions
'write(LP,4) DC,DR,IC!,OC,DW,DD,IH!,OH,G,LB!,LR!,CR,RC,LC!,LD!,LH!
PRINT "current dimensions are:"
PRINT "DC="; DC, "DR="; DR, "IC!="; IC!, "OC="; OC
PRINT "DW="; DW, "DD="; DD, "IH!="; IH!, "OH="; OH
PRINT "G="; G, "LB!="; LB!, "LR!="; LR!, "CR="; CR
PRINT "RC="; RC, "LC!="; LC!, "LD!="; LD!, "LH!="; LH!

'write(LP,5) LI!,NC,NR,N,NH,FF,AL,CX,ME!,FE,EC,SC,SE,SR,ZZ%,ZH,KM!,ID!,
'LE!,NE,BF
PRINT "LI!="; LI!, "NC="; NC, "NR="; NR, "N="; N
PRINT "NH="; NH, "FF="; FF, "AL="; AL, "CX="; CX
PRINT "ME!="; ME!, "FE="; FE, "EC="; EC, "SC="; SC
PRINT "SE="; SE, "SR="; SR, "ZZ%="; ZZ%, "ZH="; ZH
PRINT "KM!="; KM!, "ID!="; ID!, "LE!="; LE!, "NE="; NE
PRINT "BF="; BF
PRINT

'Prints power outputs and heat inputs
'write(LP,6) BP,BH,NW,RH,RW,QS,CF,QP,IP!,SL,MF!,ZH,NP!,A
PRINT "power, watts", , "heat requirement, watts"
PRINT "basic", BP, "basic", BH
PRINT "heater f.l.", HW, "reheat", RH
PRINT "regen. f.l", RW, "shuttle", QS
PRINT "cooler f.l.", CF, "pumping", QP
PRINT "net       ", IP!, "temp,swing", SL
PRINT "mech.fric.", MF!, "conduction", ZH
PRINT "brake", NP!, "flow fric.credit"; A
PRINT "---------------------------",
'write(LP,7) QN,B,C1,D,FuelInp
PRINT "heat to engine"; QN
PRINT "indicated eff.%="; B, "furnace loss", C1
PRINT "overall   eff.%="; D, "fuel input", FuelInp
'write(LP,8) TM,TW,TH,TC
PRINT "----------------------------------------------------------------------"
PRINT "hot metal temp.   k="; TM; "  cooling water inlet; temp., k ="; TW
PRINT "effec.hot sp.temp.k="; TH; "  effec.cold sp.temp.k.="; TC
PRINT "----------------------------------------------------------------------"
INPUT I$
'Prints work diagram from data
'write(LP,9)
PRINT "final work diagram:"
'PRINT "angle", "hot vol", "cold vol.", "tot. vol,", "pressure"; "  gas inv."
PRINT "angle";
PRINT RIGHT$("             " + "HOT VOL.", 13);
PRINT RIGHT$("             " + "COLD VOL.", 13);
PRINT RIGHT$("             " + "TOT. VOL.", 13);
PRINT RIGHT$("             " + "PRESSURE", 13);
PRINT RIGHT$("             " + "GAS INV.", 13)

FOR I = 1 TO 13
  F = ND! * I - 30
  G = L!(I) / R
  'write(LP,11) F,H(i),C(i),V(i),P(i),G
'  PRINT F; H(I), C(I); "  "; V(I), p(I); "  "; G
  PRINT RIGHT$("     " + STR$(F), 5);
  PRINT RIGHT$("             " + STR$(H(I)), 13);
  PRINT RIGHT$("             " + STR$(C(I)), 13);
  PRINT RIGHT$("             " + STR$(V(I)), 13);
  PRINT RIGHT$("             " + STR$(P(I)), 13);
  PRINT RIGHT$("             " + STR$(G), 13)
'  PRINT H(I), C(I); "  "; V(I), P(I); "  "; G
NEXT I
INPUT I$

RETURN 'GOTO 300
  
END IF

'output fotmat:
'1 format(10('*')'Insufficient cooler area',10('*'))
'2 format(10('*')'Insufficient heater area',10('*'))
'3 format('Current operating conditions are:'/'SP=',f10.2,t17,'PS=',
'  1f10,2,t33 "ND!=" TF+
'4 format("current dimensions are:"
'5 format(
'6 format("power, watts","heat requirement, watts","basic","basic",
'  "heater","reheat","regen.f.l","shuttle" ,"cooler f.l.","pumping","net"
'  "temp,swing","mech.fric.","conduction","brake","flow fric.","CR","edit"
'7 format("heat to engine","indicated eff.%=";"furnace loss";"overall eff.%="
'  "fuel input"
'8 format("hot metal temp. k=","cooling water ","inlet temp., k="
'  "effec.hot sp.temp.k=","effec.cold sp.temp.k.="
'9 format("final work diagram:","angle","hot vol","cold vol."
'  "tot. vol,","pressure","gas inv."
'10 format("isothermal second order calculation--", "prog, 150"
'  "10 oct 1979","written by William R. Martini"

DATA ".null."
DATA 10.16 ,12.9 ,12.02, .115, .167, 312
DATA  .00432, 25.0, 90, 80, 1, 0
DATA 9680, 41.8, 25.58, .472, .640, 36, 4.06
DATA .0174533, .0406, 6.40, 1400, .2, .0635, .1016
DATA .0510, 2.500, 3.500, 6, .2, 13.65, 2.325
DATA 4, 90, 1200, 135, 2000, 435, .1532
DATA 0.76, 71, 6, .4, .006894, 2.54, 60
DATA 8.314, 1.341e-3, .0406, 0
'     , 179
DATA bottom

'==================================================================