Tew at al Stirling engine third order calculation Quick Basic program
The engine of the future project
The following is a translation into Quick Basic of the Computer Program for Stirling Engine Performance Calculation, written in Fortran.
I found it in the A_Stirling_Engine_Computer_Model_for_Per.pdf
DOE/NASA/1011-78/24
NASA TM-78884
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19780022051.pdf
I have translated it into Quick Basic.
In the future I will debug this program and write a version that makes more use of the input and output features that modern programming languages offer. 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.
This page is meant for progammers.
For the FORTRAN source code click here
Click here for HOME
Problems
There are still bugs in this program.
F(4) and F(9) remain zero.
The iteration counter ITER& goes up and up,
and related to this the outputs are not written as it should.
Outputs are written in the file TEWATAL.TXT
If you find the bugs, please e-mail me, hvd@haghen.nl
'====================================================================
DECLARE SUB closefilenr (filecount%, filenr%)
DECLARE SUB givefilenr (filecount%, filenr%)
DECLARE SUB HEATX (TIME!, TG!(), TGA!(), TMA!(), REYNO!(), QR!(), X!(), RODY!, JCYCLE%, diamd!, X1MIN!, RODYMX!, V!(), SIG3!, W!, MWGAS%, RHCFAC!, REALG!, CONDID!, NOCYC%)
DECLARE SUB CYCL (TIME!, UTOTAL!, STROKV!, AVGWSP!)
DECLARE SUB CNDCT (TMA!(), QCREG!, QCYL!, QCONDD!, QCAN!, QSHTTL!, MWGAS%)
DECLARE SUB XDELP (K%, AREAIN!, AREAOT!, DHX!, XLGTH!, ROX!, VISX!, FLOWIN!, PXIN!, XCOEF!, DP!, M%)
DEFINT I-N
'input variables
DATA "COEFF"
DATA "leage coefficient for flow between buffer and compression spaces"
DATA "FACT1"
DATA "FACT1 and FACT2 are convergence constants used in regenerator
metal temperature convergence scheme"
DATA "FACT2"
DATA "FACT2 and FACT1 are convergence constants used in regenerator
metal temperature convergence scheme"
DATA "FIPCYC"
DATA "number of iterations per cycle"
DATA "IPRINT"
DATA "number of iterations between printouts (if JIP%=0)"
DATA "ITMPS"
DATA "1, printout gas and metal control volume temperature at every iteration
0, do not do the above"
DATA "JIP%"
DATA ">0, short-form printout
=0, long-form printout"
DATA "MWGAS%"
DATA "4, use helium working gas
not 4, use hydrogen working gas"
DATA "NOCYC%"
DATA "number of engine cycles to be run"
DATA "NOEND"
DATA"number of cycle at wich convergence scheme is turned on"
DATA "NSTRT"
DATA "number of cycle at which regenerator temperature convergence
scheme is turned on"
DATA "OMEGA"
DATA "drive speed, Hz"
DATA "P"
DATA "initial pressure, lbf/inch^2 (N/m^2) (approximately equal to mean
pressure for GPU)"
DATA "P3"
DATA "initial buffer space pressure, lbf/inch^2 (N/m^2)"
DATA "REALGS"
DATA "1, use real-gas equation of state
0, use ideal-gas equation of state"
DATA "RHCFAC"
DATA "regenerator heat-transfer coefficient multiplication factor"
DATA "TG"
DATA "TG(12) control volume interface gas temperatures, degrees R (Kelvin)"
DATA "TGA"
DATA "TGA(13) control volume gas temperatures, degrees R (Kelvin)"
DATA "TM"
DATA "TM(13) metal boundary temperatures, degrees R (Kelvin)"
DATA "T3"
DATA "buffer space temperature (assumed constant), degrees R (Kelvin)"
'output variables
DATA "AVGWSP"
DATA "mean working-space pressure, lbf/inch^2 (N/m^2)"
DATA "EFFTOT"
DATA "efficiency, WRKTOT/QIN"
DATA "ENTH"
DATA "net enthalpy flow from working space to buffer space
over one cycle, ft-lbf (J)"
DATA "EN3PM"
DATA "rate of displacement of working-space gas, inch^3/min (cm^3/min)"
DATA "FREQ"
DATA "rotational frequency of crankshaft, Hz"
DATA "PWRHP"
DATA "indicated power -WRKTOT/cycle period, hp (kW)"
DATA "QCNDCL"
DATA "heat conduction through cylinder, ft-lbf/cycle (J/cycle)"
DATA "QCNDCN"
DATA "heat conduction through insulation container, ft-lbf/cycle (J/cycle)"
DATA "QCNDD"
DATA "heat conduction through displacer, ft-lbf/cycle (J/cycle)"
DATA "QCNDRI"
DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)"
DATA "QCNDRO"
DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)"
DATA "QCOMPN"
DATA "heat from compression-space wall to gas, ft-lbf/cycle (J/cycle)"
DATA "QCOMPP"
DATA "heat from gas to compression-space wall, ft-lbf/cycle (J/cycle"
DATA "QCOOLN"
DATA "heat from cooler to gas, ft-lbf/cycle (J/cycle)"
DATA "QCOOLP"
DATA "heat from gas to cooler, ft-lbf/cycle (J/cycle)"
DATA "QEXPN"
DATA "heat from metal wall to gas in expansion space, ft-lbf/cycle (J/cycle)"
DATA "QEXPP"
DATA "heat from gas to metal wall in expansion space ft-lbf/cycle (J/cycle)"
DATA "QHEATN"
DATA "heat from heater to gas, ft-lbf/cycle (J/cycle)"
DATA "QHEATP"
DATA "heat from gas to heater, ft-lbf/cycle (J/cycle)"
DATA "QIN"
DATA "heat into gas in expansion space and heater tubes over one cycle,
ft-lbf (J)"
DATA "QOUT"
DATA "heat out of gas in cooler tubes and compression space over
one cycle, ft-lbf (J)"
DATA "QREGEN"
DATA "net heat flow from gas to regenerator metal over one cycle,
ft-lbf/cycle (J/cycle)"
DATA "QREGN"
DATA "heat from regenerator metal to gas, ft-lbf/cycle (J/cycle)"
DATA "QREGP"
DATA "heat from gas to regenerator metal, ft-lbf/cycle (J/cycle)"
DATA "QSHTL"
DATA "shuttle heat loss, ft-lbf/cycle (J/cycle)"
DATA "REFF1"
DATA "like REFF2 a measure of regenerator effectiveness"
DATA "REFF2"
DATA "like REFF1 a measure of regenerator effectiveness"
DATA "TG"
DATA "gas temperature at interfase of two adjacent control volumes
degrees R (kelvin)"
DATA "TGA"
DATA "working bulk gas temperature of a control volume, degrees R (kelvin)"
DATA "TGACYC"
DATA "cycle-time average temperature of TGA, degrees R (kelvin)"
DATA "TGCMPA"
DATA "cycle-time average temperature in compression space, degrees R (kelvin)"
DATA "TGCYC"
DATA "cycle time average temperature of TG, degree R (kelvin)"
DATA "TGEXPA"
DATA "cycle time average temperature in expansion space, degree R (kelvin)"
DATA "TM"
DATA "metal temperature of a metal control volume, degree R (kelvin)"
DATA "TMCYC"
DATA "cycke-time average temperarure of TM, degree R (kelvin)"
DATA "UTOTAL"
DATA "internal energy content of working-space gas, ft-lbf (J)"
DATA "WDISP"
DATA "net work done on displacer piston by working-space gas, ft-lbf (J)"
DATA "WDISPN"
DATA "work done on gas by displacer, ft-lbf/cycle (J/cycle)"
DATA "WDISPP"
DATA "work done on piston by gas, ft-lbf/cycle (J/cycle)"
DATA "WPISTP"
DATA "work done on power piston by gas, ft-lbf/cycle (J/cycle"
DATA "WPIST"
DATA "work done on power piston by working-space gas ft-lbf (J)"
DATA "WPISTN"
DATA "work done on gas by piston, ft-lbf/cycle ft-lbf (J/cycle)"
DATA "WRKCMP"
DATA "work done by the gas in compression space over one cycle, ft-lbf (J)"
DATA "WRKEXP"
DATA "work done by the gas in compression space over one cycle, ft-lbf (J)"
DATA "WRKLP"
DATA "power-piston work lost due to pressure drop from regenerator
to compression space, ft-lbf/cycle J/cycle"
DATA "WRKLT"
DATA "work lost due to pressure drop over one cycle, ft-lbf (J)"
DATA "WORKTOT"
DATA "net work done by gas in working space over one cycle,
WRKEXP + WRKCMP, ft-lbf (J)"
'PRINT X(2), V(13), reyno(13), qr(5), tmax(12)
'PRINT tg(12), tga(13), tm(13), tmavg(13), DTM(13)
'PRINT dtime, f(1), wg(1), delp(1), p, pe, pC
'PRINT cyc, qexp, qheatr, qregen, qcoolr, qcomp, enfrth, enfhtr
'PRINT enfctr, enfrtc, wrkp, wrkexp
'PRINT wrkcmp, wrklex, wrklcm, wrklp, wrkh, sumnum(1), sumden(1), enth
'PRINT wdisp, wpist, wtpist, qrxpn, qexpp, qheatn, qheatp
'PRINT qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
'PRINT wpistn, qcndri, qcndro, qcndcl, qcndd, qcndcn, qshtl, tgexpa, tgcmpa
'PRINT deltap, ad, ap, ha(1)
'PRINT ctit, ncurv, iset1, iset2, iset3
'PRINT jndex, JIP%, NOCYC%
'PRINT strlng, p, realgs, fipcyc, iprint, itmps, coeff, t3, p3, omega
'PRINT fact1, fact2, NOCYC%, nstrt, noend, MWGAS%, tg, tga, tm, rhfac, JIP%
' PAGE 47 (53 IN MANUAL PDF)----------------
' MAIN PROGRAM
' DOES INITIALIZATION. UPDATES PISTON POSITIONS--EXPANSION, COMPRESSION
' AND BUFFER SPACE VOLUMES. INTEGRATES P*DELTA(V) TO DETERMINE WORK.
' DETERMINES WHEN CYCLES IS COMPLETE.
DIM X(2), V(13), REYNO(13), QR(5), tmax(12)
DIM TG(12), TGA(13), TM(13), tmavg(13), DTM(13)
' COMMON DTIME, F(), WG(), DELP(), P, PE, PC
' COMMON /cyc/ qexp, qheatr, qregen, qcoolr, qcomp, enfrth, enfhtr
' COMMON enfctr, enfrtc, wrkp, wrkexp
' COMMON wrkcmp, wrklex, wrklcm, wrklp, wrkh, sumnum(), sumden(), enth
' COMMON wdisp, wpist, wtpist, QEXPN, qexpp, qheatn, qheatp
' COMMON qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
' COMMON wpistn, qcndri, qcndro, qcndcl, qcndd, qcndcn, qshtl, tgexpa, tgcmpa
' COMMON /deltap/ ad, ap, ha()
' COMMON /TAVCYC/ TGCYC(), TGACYC(), TMCYC()
' COMMON /ctit/ ncurv, iset1, iset2, iset3
' COMMON /jndex/ JIP%, NOCYC%
DIM F(12), WG(13), DELP(13)
COMMON DTIME, F(), WG(), DELP(), P, PE, PC
COMMON /CYC/ QEXP, qheatr, qregen, qcoolr, qcomp, ENFRTH, ENFHTR
COMMON enfctr, enfrtc, wrkp, wrkexp
DIM SUMNUM(13), SUMDEN(13)
COMMON wrkcmp, WRKLEX, wrklcm, wrklp, wrkh, SUMNUM(), SUMDEN(), enth
COMMON wdisp, wpist, wtpist, QEXPN, qexpp, qheatn, qheatp
COMMON qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
COMMON wpistn, qcndri, qcndro, qcndcl, qcndd, qcndcn, QSHTL, tgexpa, tgcmpa
DIM HA(13)
COMMON /deltap/ AD, AP, HA()
DIM TGCYC(12), TGACYC(13), TMCYC(13)
COMMON /TAVCYC/ TGCYC(), TGACYC(), TMCYC()
COMMON /ctit/ ncurv, iset1, iset2, iset3
' COMMON /jndex/ JIP%, NOCYC%
COMMON JIP%, NOCYC%
DIM VAR(106), STORE(11, 86)
COMMON VAR(), STORE()
DIM vin$(21)
DIM vit$(21)
DIM vi(55) ' 15 + 12 + 13 + 13 + 2)
DIM von$(48)
DIM vot$(48)
COMMON von$(), vot$()
COMMON outfilenr%
' NAMELIST /STRLNG/P,REALGS,FIPCYC,IPRINT,ITMPS,COEFF,T3,P3,OMEGA,
'FACT1,FACT2,NOCYC%,NSTRT,NOEND,MWGAS%,TG,TGA,TM,RHCFAC,JIP%
DATA AD , AP , AR , APISTR
DATA 6.0, 6.0, 0.11, 0.60
DATA RCRANK, E , RODL, V3CL , DIAMD, DIAMP, DIAMDR, DIAMPR
DATA 0.543 , 0.82, 1.81, 24.84, 2.75 , 2.75 , 0.375, 0.875
DATA 0.0, 2*1.5948, .933, 5*0.624, 3*0.240, 0.0
'cls
TRUE = 1
FALSE = -1
CLS
GOSUB 1900 'initialize input vars
GOSUB 1950 'initialize output vars
'48
RESTORE 4000
CALL givefilenr(filecount%, outfilenr%)
outfilename$ = "TEWATAL.TXT"
OPEN outfilename$ FOR OUTPUT AS outfilenr%
' GOSUB 1910 'print INPUTS TO output file, input vars from aray
10 ' CONTINUE loop 1 <<<<<<<<<<<<<<<<<<<<<<<<<<<
loop1 = loop1 + 1
loop2 = 0
loop3 = 0
LOCATE 1, 1
PRINT "loop1="; loop1; " loop2=";
IF loop2 < 100 THEN PRINT " ";
IF loop2 < 10 THEN PRINT " ";
PRINT loop2; " loop3=";
IF loop3 < 100 THEN PRINT " ";
IF loop3 < 10 THEN PRINT " ";
PRINT loop3; " NOCYC%="; NOCYC%; " JCOUNT%="; JCOUNT%;
PRINT " ITER&="; ITER&; " "
' READ IN INITIAL PRESSURES, TEMPERATURES AND OTHER PARAMETERS WHICH
' DEFINE THE NATURE OF THE RUN TO BE MADE
'READ (5,STRLNG,END=1000)
READ I$
IF I$ = "END" OR I$ = "end" THEN
GOTO 1000' end OF PROGRAM
END IF
P = VAL(I$)
READ REALGS, FIPCYC, IPRINT, ITMPS, COEFF, T3, P3, OMEGA
READ FACT1, FACT2, NOCYC%, NSTRT, NOEND, MWGAS%
FOR I = 1 TO 12
READ TG(I)
NEXT I
FOR I = 1 TO 13
READ TGA(I)
NEXT I
FOR I = 1 TO 13
READ TM(I)
NEXT I
READ RHCFAC, JIP%
'IF PRINTINP = FALSE THEN
GOSUB 1920 'print output file, inputs from variables
' PRINTINP = TRUE
'END IF
AD = 6!
AP = 6!
AR = .11
APISTR = .6
rcrank = .543' crank radius in inch
E = .82 ' eccentricity Rhombic drive in inch
rodl = 1.81 ' connecting rod length in inch
v3cl = 24.84
diamd = 2.75 ' displacer diameter (inch)
diamp = 2.75 ' piston diameter (inch)
diamr = .375 ' displacer-rod diameter (inch)
diampr = .875 ' piston-rod diameter (inch)
'volumes
V(1) = 0! ' expansion space
V(2) = 1.5948' heater tubes part 1
V(3) = 1.5948' heater tubes part 2
V(4) = .933 ' heater tubes part 3
V(5) = .624' hot connecting duct?
V(6) = .624' regenerator
V(7) = .624' regenerator
V(8) = .624' regenerator
V(9) = .624' cold connecting duct ?
V(10) = .24' cooler tubes part 1
V(11) = .24' cooler tubes part 2
V(12) = .24' cooler tubes part 3
V(13) = 0!' compression space
DTIME = 1! / (OMEGA * FIPCYC)
' PRINT ' ' (6,STRLNG)
' INITIALIZATION
FMULT = 1!
ITER& = 0
TIME = 0!
PSI = 0!
PSIDEG = 0!
SAVET = 0!
IPLOT = 0
PRSUM = 0!
NOITPC = 0
' DO 11 I=1,13
FOR I = 1 TO 13
IF (I) = (13) THEN GOTO 12
' "I="; I
' PRINT "TGCYC(I)="; TGCYC(I)
' PRINT "TG(I)="; TG(I)
TGCYC(I) = TG(I)
12 TGACYC(I) = TGA(I)
11 TMCYC(I) = TM(I)
NEXT I
CPCYC = 1! / (OMEGA * DTIME)
' R,CP,CV IN UNITS OF IN-LBF/ LBM-DEG. R
GAMMA = 1.394
R = 9197!
CP = 32557!
CV = 23360!
B = .02358
IF (MWGAS%) <> (4) THEN GOTO 15
B = .01613
'PAGE 48 (54 IN MANUAL PDF)-----------------------
GAMMA = 1.667
R = 4632!
CP = 11580!
CV = 6948!
15 ' CONTINUE
IP = IPRINT
IPRNT2 = IPRINT / 25
ADR = AD - AR
APR = AP - APISTR
X1MIN = SQR((rodl - rcrank) ^ 2 - E ^ 2)
RODYMX = SQR(rodl ^ 2 - (E - rcrank) ^ 2)
RODYMN = SQR(rodl ^ 2 - (E + rcrank) ^ 2)
STROKV = 2! * (RODYMX - RODYMN)
RODY = RODYMX
X(1) = RODYMX
X(2) = -RODYMX
V(1) = AD * (X(1) - X1MIN) + .965
' .965=expansion space dead volume (inch^3)?
' expansion space clearence=.163 cm= .064 inch
V(13) = 2! * ADR * (RODYMX - RODY) + .353
' .353=compression space dead volume (inch^3) ?
' compression space clearance=.030 cm= 0.012 inch
V3 = APR * (-X(2) - X1MIN) + v3cl
V3STAR = V3
VOVRT = 0!
' DO 30 I=1,13
FOR I = 1 TO 13
30 VOVRT = VOVRT + V(I) / (TGA(I) + B * P * REALGS)
NEXT I
W3 = P3 * V3STAR / (R * (T3 + B * P3 * REALGS))
ALPHA = W3 * R / VOVRT
SIG3 = ALPHA / (P + ALPHA)
W = W3 / SIG3
PE = P
PC = P
' DO 50 I=1,13
FOR I = 1 TO 13
50 WG(I) = P * V(I) / (R * (TGA(I) + B * P * REALGS))
IF WG(I) = 0 THEN
PRINT " V(I)="; V(I); " TGA(I)="; TGA(I); " B="; B; " P="; P; " REALGS="; REALGS
INPUT I$
END IF
NEXT I
ISTART = 1
JCOUNT% = 1
100 ' CONTINUE 'no goto 100 in main
IF (JCOUNT%) < (ISTART) THEN GOTO 140 'LOOP3 >>>>>>>>>>>>>>>>>>>>
125 'loop2 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
loop2 = loop2 + 1
loop3 = 0
LOCATE 1, 1
PRINT "loop1="; loop1; " loop2=";
IF loop2 < 100 THEN PRINT " ";
IF loop2 < 10 THEN PRINT " ";
PRINT loop2; "loop3=";
IF loop3 < 100 THEN PRINT " ";
IF loop3 < 10 THEN PRINT " ";
PRINT loop3; " NOCYC%="; NOCYC%; " JCOUNT%="; JCOUNT%;
PRINT " ITER&="; ITER&; " "
IF (JIP%) > (0) THEN GOTO 140 'LOOP3 >>>>>>>>>>>>>>>
IF (JCOUNT%) < ((NOCYC% - 1)) THEN ' (6,130)
130 'FORMAT ('2 TIME ANGLE X(1) X(2) T4',7X,'T9',7X,'T1',4X,'T2
'1 ', 5X,' PE PC PRIN PROUT ',3X,'F(4)',5X,'F(9) ')
PRINT "2 TIME ANGLE X(1)="; X(1); " X(2)="; X(2); " T4="; T4;
PRINT " PE="; PE; " PC="; PC; " PRIN="; PRIN; " PROUT="; PROUT;
PRINT #outfilenr%, "2 TIME ANGLE X(1)="; X(1); " X(2)="; X(2); " T4="; T4;
PRINT #outfilenr%, " PE="; PE; " PC="; PC; " PRIN="; PRIN; " PROUT="; PROUT;
END IF
IF (JCOUNT%) = ((NOCYC% - 1)) THEN ' (6,132)
132 'FORMAT (' TIME ANGLE RE1 RE2 RE3 RE4 RE5'
'1' RE6 RE7 RE8 RE9 RE10 RE11 RE12
'2 'RE13 ')
'PRINT " TIME ANGLE RE1 RE2 RE3 RE4 RE5";
'PRINT " RE6 RE7 RE8 RE9 RE10 RE11 RE12 ";
'PRINT "RE13 "
PRINT "TIME="; TIME; " ANGLE="; ANGLE; " RE1="; RE1; " RE2="; RE2;
PRINT " RE3="; RE3; " RE4="; RE4; " RE5="; RE5; " RE6="; RE6;
PRINT " RE7="; RE7; " RE8="; RE8; " RE9="; RE9; " RE10="; RE10;
PRINT " RE11="; RE11; " RE12="; RE12; " RE13= "; RE13
PRINT #outfilenr%, "TIME="; TIME; " ANGLE="; ANGLE; " RE1="; RE1; " RE2="; RE2;
PRINT #outfilenr%, " RE3="; RE3; " RE4="; RE4; " RE5="; RE5; " RE6="; RE6;
PRINT #outfilenr%, " RE7="; RE7; " RE8="; RE8; " RE9="; RE9; " RE10="; RE10;
PRINT #outfilenr%, " RE11="; RE11; " RE12="; RE12; " RE13= "; RE13
END IF
IF (JCOUNT%) = (NOCYC%) THEN ' (6,133)
133 'FORMAT (' TIME ANGLE F1 F2 F3 F4 F5'
'1' F6 F7 F8 F9 F10 F11 F12')
'PRINT " TIME ANGLE F1 F2 F3 F4 F5";
'PRINT " F6 F7 F8 F9 F10 F11 F12"
PRINT " TIME="; TIME; " ANGLE="; ANGLE;
PRINT " F1="; F1; " F2="; F2; " F3="; F3; " F4="; F4; " F5="; F5;
PRINT " F6="; F6; " F7="; F7; " F8="; F8; " F9="; F9;
PRINT " F10="; F10; " F11="; F11; " F12="; F12
PRINT #outfilenr%, " TIME="; TIME; " ANGLE="; ANGLE;
PRINT #outfilenr%, " F1="; F1; " F2="; F2; " F3="; F3; " F4="; F4; " F5="; F5;
PRINT #outfilenr%, " F6="; F6; " F7="; F7; " F8="; F8; " F9="; F9;
PRINT #outfilenr%, " F10="; F10; " F11="; F11; " F12="; F12
END IF
' BEGINNING OF LOOP
140 'loop 3 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
loop3 = loop3 + 1
LOCATE 1, 1
PRINT "loop1="; loop1; " loop2=";
IF loop2 < 100 THEN PRINT " ";
IF loop2 < 10 THEN PRINT " ";
PRINT loop2; " loop3=";
IF loop3 < 100 THEN PRINT " ";
IF loop3 < 10 THEN PRINT " ";
PRINT loop3; " NOCYC%="; NOCYC%; " JCOUNT%="; JCOUNT%;
PRINT " ITER&="; ITER&; " "
ITER& = ITER& + 1
IF (JCOUNT%) >= ((NOCYC% - 4)) THEN IPRINT = IPRNT2
IF ((TIME + SAVET)) > (5!) THEN GOTO 700
IP = IP + 1
VTOTL = 0!
' DO 142 IJK=1,13
FOR IJK = 1 TO 13
142 VTOTL = VTOTL + V(IJK)
NEXT IJK
IF (JIP%) > (0) THEN GOTO 170
IF (IP) < (IPRINT) THEN GOTO 180
IF (JCOUNT%) < (ISTART) THEN GOTO 170
IF (JCOUNT%) < ((NOCYC% - 1)) THEN ' PRINT ' ' (6,160) TIME,PSIDEG,X,TG(4),TG(9),
PRINT " TIME="; TIME; " PSIDEG="; PSIDEG; " X="; X
PRINT " TG(4)="; TG(4); " TG(9)="; TG(9); " TGA(1)="; TGA(1)
PRINT " TGA(13)="; TGA(13); " PE="; PE; " PC="; PC
PRINT " PRIN="; PRIN; " PROUT="; PROUT; " F(4)="; F(4); " F(9)="; F(9)
PRINT #outfilenr%, " TIME="; TIME; " PSIDEG="; PSIDEG; " X="; X
PRINT #outfilenr%, " TG(4)="; TG(4); " TG(9)="; TG(9); " TGA(1)="; TGA(1)
PRINT #outfilenr%, " TGA(13)="; TGA(13); " PE="; PE; " PC="; PC
PRINT #outfilenr%, " PRIN="; PRIN; " PROUT="; PROUT; " F(4)="; F(4); " F(9)="; F(9)
END IF
160 'FORMAT (1X,F7.4,F8.1,2F6.2,1X,2F8.2,F8.0,F6.0,1X,4F8.2,2F8.5)
IF (JCOUNT%) = ((NOCYC% - 1)) THEN
PRINT " TIME="; TIME; ",PSIDEG="; PSIDEG; ",REYNO="; REYNO
PRINT #outfilenr%, " TIME="; TIME; ",PSIDEG="; PSIDEG; ",REYNO="; REYNO
END IF
162 'FORMAT ' (1X,F7.4,F8.1,13F8.0)
IF (JCOUNT%) = (NOCYC%) THEN
PRINT " TIME="; TIME; ",PSIDEG="; PSIDEG; ",F="; F
PRINT #outfilenr%, " TIME="; TIME; ",PSIDEG="; PSIDEG; ",F="; F
END IF
163 'FORMAT (1X,F7.4,F8.1,12F8.5)
IF (ITMPS) = (1) THEN
PRINT " TG", "TGA", "TM"
FOR NBR% = 1 TO 12
PRINT NBR%; " "; TG(NBR%), TGA(NBR%), TM(NBR%)
NEXT NBR%
PRINT , TGA(13), TM(13)
PRINT #outfilenr%, " TG", "TGA", "TM"
FOR NBR% = 1 TO 12
PRINT #outfilenr%, NBR%; " "; TG(NBR%), TGA(NBR%), TM(NBR%)
NEXT NBR%
PRINT #outfilenr%, , TGA(13), TM(13)
END IF
660 'FORMAT (' TG=',4X,12F8.0/,' TGA=',13F8.0/,' TM=',13F8.0/)
' PAGE 49 (55 IN MANUAL PDF)--------------------------------
170 IP = 0
'180 TIME = DTIME * FLOAT(ITER&)
180 TIME = DTIME * (ITER&)
' IF (ITER&) = (IFIX(.3 / DTIME)) AND (JCOUNT%) < (6) THEN GOTO 700
' IF (ITER& = FIX(.3 / DTIME)) AND (JCOUNT% < 6) THEN GOTO 700
IF (ITER& >= FIX(.3 / DTIME)) AND (JCOUNT% < 6) THEN GOTO 700
' CALL HEAT EXCHANGER SUBROUTINE TO UPDATE GAS LUMP MASSES, FLOW RATES,
' TEMPERATURES, PRESSURE + REGENERATOR METAL TEMPERATURES.
J% = JCOUNT% 'JCOUNT% is used in stead of J%
'J% is used a lot as loop counter in HEATX (Harrie )
' CONDID IS NOT SET AND NOT USED ??????????????????????????????
CALL HEATX(TIME, TG(), TGA(), TM(), REYNO(), QR(), X(), RODY, J%, diamd, X1MIN, RODYMX, V(), SIG3, W, MWGAS%, RHCFAC, REALGS, CONDID, NOCYC%)
' JCOUNT% = J% '???? REPAIR DEBUG
PRSUM = PRSUM + P
NOITPC = NOITPC + 1
' UPDATE EXPANSION AND COMPRESSION SPACE PRESSURES. (LBF/IN2)
PEOLD = PE
PCOLD = PC
POLD = P
IF (JCOUNT%) < ((NOCYC% - 5)) THEN GOTO 240
PRIN = (DELP(5) + DELP(6) + .5 * DELP(7)) * FMULT + P
PROUT = P - (.5 * DELP(7) + DELP(8) + DELP(9)) * FMULT
PE = (DELP(2) + DELP(3) + DELP(4)) * FMULT + PRIN
PC = PROUT - (DELP(10) + DELP(11) + DELP(12)) * FMULT
GOTO 245
240 PE = P
PC = P
PRIN = P
PROUT = P
245 ' CONTINUE
AE = (PE + PEOLD) / 2!
AC = (PC + PCOLD) / 2!
a = (P + POLD) / 2!
DELPRG = PRIN - PROUT
DELPHT = PE - PRIN
DELPCL = PROUT - PC
' CALCULATE INTERNAL ENERGY OF WORKING SPACE GAS. (FT-LBF)
UTOTAL = 0!
' DO 250 I=1,13
FOR I = 1 TO 13
250 UTOTAL = WG(I) * CV * TGA(I) / 12! + UTOTAL
NEXT I
260 ' CONTINUE
' BUFFER SPACE BLEED FLOW,FLOTOB (LBM/SEC)
IF ((PC - P3)) >= (0!) THEN GOTO 270
SIG3P = -COEFF * (P3 - PC) ^ .5 / W
GOTO 280
270 SIG3P = COEFF * (PC - P3) ^ .5 / W
280 ' CONTINUE
FLOTOB = SIG3P * W
IF (FLOTOB) >= (0) THEN ENTHWB = CP * TGA(13) * FLOTOB * DTIME / 12!
IF (FLOTOB) < (0) THEN ENTHWB = CP * T3 * FLOTOB * DTIME / 12!
enth = enth + ENTHWB
SIG3 = SIG3 + SIG3P * DTIME
WG(13) = WG(13) - FLOTOB * DTIME
IF (FLOTOB) < (0!) THEN
TGA(13) = TGA(13) - FLOTOB * DTIME * (T3 - TGA(13)) / WG(13)
END IF
520 P3OLD = P3
P3 = T3 / (V3 / (SIG3 * W * R) - B * REALGS)
A3 = (P3OLD + P3) / 2!
' UPDATE CRANK ANGLE, PSI(RADIANS), PSIDEG(DEGREES)
PSIOLD = PSI
PSI = PSI + OMEGA * DTIME * 6.28318
IF PSI >= 6.28318 THEN PSI = PSI - 6.28318
PSIDEG = PSI * 360! / 6.28318
' UPDATE PISTON POSITIONS (IN)
RODY = SQR(rodl ^ 2 - (E - rcrank * COS(PSI)) ^ 2)
X1OLD = X(1)
X2OLD = X(2)
X(1) = rcrank * SIN(PSI) + RODY
X(2) = rcrank * SIN(PSI) - RODY
' UPDATE EXPANSION, COMPRESSION, AND BUFFER SPACE VOLUMES. (IN3)
V1OLD = V(1)
V13OLD = V(13)
' PAGE 50 (56 IN MANUAL PDF)-----------------------------------
V(1) = AD * (X(1) - X1MIN) + .965
V(13) = 2! * ADR * (RODYMX - RODY) + .353
V3 = APR * (-X(2) - X1MIN) + v3cl
' CALCULATE WORK (FT-LBF)--EXPANSION SPACE WORK,WRKEXP; COMPRESSION
' SPACE WORK,WRKCMP; TOTAL WORK; DISPLACER WORK,WDISP; POWER PISTON
' WORK,WPIST; PRESSURE DROP WORK LOSS IN EXPANSION AND COMPRESSION
' SPACE, WRKLEX AND WRKLCM
wrkp = wrkp + AC * (V(13) - V13OLD) / 12!
DWEXP = AE * (V(1) - V1OLD) / 12!
wrkexp = wrkexp + DWEXP
IF (DWEXP) > (0!) THEN WEXPP = WEXPP + DWEXP
IF (DWEXP) < (0!) THEN WEXPN = WEXPN + DWEXP
DWCMP = AC * (V(13) - V13OLD) / 12!
wrkcmp = wrkcmp + DWCMP
IF (DWCMP) > (0!) THEN WCMPP = WCMPP + DWCMP
IF (DWCMP) < (0!) THEN WCMPN = WCMPN + DWCMP
WRKLEX = WRKLEX + (a - AE) * (V(1) - V1OLD) / 12!
wrklcm = wrklcm + (a - AC) * (V(13) - V13OLD) / 12!
wrklp = wrklp + (a - AC) * (V(13) - V13OLD) / 12!
wrkh = wrkh + (AE * (V(1) - V1OLD) + AC * (V(13) - V13OLD)) / 12!
DWDISP = (AE * AD - AC * ADR) * (X(1) - X1OLD) / 12!
wdisp = wdisp + DWDISP
IF (DWDISP) > (0!) THEN wdispp = wdispp + DWDISP
IF (DWDISP) < (0!) THEN wdispn = wdispn + DWDISP
DWPIST = AC * ADR * (X(2) - X2OLD) / 12!
wpist = wpist + DWPIST
IF (DWPIST) > (0!) THEN wpistp = wpistp + DWPIST
IF (DWPIST) < (0!) THEN wpistn = wpistn + DWPIST
wtpist = wtpist + AC * ADR * (X(2) - X2OLD) / 12!
' DO 530 I=1,13
FOR I = 1 TO 13
IF (I) = (13) THEN GOTO 531
TGCYC(I) = TGCYC(I) + TG(I)
531 TGACYC(I) = TGACYC(I) + TGA(I)
530 TMCYC(I) = TMCYC(I) + TM(I)
NEXT I
IF (PSIOLD > 3.14159) AND (PSI < 3.14159) THEN
GOTO 650' jump out of the loop3 when cycle is ready
END IF
GOTO 140 ' LOOP3 >>>>>>>>>>>>>>>>>>>>>>
' END OF LOOP
650 ' CONTINUE
'PRINT "AVERAGE TEMPERATURES"
'INPUT I$
' CALCULATE AVERAGE GAS TEMPERATURES OVER THE CYCLE FOR EACH CONTROL VOLUME
' DO 533 I=1,13
FOR I = 1 TO 13
IF (I) = (13) THEN GOTO 534
TGCYC(I) = TGCYC(I) / CPCYC
534 TGACYC(I) = TGACYC(I) / CPCYC
533 TMCYC(I) = TMCYC(I) / CPCYC
NEXT I
tgexpa = TGACYC(1)
tgcmpa = TGACYC(13)
' DO 667 II=1,6
FOR II = 1 TO 6
JJ = II + 6
tmax(II) = 0!
667 tmax(JJ) = 5000!
NEXT II
IF (JIP% > 0) AND (JCOUNT% <> (NOCYC% - 1)) THEN GOTO 668
IF (JCOUNT% < (NOCYC% - 1)) THEN
'WRITE (6,160) TIME,PSIDEG,X,P,P3,TGA(1),
'1TGA(13),PE,PC,PRIN,PROUT,V(1),V(13)
PRINT " TIME="; TIME; ",PSIDEG="; PSIDEG
PRINT ",X="; X; ",P="; P; ",P3="; P3
PRINT ",TGA(1)="; TGA(1); ", TGA(13)="; TGA(13)
PRINT ", PE="; PE; ", PC="; PC
PRINT ", PRIN="; PRIN; ", PROUT="; PROUT
PRINT ", V(1)="; V(1); ", V(13)="; V(13)
PRINT #outfilenr%, " TIME="; TIME; ",PSIDEG="; PSIDEG
PRINT #outfilenr%, ",X="; X; ",P="; P; ",P3="; P3
PRINT #outfilenr%, ",TGA(1)="; TGA(1); ", TGA(13)="; TGA(13)
PRINT #outfilenr%, ", PE="; PE; ", PC="; PC
PRINT #outfilenr%, ", PRIN="; PRIN; ", PROUT="; PROUT
PRINT #outfilenr%, ", V(1)="; V(1); ", V(13)="; V(13)
END IF
IF (JIP% = 0) THEN
PRINT " TG, TGA, TM"
FOR NBR% = 1 TO 12
PRINT NBR%; " "; TG(NBR%), TGA(NBR%), TM(NBR%)
NEXT NBR%
PRINT , TGA(13), TM(13)
PRINT #outfilenr%, " TG, TGA, TM"
FOR NBR% = 1 TO 12
PRINT #outfilenr%, NBR%; " "; TG(NBR%), TGA(NBR%), TM(NBR%)
NEXT NBR%
PRINT #outfilenr%, , TGA(13), TM(13)
END IF
668 ' CONTINUE
IF (JCOUNT% < NSTRT) OR (JCOUNT% > NOEND) THEN GOTO 675
' CONVERGENCE SCHEME FOR REGENERATOR METAL TEMPERATURES.
' DO 670 JI=5.9
FOR JI = 5 TO 9
670 DTM(JI) = FACT1 * DTM(JI) + FACT2 * SUMNUM(JI) / SUMDEN(JI)
NEXT JI
TM(5) = TM(5) - (5! * DTM(5) + 4! * DTM(6) + 3! * DTM(7) + 2! * DTM(8) + DTM(9)) / 3!
TM(6) = TM(6) - (4! * DTM(5) + 8! * DTM(6) + 6! * DTM(7) + 4! * DTM(8) + 2! * DTM(9)) / 3!
TM(7) = TM(7) - (1! * DTM(5) + 2! * DTM(6) + 3! * DTM(7) + 2! * DTM(8) + DTM(9))
TM(8) = TM(8) - (2! * DTM(5) + 4! * DTM(6) + 6! * DTM(7) + 8! * DTM(8) + 4! * DTM(9)) / 3!
TM(9) = TM(9) - (DTM(5) + 2! * DTM(6) + 3! * DTM(7) + 4! * DTM(8) + 5! * DTM(9)) / 3!
GOTO 675
675 ' CONTINUE
' PAGE 51 (57 IN MANUAL PDF)-------------
' DO 673 I=1,5
FOR I = 1 TO 5
673 QR(I) = 0!
NEXT I
AVGWSP = PRSUM / NOITPC
CALL CYCL(TIME, UTOTAL, STROKV, AVGWSP)
' DO 532 I=1,13
FOR I = 1 TO 13
IF (I) = (13) THEN GOTO 535
TGCYC(I) = 0!
535 TGACYC(I) = 0!
532 TMCYC(I) = 0!
NEXT I
PRSUM = 0!
NOITPC = 0
JCOUNT% = JCOUNT% + 1
' PRINT "########## JCOUNT%="; JCOUNT%
' INPUT I$
IF (JCOUNT% <> (NOCYC% - 1)) THEN GOTO 680
PRINT "### NOCYC% - 1 ####### JCOUNT%="; JCOUNT%
'INPUT I$
IPLOT = 0
ITER& = 0
SAVET = TIME
TIME = 0!
680 IF (JCOUNT% <= NOCYC%) THEN
GOTO 125 'loop2
END IF
700 ' CONTINUE
GOTO 10 'LOOP 1
1000 ' CONTINUE
' STOP
CALL closefilenr(filecount%, outfilenr%)
CLOSE outfilenr%
LOCATE 24, 1
PRINT " ######## END OF PROGRAM #########";
PRINT " OUTPUTFILE="; outfilename$; '" PRESS ANEY KEY ";
' INPUT I$
END
1188 'Hit key
IF INKEY$ <> "" THEN skipit = FALSE
IF skipit = FALSE THEN
PRINT
PRINT "#### Hit any key to continue. ##### (spacebar to skip most outputs)"
1189 I$ = INKEY$
IF I$ = "" THEN GOTO 1189
IF I$ = " " THEN skipit = TRUE
END IF
RETURN
1900 'initialize input vars
RESTORE 5000
FOR NBR% = 1 TO 20
READ vin$(NBR%), vit$(NBR%)
NEXT NBR%
RESTORE 4000
FOR NBR% = 1 TO 15 + 12 + 13 + 13 + 2
READ vi(NBR%)
NEXT NBR%
PRINT
PRINT "$STRLNG"
FOR NBR% = 1 TO 15
PRINT vin$(NBR%), " =", vi(NBR%) '; " ";
' PRINT vit$(NBR%)
NEXT NBR%
'PRINT vin$(19); " =",
INPUT I$
CLS
RETURN
1910 'print output file input vars from aray
PRINT #outfilenr%, "$STRLNG"
FOR NBR% = 1 TO 15
PRINT #outfilenr%, vin$(NBR%), " = "; vi(NBR%)
' PRINT vit$(NBR%)
NEXT NBR%
PRINT #outfilenr%, vin$(16),
FOR NBR% = 16 TO 19
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 20 TO 23
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 24 TO 27
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, vin$(17),
FOR NBR% = 28 TO 31
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 32 TO 35
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 36 TO 39
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
PRINT #outfilenr%, vi(40)
PRINT #outfilenr%, vin$(18),
FOR NBR% = 41 TO 44
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 45 TO 48
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 49 TO 52
PRINT #outfilenr%, vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
PRINT #outfilenr%, vi(53)
PRINT #outfilenr%, vin$(19), " = "; vi(54)
PRINT #outfilenr%, vin$(20), " = "; vi(55)
PRINT #outfilenr%,
'CALL closefilenr(filecount%, outfilenr%)
'CLOSE outfilenr%
RETURN
1920 'print output file, inputs from variables
PRINT #outfilenr%, "$STRLNG"
FOR NBR% = 1 TO 15
PRINT #outfilenr%, vin$(NBR%), " = "; ' vi(NBR%)
IF vin$(NBR%) = "P" THEN
PRINT #outfilenr%, P
ELSEIF vin$(NBR%) = "REALGS" THEN
PRINT #outfilenr%, REALGS
ELSEIF vin$(NBR%) = "FIPCYC" THEN
PRINT #outfilenr%, FIPCYC
ELSEIF vin$(NBR%) = "IPRINT" THEN
PRINT #outfilenr%, IPRINT
ELSEIF vin$(NBR%) = "ITMPS" THEN
PRINT #outfilenr%, ITMPS
ELSEIF vin$(NBR%) = "COEFF" THEN
PRINT #outfilenr%, COEFF
ELSEIF vin$(NBR%) = "T3" THEN
PRINT #outfilenr%, T3
ELSEIF vin$(NBR%) = "P3" THEN
PRINT #outfilenr%, P3
ELSEIF vin$(NBR%) = "OMEGA" THEN
PRINT #outfilenr%, OMEGA
ELSEIF vin$(NBR%) = "FACT1" THEN
PRINT #outfilenr%, FACT1
ELSEIF vin$(NBR%) = "FACT2" THEN
PRINT #outfilenr%, FACT2
ELSEIF vin$(NBR%) = "NOCYC%" THEN
PRINT #outfilenr%, NOCYC%
ELSEIF vin$(NBR%) = "NSTRT" THEN
PRINT #outfilenr%, NSTRT
ELSEIF vin$(NBR%) = "NOEND" THEN
PRINT #outfilenr%, NOEND
ELSEIF vin$(NBR%) = "MWGAS%" THEN
PRINT #outfilenr%, MWGAS%
END IF
NEXT NBR%
PRINT #outfilenr%, vin$(16),
FOR NBR% = 16 TO 19
PRINT #outfilenr%, TG(NBR% - 15), ')vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 20 TO 23
PRINT #outfilenr%, TG(NBR% - 15), 'vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 24 TO 27
PRINT #outfilenr%, TG(NBR% - 15), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, vin$(17),
FOR NBR% = 28 TO 31
PRINT #outfilenr%, TGA(NBR% - 27), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 32 TO 35
PRINT #outfilenr%, TGA(NBR% - 27), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 36 TO 39
PRINT #outfilenr%, TGA(NBR% - 27), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
'PRINT #outfilenr%, TGA(NBR% - 27)' + 1)'vi(40)
PRINT #outfilenr%, TGA(13)'vi(40)
PRINT #outfilenr%, vin$(18),
FOR NBR% = 41 TO 44
PRINT #outfilenr%, TM(NBR% - 40), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 45 TO 48
PRINT #outfilenr%, TM(NBR% - 40), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
FOR NBR% = 49 TO 52
PRINT #outfilenr%, TM(NBR% - 40), ' vi(NBR%),
NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, ,
'PRINT #outfilenr%, TM(NBR% - 40 + 1)'vi(53)
PRINT #outfilenr%, TM(13)'vi(53)
PRINT #outfilenr%, vin$(19), " = "; RHCFAC' vi(54)
PRINT #outfilenr%, vin$(20), " = "; JIP%'vi(55)
PRINT #outfilenr%,
RETURN
1950 'initialize output vars
RESTORE 5010
FOR NBR% = 1 TO 48
READ von$(NBR%), vot$(NBR%)
NEXT NBR%
RETURN
1960 'print outputs not used here, is in CYCLE
FOR NBR% = 1 TO 48
PRINT von$(NBR%); "="; vot$(NBR%); " "; ', STORE$(11,NBR%)
NEXT NBR%
FOR NBR% = 1 TO 44
PRINT #outfilenr%, von$(NBR%); "="; vot$(NBR%); " "; ', STORE$(11,NBR%)
NEXT NBR%
RETURN
2000 ' DATA FOR SUB CNDCT not used
' TR0 ,TR1 ,TR2 ,TCYL(1),TCYL(2),TCYL(3)
DATA 1760.,1172.,585.,1760. ,1255. ,975.
' RI2 ,R2 , R3 , R4 , XL3, XL4
DATA 1.38,1.62,1.56,1.53,1.126,0.400
' RI1 ,R0 , R1 , XL1 ,XL2
DATA 0.454,0.521,0.495,0.400,0.470
' ACONDD, DCONDD/
DATA 0.5875, 1.716
' TCAN(1),TCAN(2),TCAN(3)
DATA 1400. ,1280. ,1100.
' ACAN(1),ACAN(2),DCAN(1),DCAN(2)
DATA 1.154 ,1.196 ,0.5 ,0.75
' CLEARDIS,DISPD,STROKE/
DATA .010 ,2.75 ,1.23
' ACONDR,ACYL,ACONDD,ACAN---AREA S IN IN2
' DCONDR,DCYL,DCONDD,DCAN---DISTANCE IN IN
3000 ' DATA FOR SUB CYCL
DATA J%,K%,STORE,OLDTIM/2*0.947*0./
DATA JIP1/0/
DATA DH/1.437,3*0.1190,5*0.003800,3*0.0409,1.487/
DATA ACS/3.250,3*0.05561,5*0.4380,3*0.05124,3.150/
DATA AHT/0.0,2*5.720,3.922,5*82.56,3*2.288,0.0/
DATA CPM/4*10.0,5*.011, 4*10.0/
DATA XL/0.0,2*3.06,2.098,5*0.178,3*0.587,0.0/
DATA IDEX/1/
DATA GRAV/32.2/, PSI/0.006945/
4000 ' data input
' NAMELIST /STRLNG/P,REALGS,FIPCYC,IPRINT,ITMPS,COEFF,T3,P3,OMEGA,
'FACT1,FACT2,NOCYC%,NSTRT,NOEND,MWGAS%,TG,TGA,TM,RHCFAC,JIP%
' P ,REALGS,FIPCYC,IPRINT,ITMPS,COEFF,T3 ,P3 ,OMEGA
DATA .85e+3,.1e+1 , .1e4 , 500 ,0 ,.1e-3,.7e3,.82e3,.5e2
' FACT1,FACT2,NOCYC%,NSTRT,NOEND,MWGAS%
DATA .4 ,.1e2 ,50 ,5 ,20 ,4
'TG(12)
DATA .186e4, .186e4, .186e4, .186e4
DATA .1658e4, .1494e4, .1304e4, .1074e4
DATA .775e3, .54e3, .54e3, .54e3
'TGA(13)
DATA .19e4, .186e4, .186e4, .186e4
DATA .1658e4, .1494e4, .1304e4, .1074e4
DATA .775e3, .54e3, .54e3, .54e3
DATA .54e3
'TM(13)
DATA .176e4, .193e4, .193e4, .1625e4
DATA .1658e4, .1494e4, .1304e4, .1074e4
DATA .775e3, .54e3, .54e3, .54e3
DATA .54e3
'RHCFAC
DATA .1e1
'JIP%
DATA 1
DATA END
5000 ' input variables sorted as in input file
DATA "P"
DATA "initial pressure, lbf/inch^2 (N/m^2) (approximately equal to mean
pressure for GPU)"
DATA "REALGS"
DATA "1, use real-gas equation of state
0, use ideal-gas equation of state"
DATA "FIPCYC"
DATA "number of iterations per cycle"
DATA "IPRINT"
DATA "number of iterations between printouts (if JIP%=0)"
DATA "ITMPS"
DATA "1, printout gas and metal control volume temperature at every iteration
0, do not do the above"
DATA "COEFF"
DATA "leage coefficient for flow between buffer and compression spaces"
DATA "T3"
DATA "buffer space temperature (assumed constant), degrees R (Kelvin)"
DATA "P3"
DATA "initial buffer space pressure, lbf/inch^2 (N/m^2)"
DATA "OMEGA"
DATA "drive speed, Hz"
DATA "FACT1"
DATA "FACT1 and FACT2 are convergence constants used in regenerator
metal temperature convergence scheme"
DATA "FACT2"
DATA "FACT2 and FACT1 are convergence constants used in regenerator
metal temperature convergence scheme"
DATA "NOCYC%"
DATA "number of engine cycles to be run"
DATA "NSTRT"
DATA "number of cycle at which regenerator temperature convergence
scheme is turned on"
DATA "NOEND"
DATA"number of cycle at wich convergence scheme is turned on"
DATA "MWGAS%"
DATA "4, use helium working gas
not 4, use hydrogen working gas"
DATA "TG"
DATA "12 control volume interface gas temperatures, degrees R (Kelvin)"
DATA "TGA"
DATA "13 control volume gas temperatures, degrees R (Kelvin)"
DATA "TM"
DATA "13 metal boundary temperatures, degrees R (Kelvin)"
DATA "RHCFAC"
DATA "regenerator heat-transfer coefficient multiplication factor"
DATA "JIP%"
DATA ">0, short-form printout
=0, long-form printout"
5010 'output variables sorted as in output file
DATA "QIN"
DATA "heat into gas in expansion space and heater tubes over one cycle,
ft-lbf (J)"
DATA "QOUT"
DATA "heat out of gas in cooler tubes and compression space over
one cycle, ft-lbf (J)"
DATA "WRKEXP"
DATA "work done by the gas in compression space over one cycle, ft-lbf (J)"
DATA "WRKCMP"
DATA "work done by the gas in compression space over one cycle, ft-lbf (J)"
DATA "WORKTOT"
DATA "net work done by gas in working space over one cycle,
WRKEXP + WRKCMP, ft-lbf (J)"
DATA "EFFTOT"
DATA "efficiency, WRKTOT/QIN"
DATA "REFF1"
DATA "like REFF2 a measure of regenerator effectiveness"
DATA "REFF2"
DATA "like REFF1 a measure of regenerator effectiveness"
DATA "QREGEN"
DATA "net heat flow from gas to regenerator metal over one cycle,
ft-lbf/cycle (J/cycle)"
DATA "PWRHP"
DATA "indicated power -WRKTOT/cycle period, hp (kW)"
DATA "FREQ"
DATA "rotational frequency of crankshaft, Hz"
DATA "WRKLT"
DATA "work lost due to pressure drop over one cycle, ft-lbf (J)"
DATA "ENTH"
DATA "net enthalpy flow from working space to buffer space
over one cycle, ft-lbf (J)"
DATA "UTOTAL"
DATA "internal energy content of working-space gas, ft-lbf (J)"
DATA "EN3PM"
DATA "rate of displacement of working-space gas, inch^3/min (cm^3/min)"
DATA "AVGWSP"
DATA "mean working-space pressure, lbf/inch^2 (N/m^2)"
DATA "WDISP"
DATA "net work done on displacer piston by working-space gas, ft-lbf (J)"
DATA "WPIST"
DATA "work done on power piston by working-space gas ft-lbf (J)"
DATA "WTPIST"
DATA "?????"
DATA "WRKLP"
DATA "power-piston work lost due to pressure drop from regenerator
to compression space, ft-lbf/cycle J/cycle"
DATA "QEXPN"
DATA "heat from metal wall to gas in expansion space, ft-lbf/cycle (J/cycle)"
DATA "QEXPP"
DATA "heat from gas to metal wall in expansion space ft-lbf/cycle (J/cycle)"
DATA "QHEATN"
DATA "heat from heater to gas, ft-lbf/cycle (J/cycle)"
DATA "QHEATP"
DATA "heat from gas to heater, ft-lbf/cycle (J/cycle)"
DATA "QCOOLN"
DATA "heat from cooler to gas, ft-lbf/cycle (J/cycle)"
DATA "QCOOLP"
DATA "heat from gas to cooler, ft-lbf/cycle (J/cycle)"
DATA "QCOMPN"
DATA "heat from compression-space wall to gas, ft-lbf/cycle (J/cycle)"
DATA "QCOMPP"
DATA "heat from gas to compression-space wall, ft-lbf/cycle (J/cycle"
DATA "QREGN"
DATA "heat from regenerator metal to gas, ft-lbf/cycle (J/cycle)"
DATA "QREGP"
DATA "heat from gas to regenerator metal, ft-lbf/cycle (J/cycle)"
DATA "WDISPP"
DATA "work done on piston by gas, ft-lbf/cycle (J/cycle)"
DATA "WDISPN"
DATA "work done on gas by displacer, ft-lbf/cycle (J/cycle)"
DATA "WPISTP"
DATA "work done on power piston by gas, ft-lbf/cycle (J/cycle"
DATA "WPISTN"
DATA "work done on gas by piston, ft-lbf/cycle ft-lbf (J/cycle)"
DATA "QCNDRI"
DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)"
DATA "QCNDRO"
DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)"
DATA "QCNDCL"
DATA "heat conduction through cylinder, ft-lbf/cycle (J/cycle)"
DATA "QCNDD"
DATA "heat conduction through displacer, ft-lbf/cycle (J/cycle)"
DATA "QCNDCN"
DATA "heat conduction through insulation container, ft-lbf/cycle (J/cycle)"
DATA "TGEXPA"
DATA "cycle time average temperature in expansion space, degree R (kelvin)"
DATA "TGCMPA"
DATA "cycle-time average temperature in compression space, degrees R (kelvin)"
DATA "QSHTL"
DATA "shuttle heat loss, ft-lbf/cycle (J/cycle)"
DATA "QINB"
DATA "QINB= -qexp - qheatr"
DATA "WRKBAS"
DATA "WRKBAS = wrktot + wrklt"
DATA "TGCYC"
DATA "cycle time average temperature of TG= gas temperature at interfase
of two adjacent control volumes , degree R (kelvin)"
DATA "TGACYC"
DATA "cycle-time average temperature of TGA= working bulk gas temperature
of a control volume, degrees R (kelvin)"
DATA "TMCYC"
DATA "cycke-time average temperarure of TM= metal temperature of a metal
control volume, degree R (kelvin)"
'These are input arays?
DATA "TG"
DATA "gas temperature at interfase of two adjacent control volumes
degrees R (kelvin)"
DATA "TGA"
DATA "working bulk gas temperature of a control volume, degrees R (kelvin)"
DATA "TM"
DATA "metal temperature of a metal control volume, degree R (kelvin)"
DEFSNG I-N
SUB closefilenr (filecount%, filenr%)
'CALL closefilenr (filecount%,filenr%)
DEFINT I-N
IF filecount% AND 2 ^ filenr% = 0 THEN
PRINT "<< ERROR >> filenr"; filenr%; " is already closed."
PRINT "Hit the 'Enter' key to continue.";
INPUT I$
ELSE
CLOSE #filenr%
filecount% = filecount% - 2 ^ filenr%
END IF
END SUB
DEFSNG I-N
SUB CNDCT (TMA(), QCREG, QCYL, QCONDD, QCAN, QSHTTL, MWGAS%)
' SUBROUTINE CNDCT
DEFINT I-N
' DIM TMA(13), TCYL(3), TCAN(3)
DIM TCYL(3), TCAN(3)
DIM ACAN(2), DCAN(2)
' COMMON /CYC/ QEXP, DHEATR, QREGEN, QCOOLR, QCOMP, ENFRTH, ENFHTR
' COMMON ENFCTR, ENFRTC, WRKP, WRKEXP
' COMMON WRKCMP, WRKLEX, WRKLCM, WRKLP, WRKH, SUMNUM(), SUMDEN(), ENTH
' COMMON WDISP, WPIST, WTPIST, QEXPN,QEXPP, QHEATN, QHEATP
' COMMON QCOOLN, QCOOLP, QCOMPN, QCOMPP, QREGN, QREGP, WDISPP, WDISPN, WPISTP
' COMMON WPISTN, QCNDRI, QCNDRO, QCNDRCL, QCNDD, QCNDCN, QSHTL, TGEXPA, TGCMPA
SHARED QEXP, DHEATR, qregen, qcoolr, qcomp, ENFRTH, ENFHTR
SHARED enfctr, enfrtc, wrkp, wrkexp
SHARED wrkcmp, WRKLEX, wrklcm, wrklp, wrkh, SUMNUM(), SUMDEN(), enth
SHARED wdisp, wpist, wtpist, QEXPN, qexpp, qheatn, qheatp
SHARED qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
SHARED wpistn, qcndri, qcndro, QCNDRCL, qcndd, qcndcn, QSHTL, tgexpa, tgcmpa
' DATA TR0,TR1,TR2,TCYL/1760.,1172.,585.,1760.,1255.,975./,
' 1 RI2,R2,R3,R4,XL3,XL4/1.38,1.62,1.56,1.53,1.126,0.400/,
' 1 RI1,R0,R1,XL1,XL2,/0.454,0.521,0.495,0.400,0.470/,
' 1 ACONDD,DCONDD/0.5875,1.716/,
' 1 TCAN/1400.,1280.,1100./,
' 1 ACAN,DCAN/1.154,1.196,0.5,0.75/,
' 1 CLEARDIS,DISPD,STROKE/.010,2.75,1.23/
' ACONDR,ACYL,ACONDD,ACAN---AREA S IN IN2
' DCONDR,DCYL,DCONDD,DCAN---DISTANCE IN IN
'RESTORE 2000 ' DATA FOR SUB CNDCT
'READ TR0, TR1, TR2, TCYL(1), TCYL(2), TCYL(3)
'READ RI2, R2, R3, R4, XL3, XL4
'READ RI1, R0, R1, XL1, XL2
'READ ACONDD, DCONDD
'READ TCAN(1), TCAN(2), TCAN(3)
'READ ACAN(1), ACAN(2), DCAN(1), DCAN(2)
'READ CLEARDIS, DISPD, STROKE
TR0 = 1760!: TR1 = 1172!: TR2 = 585!
TCYL(1) = 1760!: TCYL(2) = 1255!: TCYL(3) = 975!
RI2 = 1.38: R2 = 1.62: R3 = 1.56: R4 = 1.53: XL3 = 1.126: XL4 = .4
RI1 = .454: R0 = .521: R1 = .495: XL1 = .4: XL2 = .47
ACONDD = .5875: DCONDD = 1.716
TCAN(1) = 1400!: TCAN(2) = 1280!: TCAN(3) = 1100!
ACAN(1) = 1.154: ACAN(2) = 1.196: DCAN(1) = .5: DCAN(2) = .75
CLEARDIS = .01: DISPD = 2.75: STROKE = 1.23
'CLEAR is reserved word in Quick Basic
' PAGE 59 (65 IN PDF) ------------------
' TCYL,TCAN--- TEMPERATURES IN DEGREES R
' CLEARDIS--CLEARANCE BETEEN DISPLACER AND CYLINDER, IN.
' DISPD---DISPLACER DIAMETER, IN
' STROKE---DISPLACER STROKE, IN.
' CALCULATE HEAT CONDUCTED THROUGH CYLINDER,BTU/SEC
TCAVG1 = (TCYL(1) + TCYL(2)) / 2!
TCAVG2 = (TCYL(2) + TCYL(3)) / 2!
THCND1 = (.00387 * TCAVG1 + 6.793) / (3600! * 12!)
IF (TCAVG1 > 1060!) THEN
THCND1 = (.00467 * TCAVG1 + 5.932) / (3600! * 12!)
END IF
THCND2 = (.00387 * TCAVG2 + 6.793) / (3600! * 12!)
IF (TCAVG2 > 1060!) THEN
THCND2 = (.00467 * TCAVG2 + 5.932) / (3600! * 12!)
END IF
QCYL1 = 6.283 * THCND1 * RI2 * (R3 - R2) * (TCYL(1) - TCYL(2)) / (XL3 * LOG(((R3 - RI2) * (R2 + RI2)) / ((R3 + RI2) * (R2 - RI2))))
QCYL2 = 6.283 * THCND2 * RI2 * (R4 - R3) * (TCYL(2) - TCYL(3)) / (XL4 * LOG(((R4 - RI2) * (R3 + RI2)) / ((R4 + RI2) * (R3 - RI2))))
QCYL = (QCYL1 + QCYL2) / 2!
' CALCULATE HEAT CONDUCTED THROUGH DISPLACER,BTU/SEC
TGAVG = (tgexpa + tgcmpa) / 2!
THCOND = (.00387 * TGAVG + 6.793) / (3600! * 12!)
IF (TGAVG > 1060!) THEN THCOND = (.00467 * TGAVG + 5.932) / (3600! * 12!)
QCONDD = ACONDD * THCOND * (tgexpa - tgcmpa) / DCONDD
' CALCULATE HEAT CONDUCTED THROUGH EXTERNAL CAN,BTU/SEC
QCAN = 0!
' DO 300 I=1,2
FOR I = 1 TO 2
TCAVG = (TCAN(I) + TCAN(I + 1)) / 2!
THCOND = (.00387 * TCAVG + 6.793) / (3600! * 12!)
IF (TCAVG > 1060!) THEN THCOND = (.00467 * TCAVG + 5.932) / (3600! * 12!)
' 300 QCAN=QCAN+ACAN(I)*THCOND*(TCAN(I)-TCAN(I+1))/(DCAN(I)*2.)
QCAN = QCAN + ACAN(I) * THCOND * (TCAN(I) - TCAN(I + 1)) / (DCAN(I) * 2!)
NEXT I
' CALCULATE HEAT CONDUCTION THROUGH REGENERATOR,BTU/SEC
TRAVG1 = (TR0 + TR1) / 2!
TRAVG2 = (TR1 + TR2) / 2!
THCND1 = (.00387 * TRAVG1 + 6.793) / (3600! * 12!)
IF (TRAVG1) > (1060!) THEN
THCND1 = (.00467 * TRAVG1 + 5.932) / (3600! * 12!)
END IF
THCND2 = (.00387 * TRAVG2 + 6.793) / (3600! * 12!)
IF (TRAVG2) > (1060!) THEN
THCND2 = (.00467 * TRAVG2 + 5.932) / (3600! * 12!)
END IF
QCRIN = 3.142 * (R0 ^ 2 - RI1 ^ 2) * THCND1 * (TR0 - TR1) / XL1
QCROUT = 6.283 * THCND2 * RI1 * (R1 - R0) * (TR0 - TR1) / (XL2 * LOG(((R1 - RI1) * (R0 + RI1)) / ((R1 + RI1) * (R0 - RI1))))
QCREG = (QCRIN + QCROUT) / 2! * 8!
' CALCULATE SHUTTLE HEAT TRANSFER, BTU/SEC
THCNDG = 2.481E-09 * (TCYL(1) + TCYL(3)) / 2! + 1.263E-06
IF (MWGAS%) = (4) THEN
THCNDG = 1.962E-09 * (TCYL(1) + TCYL(3)) / 2! + 1.115E-06
END IF
QSHTTL = THCNDG * 3.1416 * DISPD * STROKE ^ 2 * (TCYL(1) - TCYL(3)) / (8! * CLEARDIS * (XL3 + XL4))
' RETURN
' END
END SUB
DEFSNG I-N
SUB CYCL (TIME, UTOTAL, STROKV, AVGWSP)
' SUBROUTINE CYCL
'
' THIS SUBROUTINE IS CALLED JUST ONCE PER CYCLE--WHEN THE CYCLE IS COMLETE.
' IT CALCULATES--NET HEAT IN AND OUT, NET WORK PER CYCLE, INDICATED POWERS
' AND EFFICIENCY, ETC. HEAT AND WORK QANTITIES ARE STORED FOR CALCULATION
' QUANTITIES ARE RESET TO ZERO IN PREPARATION FOR NEXT CYCLE.
DEFINT I-N
SHARED VAR(), STORE()
' COMMON DTIME, F(), WG(), DELP(), P, PE, PC
SHARED DTIME, F(), WG(), DELP(), P, PE, PC
SHARED outfilenr%
'DATA J%,K%,STORE,OLDTIM/2*0,947*0./
'DATA JIP1/0/
SHARED von$(), vot$() 'variables outputs names and explanation text strings
STATIC J%, K%, OLDTIME, JIP1
'J% = 0
'K% = 0
'FOR LCOUNT% = 1 TO 11
' FOR KCOUNT% = 1 TO 86
' STORE(LCOUNT%, KCOUNT%) = 0!
' NEXT KCOUNT%
'NEXT LCOUNT%
'OLDTIME = 0!
'JIP1 = 0
' COMMON /CYC/ QEXP, QHEATR, QREGEN, QCOOLR, QCOMP, ENFRTH, ENFHTR
' COMMON ENFCTR, ENFRTC, WRKP, WRKEXP
' COMMON WRKCMP, WRKLEX, WRKLCM, WRKLP, WRKH, SUMNUM(), SUMDEN(), ENTH
' COMMON WDISP, WPIST, WTPIST, QEXPN, QEXPP, QHEATN, QHEATP
' COMMON QCOOLN, QCOOLP, QCOMPN, QCOMPP, QREGN, QREGP, WDISPP, WDISPN, WPISTP
' COMMON WPISTN, QCNDRI, QCNDRO, QCNDCL, QCNDD, QCNDCN, QSHTL, TGEXPA, TGCMPA
' COMMON /TAVCYC/ TGCYC(), TGACYC(), TMCYC()
' COMMON /JNDEX/ JIP%, NOCYC%
SHARED QEXP, qheatr, qregen, qcoolr, qcomp, ENFRTH, ENFHTR
SHARED enfctr, enfrtc, wrkp, wrkexp
SHARED wrkcmp, WRKLEX, wrklcm, wrklp, wrkh, SUMNUM(), SUMDEN(), enth
SHARED wdisp, wpist, wtpist, QEXPN, qexpp, qheatn, qheatp
SHARED qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
SHARED wpistn, qcndri, qcndro, qcndcl, qcndd, qcndcn, QSHTL, tgexpa, tgcmpa
SHARED TGCYC(), TGACYC(), TMCYC()
SHARED JIP%, NOCYC%
' DIM TGCYC(13), TGACYC(13), TMCYC(13)
' EQUIVALENCE (VAR(1),QEXP)
VAR(1) = QEXP ' ?????????????????????????? IS THIS RIGHT ?
J% = J% + 1
K% = K% + 1
JIP1 = JIP1 + 1
IF (OLDTIM) > (TIME) THEN OLDTIME = 0!
DELTIM = TIME - OLDTIM
OLDTIM = TIME
QIN = -QEXP - qheatr + qcndri + qcndcl + qcndd + qcndcn + QSHTL - wrklt / 2!
QINB = -QEXP - qheatr
qout = qcoolr + qcomp + qcndro + qcndcl + qcndd + qcndcn + QSHTL
wrktot = wrkexp + wrkcmp
wrklt = WRKLEX + wrklcm
WRKBAS = wrktot + wrklt
efftot = wrktot / QIN
EFFP = wrkp / QIN
PRINT "ENFHTR="; ENFHTR
PRINT "ENFRTH="; ENFRTH
PRINT "F(4)="; F(4)
' reff1 = ENFRTH / ENFHTR !!!!!!!!!!!!!!!! BLOCKED !!!!!!!!!!!!!
' reff2 = enfctr / enfrtc
pwrhp = wrkh / (550! * DELTIM)
' PAGE 57 (63 IN PDF) ----------------------------------------
PWRKW = pwrhp * .7457
freq = 1! / DELTIM
EN3PM = STROKV * freq * 60!
' STORE VARIABLES FOR CALCULATION OF AVERAGES OVER 5 CYCLES.
STORE(K%, 1) = QIN
STORE(K%, 2) = qout
STORE(K%, 3) = wrkexp
STORE(K%, 4) = wrkcmp
STORE(K%, 5) = wrktot
STORE(K%, 6) = efftot
STORE(K%, 7) = reff1
STORE(K%, 8) = reff2
STORE(K%, 9) = qregen
STORE(K%, 10) = pwrhp
STORE(K%, 11) = freq
STORE(K%, 12) = wrklt
STORE(K%, 13) = enth
STORE(K%, 14) = UTOTAL
STORE(K%, 15) = EN3PM
STORE(K%, 16) = AVGWSP
STORE(K%, 17) = wdisp
STORE(K%, 18) = wpist
STORE(K%, 19) = wtpist
STORE(K%, 20) = wrklp
STORE(K%, 21) = QEXPN
STORE(K%, 22) = qexpp
STORE(K%, 23) = qheatn
STORE(K%, 24) = qheatp
STORE(K%, 25) = qcooln
STORE(K%, 26) = qcoolp
STORE(K%, 27) = qcompn
STORE(K%, 28) = qcompp
STORE(K%, 29) = qregn
STORE(K%, 30) = qregp
STORE(K%, 31) = wdispp
STORE(K%, 32) = wdispn
STORE(K%, 33) = wpistp
STORE(K%, 34) = wpistn
STORE(K%, 35) = qcndri
STORE(K%, 36) = qcndro
STORE(K%, 37) = qcndcl
STORE(K%, 38) = qcndd
STORE(K%, 39) = qcndcn
STORE(K%, 40) = tgexpa
STORE(K%, 41) = tgcmpa
STORE(K%, 42) = QSHTL
STORE(K%, 43) = QINB
STORE(K%, 44) = WRKBAS
' DO 90 L=1,13
FOR L = 1 TO 13
IF (L) = (13) THEN GOTO 91
M = L + 48
STORE(K%, M) = TGCYC(L)
91 M = L + 60
STORE(K%, M) = TGACYC(L)
M = L + 73
90 STORE(K%, M) = TMCYC(L)
NEXT L
' 100 ' CONTINUE
IF (JIP% > 0) AND (JIP1 <> NOCYC%) THEN
'PRINT "(JIP% > 0) AND (JIP1 <> NOCYC%)"
'INPUT I$
GOTO 301
END IF
IF (JIP% > 0) THEN
PRINT "JIP%>0 WRITE (6,202) LAST CYCLE"
202 'FORMAT ('LAST CYCLE')
PRINT
PRINT "LAST CYCLE"
PRINT
PRINT #outfilenr%,
PRINT #outfilenr%, "LAST CYCLE"
PRINT #outfilenr%,
END IF
' WRITE (6,200)
GOSUB 2200 '200 FORMAT IS IN GOSUB 2200
' PAGE 58 (64 IN PDF) ----------------------------------
'200 'FORMAT (/' QIN=',F8.3.1X,'QOUT=',F8.3,1X,'WRKEXP='F8.3,1X,
'1'WRKCMP=',F8.3,1X,'WRKTOT=',F8.3,1X,
'2'EFFTOT=',F8.3,1X,'REFF1=',F8.3,1X,'REFF2=',F8.3,1X/
'3 ' QREGEN=',F8.3,' PWRHP=',F8.3,' FREQ=',F8.3,
'4 ' WRKLT=',F8.3, 'ENTH=',F8.3,' UTOTAL=',
' 5F10.3/' EN3PN=',F10.3,' AVGWSP=',F9.3,
' 6' WDISP=',F8.3,' WPIST=',F8.3,' WTPIST=',F8.3,' WRKLP='F8.3/,
'7 ' QEXPN=',F8.3,' QEXPP=',F8.3,' QHEATN=',F8.3,' QHEATP='F8.3,
' 8' QCOOLN=',F8.3,' QCOOLP=',F8.3,' QCONPN=',F8.3,' QCOMPP=',F8.3/,
'9 ' QREGN=',F9.3,' QREGP=',F8.3,' WDISPP=',F8.3,' WDISPN=',F8.3,
' 1' WPISTP=',F8.3,' WPISTN=',F8.3,' QCNDRI=',F7.3,' QCNDRO=',F7.3,/
' 2' QCNDCL=',F8.3,' QCNDD=',F8.3,' QCNDCN=',F8.3,' TGEXPA=',F8.3,
' 3' TGCMPA=',F8.3,' QSHTL=',F8.3/,'QINB =',F8.3,' WRKBAS ='F8.3)
301 ' CONTINUE
' DO 300 I=1,66
FOR I = 1 TO 66
' 300 VAR(I)=0.0)
VAR(I) = 0!
NEXT I
IF (J% = 5) THEN GOTO 400
GOTO 601' RETURN
400 K% = 0
J% = 0
' DO 450 KK=1,86
FOR KK = 1 TO 86
STORE(11, KK) = STORE(1, KK) / 5!
' DO 450 IK=2,5
FOR IK = 2 TO 5
450 STORE(11, KK) = STORE(11, KK) + STORE(IK, KK) / 5!
NEXT IK
NEXT KK
575 ' CONTINUE
IF (JIP% > 0) AND (JIP1 <> NOCYC%) THEN GOTO 601
' WRITE (6,600)
' WRITE ' (6,200) (STORE(11,KK),KK=1,44)
'1960 'print outputs OF LAST 5 CYCLES
PRINT "AVERAGE VALUE OF LAST 5 CYCLES"
FOR NBR% = 1 TO 48
' PRINT von$(NBR%); "="; STORE(11,NBR%), vot$(NBR%)
PRINT von$(NBR%); "="; STORE(11, NBR%)
NEXT NBR%
'PRINT " TGCYC()", , " TGACYC()", " TMCYC() "
'FOR NBR% = 1 TO 13
' PRINT NBR%;
' IF NBR% < 10 THEN PRINT " ";
' IF NBR% <= 12 THEN
' M = NBR% + 48
' PRINT STORE(1, M),
' ELSE
' PRINT , ,
' END IF
' M = NBR% + 60
' PRINT STORE(1, M),
' M = NBR% + 73
' PRINT STORE(1, M)
'NEXT NBR%
PRINT #outfilenr%,
PRINT #outfilenr%, "AVERAGE VALUE OF LAST 5 CYCLES"
FOR NBR% = 1 TO 47
IF NBR% < 45 THEN
PRINT #outfilenr%, von$(NBR%); "="; STORE(11, NBR%),
ELSE
IF (NBR% = 45) OR (NBR% = 48) THEN
PRINT #outfilenr%, von$(NBR%); "(12)", ,
ELSE
PRINT #outfilenr%, von$(NBR%); "(13)", ,
END IF
END IF
IF INSTR(vot$(NBR%), "
") THEN
PRINT #outfilenr%, LEFT$(vot$(NBR%), INSTR(vot$(NBR%), "
") - 1)
PRINT #outfilenr%, , , RIGHT$(vot$(NBR%), LEN(vot$(NBR%)) - 4 - INSTR(vot$(NBR%), "
"))
ELSE
PRINT #outfilenr%, vot$(NBR%)
END IF
' PRINT #outfilenr%, von$(NBR%); "="; STORE(1, NBR%)
NEXT NBR%
'PRINT " TGCYC()", " TGACYC()", " TMCYC() "
'FOR NBR% = 1 TO 13
' PRINT #outfilenr%, NBR%;
' IF NBR% < 10 THEN PRINT #outfilenr%, " ";
' IF NBR% <= 12 THEN
' M = NBR% + 48
' PRINT #outfilenr%, STORE(1, M),
' ELSE
' PRINT #outfilenr%, , ,
' END IF
' M = NBR% + 60
' PRINT #outfilenr%, STORE(1, M),
' M = NBR% + 73
' PRINT #outfilenr%, STORE(1, M)
'NEXT NBR%
'GOSUB 2200 ' PRINT FORMAT 200
IF (JIP%) > (0) THEN 'WRITE (6,204)
204 'FORMAT (' AVERAGE TEMPERATURES OVER LAST CYCLE')
PRINT
PRINT " AVERAGE TEMPERATURES OVER LAST CYCLE"
PRINT
PRINT #outfilenr%,
PRINT #outfilenr%, " AVERAGE TEMPERATURES OVER LAST CYCLE"
PRINT #outfilenr%,
END IF
'write (6,201) TGCYC,TGACYC,TMCYC
201 'FORMAT (' TGCYC=',4X,12F8.0/,' TGACYC=',13F8.0/,' TMCYC='13F8.0
' X/)
PRINT " TGCYC()"; , , ; " TGACYC()"; , ; " TMCYC() "
FOR NBR% = 1 TO 13
IF NBR% < 10 THEN PRINT " ";
IF NBR% <> 13 THEN
PRINT NBR%; " "; TGCYC(NBR%),
ELSE
PRINT NBR%; " "; , ,
END IF
PRINT TGACYC(NBR%), TMCYC(NBR%)
NEXT NBR%
PRINT #outfilenr%, " TGCYC()"; , , ; " TGACYC() "; , ; " TMCYC()"
FOR NBR% = 1 TO 13
IF NBR% < 10 THEN PRINT #outfilenr%, " ";
IF NBR% <> 13 THEN
PRINT #outfilenr%, NBR%; " "; TGCYC(NBR%),
ELSE
PRINT #outfilenr%, NBR%; " "; , ,
END IF
PRINT #outfilenr%, TGACYC(NBR%), TMCYC(NBR%)
NEXT NBR%
600 'FORMAT (1HD,'AVERAGE VALUES OF LAST 5 CYCLES')
601 ' RETURN
' END
GOTO 2500'END SUB CYCL
2200 ' PRINT FORMAT 200
' PRINT ' (6,200) QIN,QOUT,WRKEXP,WRKCMP,WRKTOT,EFFTOT,REFF1,REFF2,
PRINT "QIN ="; QIN; " QOUT ="; qout; " WRKEXP="; wrkexp;
PRINT " WRKCMP="; wrkcmp; " WRKTOT = "; wrktot; " EFFTOT = "; efftot;
PRINT " REFF1="; reff1; " REFF2="; reff2
PRINT "QREGEN="; qregen; " PWRHP="; pwrhp; " FREQ="; freq;
PRINT " WRKLT="; wrklt; " ENTH="; enth; " UTOTAL="; UTOTAL
PRINT "EN3PN="; EN3PN; " AVGWSP="; AVGWSP; " WDISP="; wdisp;
PRINT " WPIST="; wpist; " WTPIST="; wtpist; " WRKLP="; wrklp
PRINT "QEXPN="; QEXPN; " QEXPP="; qexpp; " QHEATN="; qheatn;
PRINT " QHEATP="; qheatp; " QCOOLN="; qcooln; " QCOOLP="; qcoolp;
PRINT " QCOMPN="; qcompn; " QCOMPP="; qcompp
PRINT "QREGN="; qregn; " QREGP="; qregp; " WDISPP="; wdispp;
PRINT " WDISPN="; wdispn; " WPISTP="; wpistp; " WPISTN="; wpistn;
PRINT " QCNDRI="; qcndri; " QCNDRO = "; qcndro
PRINT "QCNDCL="; qcndcl; " QCNDD="; qcndd; " QCNDCN="; qcndcn;
PRINT " TGEXPA="; tgexpa; " TGCMPA="; tgcmpa; " QSHTL="; QSHTL;
PRINT " QINB="; QINB; " WRKBAS="; WRKBAS
'print to output file
'PRINT #outfilenr%,
'PRINT #outfilenr%, "LAST CYCLE"
'PRINT #outfilenr%,
PRINT #outfilenr%, "QIN ="; QIN; " QOUT ="; qout; " WRKEXP="; wrkexp;
PRINT #outfilenr%, " WRKCMP="; wrkcmp; " WRKTOT = "; wrktot; " EFFTOT = "; efftot;
PRINT #outfilenr%, " REFF1="; reff1; " REFF2="; reff2
PRINT #outfilenr%, "QREGEN="; qregen; " PWRHP="; pwrhp; " FREQ="; freq;
PRINT #outfilenr%, " WRKLT="; wrklt; " ENTH="; enth; " UTOTAL="; UTOTAL
PRINT #outfilenr%, "EN3PN="; EN3PN; " AVGWSP="; AVGWSP; " WDISP="; wdisp;
PRINT #outfilenr%, " WPIST="; wpist; " WTPIST="; wtpist; " WRKLP="; wrklp
PRINT #outfilenr%, "QEXPN="; QEXPN; " QEXPP="; qexpp; " QHEATN="; qheatn;
PRINT #outfilenr%, " QHEATP="; qheatp; " QCOOLN="; qcooln; " QCOOLP="; qcoolp;
PRINT #outfilenr%, " QCOMPN="; qcompn; " QCOMPP="; qcompp
PRINT #outfilenr%, "QREGN="; qregn; " QREGP="; qregp; " WDISPP="; wdispp;
PRINT #outfilenr%, " WDISPN="; wdispn; " WPISTP="; wpistp; " WPISTN="; wpistn;
PRINT #outfilenr%, " QCNDRI="; qcndri; " QCNDRO = "; qcndro
PRINT #outfilenr%, "QCNDCL="; qcndcl; " QCNDD="; qcndd; " QCNDCN="; qcndcn;
PRINT #outfilenr%, " TGEXPA="; tgexpa; " TGCMPA="; tgcmpa; " QSHTL="; QSHTL;
PRINT #outfilenr%, " QINB="; QINB; " WRKBAS="; WRKBAS
RETURN
2500 'END SUB CYCL
END SUB
DEFSNG I-N
SUB dosfilename (filename$, mes$)
DEFINT I-N
mes$ = ""
IF filename$ = "" THEN
mes$ = "NO FILENAME GIVEN, FILE NAME NOW IS THE DEFAULD ONE"
EXIT SUB
END IF
WHILE LEFT$(filename$, 1) = " "
filename$ = RIGHT$(filename$, LEN(filename$) - 1)
WEND
WHILE RIGHT$(filename$, 1) = " "
filename$ = LEFT$(filename$, LEN(filename$) - 1)
WEND
IF INSTR(filename$, ".") = 0 THEN
IF LEN(filename$) > 8 THEN
mes$ = "FILE NAME IS TO LONG"
EXIT SUB
END IF
FOR ptr = 1 TO LEN(filename$)
IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
mes$ = mes$ + " IN FILE NAME"
EXIT SUB
END IF
END IF
NEXT ptr
ELSE
IF INSTR(filename$, ".") > 9 THEN
mes$ = "FILE NAME IS TO LONG"
EXIT SUB
ELSEIF INSTR(filename$, ".") = 1 THEN
mes$ = "WRONG FILE NAME, STARTING WITH CHARACTER '.'"
EXIT SUB
ELSEIF INSTR(filename$, ".") = LEN(filename$) THEN
mes$ = "WRONG FILE NAME, ENDING WITH CHARACTER '.'"
EXIT SUB
ELSEIF LEN(filename$) - INSTR(filename$, ".") > 3 THEN
mes$ = "FILE NAME EXTESION IS TO LONG."
EXIT SUB
END IF
FOR ptr = 1 TO INSTR(filename$, ".") - 1
IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
mes$ = mes$ + " IN FILE NAME."
EXIT SUB
END IF
END IF
NEXT ptr
FOR ptr = INSTR(filename$, ".") + 1 TO LEN(filename$)
IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
mes$ = mes$ + " IN FILENAME EXTENSION."
EXIT SUB
END IF
END IF
NEXT ptr
END IF
END SUB
DEFSNG I-N
SUB existfile (filename$, filecount%, mes$)
'call existfile(filename$,filecount%,mes$)
DEFINT I-N
mes$ = ""
CALL givefilenr(filecount%, filenr%)
OPEN filename$ FOR APPEND AS filenr%
fileend = SEEK(filenr%)
CLOSE filenr%
IF fileend = 1 THEN
KILL filename$
mes$ = "file " + filename$ + " : does not exist."
END IF
END SUB
DEFSNG I-N
SUB givefilenr (filecount%, filenr%)
'call givefilenr (filecount%, filenr%)
DEFINT I-N
IF filecount% = 0 THEN filecount% = 1
filenr% = 0
WHILE (filecount% OR 2 ^ filenr%) = filecount%
filenr% = filenr% + 1
WEND
filecount% = filecount% + 2 ^ filenr%
END SUB
DEFSNG I-N
SUB HEATX (TIME, TG(), TGA(), TMA(), REYNO(), QR(), X(), RODY, JCYCLE%, diamd, X1MIN, RODYMX, V(), SIG3, W, MWGAS%, RHCFAC, REALGS, CONDID, NOCYC%)
' CONDID IS NOT USED ???????????????????????????/
' SUBROUTINE HEATEX ---------
' THIS SUBROUTINE CALCALATES HEAT TRANSFER AND UPDATES TEMPERATURES AND
' PRESSURE
' COMMON DTIME, F(), WG(), DELP(), P, PE, PC
' DIM F(12), WG(13), DELP(13) ', P, PE, PC
' DIM x(2), HA(13), SUMDEN(13), DH(12)
' COMMON /CYC/ QECP, QHEATR, QREGEN, QCOOLR, QCOMP, ENFRTH, ENFHTR
' COMMON ENFCTR, ENFRTC, WRKP, WRKEXP
' COMMON WRKCMP, WRKLEX, WRKLCM, WRKLP, WRKH, SUMNUM(), SUMDEN(), ENTH
' COMMON WDISP, WPIST, WTPIST, QEXPN, QEXPP, QHEATN, QHEATP
' COMMON QCOOLN, QCOOLP, QCOMPN, QCOMPP, QREGN, QREGP, WDISPP, WDISPN, WPISTP
' COMMON WPISTN, QCNDRI, QCNDRO, QCNDCL, QCNDD, QCNDCN, QSHTL, TGEXPA, TGCMPA
' COMMON /DELTAP/ AD, AP, HA()
DEFINT I-N
SHARED DTIME, F(), WG(), DELP(), P, PE, PC
SHARED QECP, qheatr, qregen, qcoolr, qcomp, ENFRTH, ENFHTR
SHARED enfctr, enfrtc, wrkp, wrkexp
SHARED wrkcmp, WRKLEX, wrklcm, wrklp, wrkh, SUMNUM(), SUMDEN(), enth
SHARED wdisp, wpist, wtpist, QEXPN, qexpp, qheatn, qheatp
SHARED qcooln, qcoolp, qcompn, qcompp, qregn, qregp, wdispp, wdispn, wpistp
SHARED wpistn, qcndri, qcndro, qcndcl, qcndd, qcndcn, QSHTL, tgexpa, tgcmpa
SHARED AD, AP, HA()
'MOT TG(), TGA(), TMA(), reyno(), qr(), X(), V()
DIM COEFX(40), DPX(40)
' DIM TGA(13), TMA(13), TGAO(13), TMAO(13), X(2)
DIM TGAO(13), TMAO(13)
DIM CPM(13), Q(13), TMIX(13), DH(13), ACS(13)
' DIM AHT(13), COND(13), VIS(13),REYNO(13), qr(5)
DIM AHT(13), COND(13), VIS(13)
' DIM V(13), WGOLD(13), XL(13), FAVG(13), DNSTY(13), FRICT(13)
DIM WGOLD(13), XL(13), FAVG(13), DNSTY(13), FRICT(13)
DH(1) = 1.437
DH(2) = .119
DH(3) = .119
DH(4) = .119
DH(5) = .0038
DH(6) = .0038
DH(7) = .0038
DH(8) = .0038
DH(9) = .0038
DH(10) = .0409
DH(11) = .0409
DH(12) = .0409
DH(13) = 1.487
ACS(1) = 3.25
ACS(2) = .05561
ACS(3) = .05561
ACS(4) = .05561
ACS(5) = .438
ACS(6) = .438
ACS(7) = .438
ACS(8) = .438
ACS(9) = .438
ACS(10) = .05124
ACS(11) = .05124
ACS(12) = .05124
ACS(13) = 3.15
AHT(1) = 0!
AHT(2) = 5.72
AHT(3) = 5.72
AHT(4) = 3.922
AHT(5) = 82.56
AHT(6) = 82.56
AHT(7) = 82.56
AHT(8) = 82.56
AHT(9) = 82.56
AHT(10) = 2.288
AHT(11) = 2.288
AHT(12) = 2.288
AHT(13) = 0!
CPM(1) = 10!
CPM(2) = 10!
CPM(3) = 10!
CPM(4) = 10!
CPM(5) = .011
CPM(6) = .011
CPM(7) = .011
CPM(8) = .011
CPM(9) = .011
CPM(10) = 10!
CPM(11) = 10!
CPM(12) = 10!
CPM(13) = 10!
XL(1) = 0!
XL(2) = 3.06
XL(3) = 3.06
XL(4) = 2.098
XL(5) = .178
XL(6) = .178
XL(7) = .178
XL(8) = .178
XL(9) = .178
XL(10) = .587
XL(11) = .587
XL(12) = .587
XL(13) = 0!
' IDEX=1 ???????????????? not used could be mdex
' JH(I)=HYDRAULIC DIAMETER,IN.; ACS(I)=FLOW AREA, IN2; AHT=HEAT TRANSFER,IN2;
' CPM(I)=METAL HEAT CAPACITY, BTU/LBM-DEG,R; XL(I)=LENGTHS FOR CALCULATING
' PRESSURE DROP,IN
' CP,CV ARE IN BTU/LBM-DEG.R;GAS CONSTANT, R, IS
' IN IN-LBF/(LBM-DEG R)
B = .02358
CP = 3.484
GAMMA = 1.394
' PAGE 52 (58 IN MANUAL PDF)---------------------
R = 9197!
IF (MWGAS%) <> (4) THEN GOTO 5
B = .01613
CP = 1.239
R = 4632!
GAMMA = 1.667
5 ' CONTINUE
CV = CP / GAMMA
' DO 10 I=1,13
FOR I = 1 TO 13
IF (MWGAS%) = (4) THEN GOTO 6
' CONDUCTIVITY, COND(I) HAS UNITS BTU/IN-SEC-DEG.R
' VISCOSITY,VIS(I) HAS UNITS LBM/IN-SEC
COND(I) = 2.481E-09 * TGA(I) + 1.263E-06
VIS(I) = 4.683E-10 * TGA(I) + 2.619E-07
GOTO 8
6 COND(I) = 1.962E-09 * TGA(I) + 1.115E-06
VIS(I) = 1.058E-09 * TGA(I) + 5.919E-07
8 TMAO(I) = TMA(I)
' 10 TGAO(I)=TGA(I)
TGAO(I) = TGA(I)
NEXT I
POLD = P
' CALCULATE PRESSURE,LBF/IN2
SUMWT = 0!
SUMV = 0!
' DO 14 I=1,13
FOR I = 1 TO 13
SUMWT = SUMWT + WG(I) * TGA(I)
14 SUMV = SUMV + V(I)
NEXT I
P = R * SUMWT / (SUMV - B * R * W * (1! - SIG3) * REALGS)
FGAMMA = (P / POLD) ^ ((GAMMA - 1!) / GAMMA)
' DO 20 J%=1,13
FOR J% = 1 TO 13
TGAO(J%) = TGAO(J%) * FGAMMA
20 TMIX(J%) = TGAO(J%)
NEXT J%
' CALCULATE MASSES, WG(I), IN EACH GAS LUMP. ,LBM
' 2.20462234 LBM= 1 Kg (Harrie)
SUM = 0!
' DO 200 I=1,13
FOR I = 1 TO 13
'200
SUM = SUM + P * V(I) / (R * (TMIX(I) + B * P * REALGS))
NEXT I
' DO 220 I=1,13
FOR I = 1 TO 13
220 WGOLD(I) = WG(I)
NEXT I
' DO 240 I=1,13
FOR I = 1 TO 13
' 240 WG(I)=W*(1.-SIG3)*P*V(I)/(R*SUM*(TMIX(I)+B*P*REALGS))
WG(I) = W * (1! - SIG3) * P * V(I) / (R * SUM * (TMIX(I) + B * P * REALGS))
IF WG(I) = 0 THEN
PRINT "ERROE WG("; I; ")="; WG(I)
PRINT "ERROE V("; I; ")="; V(I)
INPUT I$
END IF
NEXT I
IF (MDEX) = (2) THEN GOTO 243
' DO 242 I=1,13
FOR I = 1 TO 13
242 WGOLD(I) = WG(I)
NEXT I
' CALCULATES FLOW RATES, F(I) IN LBM/SEC, BETWEEN GAS LUMPS.
243 MDEX = 2
F(1) = (WGOLD(1) - WG(1)) / DTIME
' DO 300 I=2,12
FOR I = 2 TO 12
300 F(I) = (WGOLD(I) - WG(I)) / DTIME + F(I - 1)
NEXT I
' DO 24 I=1,12
FOR I = 1 TO 12
IF (F(I)) < (0!) THEN GOTO 22
TG(I) = TGAO(I)
IF (I > 4) AND (I < 10) THEN
TG(I) = TGAO(I) - (TMA(5) - TMA(9)) / 8! + F(I) * DTIME * (TMA(5) - TMA(9)) / 8! / WGOLD(I)
END IF
GOTO 24
22 TG(I) = TGAO(I + 1)
IF (I > 3) AND (I < 9) THEN
TG(I) = TGAO(I + 1) + (TMA(5) - TMA(9)) / 8! + F(I) * DTIME * (TMA(5) - TMA(9)) / 8! / WGOLD(I + 1)
END IF
24 'CONTINUE
NEXT I
' CALCULATE AVERAGE FLOW RATES FOR EACH GAS LUMP.
FAVG(1) = F(1)
' DO 30 I=2,12
FOR I = 2 TO 12
' 30 FAVG(I)=(F(I-1)+F(I))/2.
FAVG(I) = (F(I - 1) + F(I)) / 2!
NEXT I
' FAVG(13) = F(13)
' F(12) IS IN DIM AND COMMON !!!!!!!!!!!!! F(13) DOES NOT EXIST
FAVG(13) = F(12) ' REPARATION !!!!!!!!!!!!!!!!!!
' REVISE GAS TEMPERATURES TO ACCOUNT FOR MIXING.
IF (F(1)) < (0!) THEN
TMIX(1) = (WGOLD(1) * TMIX(1) - F(1) * DTIME * TG(1)) / WG(1)
END IF
' DO 100 I=2,12
FOR I = 2 TO 12
' PAGE 53 (59 IN MANUAL PDF)---------------------------
' 100 TMIX(I)=(WGOLD(I)*TMIX(I)+F(I-1)*TG(I-1)-F(I)*TG(I))*
TMIX(I) = (WGOLD(I) * TMIX(I) + (F(I - 1) * TG(I - 1) - F(I) * TG(I)) * DTIME) / WG(I)
NEXT I
IF (F(12)) > (0!) THEN
TMIX(13) = (WGOLD(13) * TMIX(13) + F(12) * DTIME * TG(12)) / WG(13)
END IF
110 ' CONTINUE
' CALCULATE REYNOLDS NO S, REYNO(I); HEAT TRANSFER AREAS, AHT(I);
' HEAT TRANSFER COEFFICIENTS, HA(I) BETWEEN GAS + METAL LUMPS.
' HA(I) HAS UNITS BTU/SEC-DEG.R
REYNO(1) = DH(1) * ABS(F(1)) / (VIS(1) * ACS(1)) + .0001
AHT(1) = 3.1416 * diamd * (X(1) - X1MIN) + 16.6
HA(1) = COND(1) / DH(1) * (.001871 * DH(1) * ABS(F(1)) / (VIS(1) * ACS(1)) + 15!) * AHT(1)
' TO MAKE EXPANSION SPACE PROCES ADIABATIC SET HA 1 =0.0
' HA(1)= 0.0
' QUENCHING IN EXP. SPACE--RADIATION .70 EFFECTIVE--HCONV=5/3600/144 MI
I = 1
XL(1) = X(1) - X1MIN + .01
XL(13) = V(13) / AP + .01
FOA = .7
IF (ABS(TMIX(1) - TMA(1))) <= (1!) THEN GOTO 42
QRAD = .173 * FOA * ((TMIX(1) / 100!) ^ 4 - (TMA(1) / 100!) ^ 4)
HRAD = QRAD / (TMIX(1) - TMA(1)) / 3600! / 144!
GOTO 43
42 HRAD = 0!
43 ' CONTINUE
PR = CP * VIS(I) / COND(I)
HCONV = .023 * REYNO(I) ^ .8 * PR ^ .4 * COND(I) / DH(I)
IF (REYNO(I) > 2100!) AND (REYNO(I) < 10000!) THEN
HCONV = HCONV * (1! + (DH(I) / XL(I)) ^ .7)
END IF
GRAETZ = REYNO(I) * PR * DH(I) / XL(I)
IF (REYNO(I) <= 2100!) AND (GRAETZ > 10!) THEN
HCONV = 1.86 * GRAETZ ^ .333 * COND(I) / DH(I)
END IF
IF (REYNO(I) <= 2100!) AND (GRAETZ <= 10!) THEN
HCONV = 5! / 3600! / 144!
END IF
IF (I) = (13) THEN GOTO 44
' ACCOUNT FOR HEAT TRANSFER TO INSULATED PORTION OF HEATER TUBES
' (ADJACENT TO EXPANSION SPACE) ALONG WITH HEAT TRANSFER TO EXPANSION
' SPACE WALL
XL1 = 1.321
REYNO1 = DH(2) * ABS(F(1)) / 8! / (VIS(1) * ACS(2)) + .0001
AHT1 = 14.955 * XL1
HCONV1 = .023 * REYNO1 ^ .8 * PR ^ .4 * COND(1) / DH(2)
HA1 = (HRAD + HCONV1) * AHT1
HA(1) = (HRAD + HCONV) * AHT(1) + HA1
' DO 40 I=2,12
FOR I = 2 TO 12
REYNO(I) = DH(I) * ABS(F(I - 1) + F(I)) / 2! / 8! / (VIS(I) * ACS(I)) + .0001
IF (I < 5) OR (I > 9) THEN
HA(I) = COND(I) / DH(I) * (.001871 * DH(I) * ABS((F(I - 1) + F(I)) / 2!) / 8! / (VIS(I) * ACS(I)) + 15!) * AHT(I) * 8!
END IF
IF (I >= 5) AND (I <= 9) THEN
HA(I) = COND(I) / DH(I) * (.06071 * DH(I) * ABS((F(I - 1) + F(I)) / 2!) / 8! / (VIS(I) * ACS(I)) + 3.7) * AHT(I) * 8!
END IF
40 ' CONTINUE
NEXT I
HA(4) = 0!
' DO 41 I=5,9
FOR I = 5 TO 9
HA(I) = RHCFAC * HA(I)
NEXT I
REYNO(13) = DH(13) * ABS(F(12)) / (VIS(13) * ACS(13)) + .0001
AHT(13) = 6.2832 * diamd * (RODYMX - RODY) + 5.5
HA(13) = COND(13) / DH(13) * (.001871 * DH(13) * ABS(F(12)) / (VIS(13) * ACS(13)) + 15!) * AHT(13)
' TO MAKE COMPRESSION SPACE PROCES ADIABATIC SET HA 13 =0.0
' HA(13)= 0.0
I = 13
IF (I = 13) THEN GOTO 42
44 HA(13) = HCONV * AHT(13)
' PAGE 54 (60 IN MANUAL PDF)-----------------
' CALCULATE RATES OF HEAT TRANSFER BETWEEN GAS AND METAL REVISE
' GAS TEMPERATURES TO ACCOUNT FOR HEAT TRANSFER TO OR FROM METAL.
' DO 120 N=1,13
FOR N = 1 TO 13
HAWC = -DTIME * HA(N) / (WG(N) * CP)
FCTR = 1! - EXP(HAWC)
DTGA = (TMAO(N) - TMIX(N)) * FCTR
Q(N) = -DTGA * WG(N) * CP / DTIME
WFCTR = FCTR * WG(N)
WDTGA = DTGA * WG(N)
SUMNUM(N) = SUMNUM(N) + WDTGA
SUMDEN(N) = SUMDEN(N) + WFCTR
120 TGA(N) = TMIX(N) + DTGA
NEXT N
'
IF (JCYCLE%) < ((NOCYC% - 5)) THEN GOTO 64
' CALCULATE DENSITY,FRICTION FACTORS, AND DELTA PRESSURES
' FOR EACH LUMP.
' DO 130 I=1,13
FOR I = 1 TO 13
DNSTY(I) = P / (R * (TGA(I) + B * P * REALGS))
NEXT I
' DO 140 I=2,12
FOR I = 2 TO 12
IF ((REYNO(I) < 1500) AND (I < 5)) OR (I > 9) THEN
FRICT(I) = 16! / REYNO(I)
END IF
IF ((REYNO(I) >= 1500!) AND (I < 5)) OR (I > 9) THEN
FRICT(I) = .046 / REYNO(I) ^ .2
END IF
IF (I >= 5) AND (I <= 9) THEN
FRICT(I) = .96 * EXP(-.019 * REYNO(I)) + .54
END IF
'140 DELP(I)=FRICT(I)*FAVG(I)*ABS(FAVG(I))*XL(I)/(32.2*DH(I)*DNSTY(I)*ACS(I)^2*6.*64.)
DELP(I) = FRICT(I) * FAVG(I) * ABS(FAVG(I)) * XL(I) / (32.2 * DH(I) * DNSTY(I) * ACS(I) ^ 2 * 6! * 64!)
NEXT I
' PRESSURE DROP CALCULATIONS
' DO 141 I=1,40
FOR I = 1 TO 40
141 DPX(I) = GAMMA
NEXT I
ACSEXP = V(1) / AD * 3.1416 * 2.75 / 8!
ACSCOM = V(13) / AP * 3.1416 * 2.75 / 8!
J% = 1
IF (F(1)) <= (0!) THEN J% = 2
Z = P
PRINT "HEATX F(1)="; F(1)
F01 = F(1)
F02 = F(2)
PRINT "HEATX F(4)="; F(4)
F04 = F(4)
F09 = F(9)
F12 = F(12)
CALL XDELP(J%, ACSEXP, ACS(2), DH(2), XL(2), DNSTY(2), VIS(2), F01, Z, COEFX(1), DPX(1), 1)
' XDELP (K%, AREAIN, AREAOT, DHX, XLGTH, ROX, VISX, FLOWIN, PXIN, XCOEF, DP, M%)
CALL XDELP(3, ACS(2), ACS(2), DH(2), XL(2), DNSTY(2), VIS(2), F01, Z, COEFX(2), DPX(2), 2)
CALL XDELP(4, ACS(2), ACS(2), DH(2), XL(2), DNSTY(3), VIS(3), F02, Z, COEFX(3), DPX(3), 3)
DPFRIC = 0!
XL2 = XL(2)
' DO 142 I=2,12
FOR I = 2 TO 12
IF (I = 2) THEN XL(2) = XL2 + 1.321
IF (I < 5) OR (I > 9) THEN KTYPE% = 10
IF (I >= 5) AND (I <= 9) THEN KTYPE% = 11
CALL XDELP(KTYPE%, ACS(I), ACS(I), DH(I), XL(I), DNSTY(I), VIS(I), FAVG(I), Z, COEFX(I + 2), DPX(I + 2), (I + 2))
' X,Z,COEFX(I+2),DPX(I+02),(I+2)) MUST BE WRONG (HARRIE) !!!!!!!!
XL(2) = XL2
DELP(I) = DPX(I + 2)
DPFRIC = DPFRIC + DELP(I)
' 142 ' CONTINUE
NEXT I
PRINT "GOTO 61 IF >0 F(4)="; F(4)
'INPUT I$
IF (F(4)) > (0!) THEN GOTO 61
PRINT " before XDELP, F(4)=0"
PRINT "ACS(4)="; ACS(4); " DH(4)="; DH(4); " XL(4)="; XL(4);
PRINT " DNSTY(4)="; DNSTY(4); " VIS(4)="; VIS(4); " F04="; F04; " Z="; Z;
PRINT " COEFX(15)="; COEFX(15); " DPX(15)="; DPX(15)
CALL XDELP(1, ACS(4), .622, DH(4), XL(4), DNSTY(4), VIS(4), F04, Z, COEFX(15), DPX(15), 15)
PRINT " AFTER XDELP, F(4)=0"
PRINT "ACS(4)="; ACS(4); " DH(4)="; DH(4); " XL(4)="; XL(4);
PRINT " DNSTY(4)="; DNSTY(4); " VIS(4)="; VIS(4); " F04="; F04; " Z="; Z;
PRINT " COEFX(15)="; COEFX(15); " DPX(15)="; DPX(15)
'INPUT I$
GOTO 62
61 'CONTINUE
CALL XDELP(1, .622, ACS(10), DH(10), XL(10), DNSTY(10), VIS(10), F09, Z, COEFX(16), DPX(16), 16)
62 ' CONTINUE
'( PAGE 55, 61 IN PDF -------------------------
CALL XDELP(5, ACS(12), .04337, DH(12), XL(12), DNSTY(12), VIS(12), F12, Z, COEFX(17), DPX(17), 17)
J% = 2
IF (F(12)) <= (0!) THEN J% = 1
CALL XDELP(J%, ACS(12), ACSCOM, DH(12), XL(12), DNSTY(12), VIS(12), F12, Z, COEFX(18), DPX(18), 18)
DPSUM = 0!
' DO 63 I=1,17
FOR I = 1 TO 17
PSUM = DPSUM + DPX(I)
63 ' CONTINUE
NEXT I
DELP(1) = DPFRIC
DELP(13) = DPSUM
DELP(2) = DELP(2) + DPX(1) + DPX(2)
DELP(3) = DELP(3) + DPX(3)
DELP(4) = DELP(4) + DPX(15)
DELP(10) = DELP(10) + DPX(16)
DELP(12) = DELP(12) + DPX(17) + DPX(18)
174 ' CONTINUE
' CALCULATE HEAT CONDUCTION (CALL HEAT CONDUCTION SUBROUTINE)
IF (NCOND) <> (3) THEN
CALL CNDCT(TMA(), QCREG, QCYL, QCONDD, QCAN, QSHTTL, MWGAS%)
END IF
NCOND = 3
64 ' CONTINUE
' REVISE REGENERATOR METAL TEMPERATURES.
' DO 177 J%=5,9
FOR J% = 5 TO 9
177 TMA(J%) = TMA(J%) + Q(J%) * DTIME / CPM(J%)
NEXT J%
' CONVERT Q(J%)'S TO F1-LBF & SUM UP FOR EACH COMPONENT.
QTOT = 0!
QEX = 0!
QHEAT = 0!
QHETN = 0!
QHETP = 0!
QREGE = 0!
QRGP = 0!
QRGN = 0!
QCOOL = 0!
QCOLP = 0!
QCOLN = 0!
QCON = 0!
' DO 179 J%=1,13
FOR J% = 1 TO 13
QADD = Q(J%) * DTIME * 777.3
IF (J% = 1) THEN QEX = QEX + QADD
IF (J% >= 2) AND (J% <= 4) THEN QHEAT = QHEAT + QADD
IF (J% >= 2) AND (J% <= 4) AND (QADD < 0!) THEN QHETN = QHETN + QADD
IF (J% >= 2) AND (J% <= 4) AND (QADD > 0!) THEN QHETP = QHETP + QADD
IF (J% >= 5) AND (J% <= 9) THEN QREGE = QREGE + QADD
IF (J% >= 5) AND (J% <= 9) AND (QADD > 0!) THEN QRGP = QRGP + QADD
IF (J% >= 5) AND (J% <= 9) AND (QADD < 0!) THEN QRGN = QRGN + QADD
IF (J% >= 10) AND (J% <= 12) THEN QCOOL = QCOOL + QADD
IF (J% >= 10) AND (J% <= 12) AND (QADD < 0!) THEN QCOLN = QCOLN + QADD
IF (J% >= 10) AND (J% <= 12) AND (QADD > 0!) THEN QCOLP = QCOLP + QADD
IF (J% = 13) THEN QCOM = QCOM + QADD
179 QTOT = QTOT + QADD
NEXT J%
' DO 185 J%=1,5
FOR J% = 1 TO 5
185 QR(J%) = QR(J%) + Q(J%) * DTIME / CPM(J%)
NEXT J%
QEXP = QEXP + QEX
IF (QEX < 0!) THEN QEXPN = QEXPN + QEX
IF (QEX > 0!) THEN qexpp = qexpp + QEX
qheatr = qheatr + QHEAT
qheatn = qheatn + QHETN
qheatp = qheatp + QHETP
qregen = qregen + QREGE
qregn = qregn + QRGN
qregp = qregp + QRGP
qcoolr = qcoolr + QCOOL
' PAGE 56 (62 IN PDF) -----------------------------
qcoolp = qcoolp + QCOLP
qcooln = qcooln + QCOLN
qcomp = qcomp + QCOM
IF (QCOM > 0!) THEN qcompp = qcompp + QCOM
IF (QCOM < 0!) THEN qcompn = qcompn + QCOM
qcndri = qcndri + QCREG * DTIME * 777.3
qcndro = qcndro + QCREG * DTIME * 777.3
qcndcl = qcndcl + QCYL * DTIME * 777.3
qcndd = qcndd + QCONDD * DTIME * 777.3
qcndcn = qcndcn + QCAN * DTIME * 777.3
QSHTL = QSHTTL * DTIME * 777.3 + QSHTL
350 ' CONTINUE
' enfrth = enfrth - AMIN1(F(4) * CP, 0!) * tg(4) * DTIME
'PRINT " ENFRTH="; ENFRTH; " F(4)="; F(4); " CP="; CP; " TG(4)="; TG(4);
'PRINT " DTIME="; DTIME
'INPUT I$
IF F(4) * CP < 0! THEN
ENFRTH = ENFRTH - (F(4) * CP) * TG(4) * DTIME
ELSE
ENFRTH = ENFRTH - (0!) * TG(4) * DTIME
END IF
' enfhtr = enfhtr + AMAX1(F(4) * CP, 0!) * tg(4) * DTIME
IF F(4) * CP > 0! THEN
ENFHTR = ENFHTR + (F(4) * CP) * TG(4) * DTIME
ELSE
ENFHTR = ENFHTR + (0!) * TG(4) * DTIME
END IF
' enfctr = enfctr - AMIN1(F(9) * CP, 0!) * tg(9) * DTIME
IF F(9) * CP < 0! THEN
enfctr = enfctr - (F(9) * CP) * TG(9) * DTIME
ELSE
enfctr = enfctr - (0!) * TG(9) * DTIME
END IF
' enfrtc = enfrtc + AMAX1(F(9) * CP, 0!) * tg(9) * DTIME
IF F(9) * CP > 0! THEN
enfrtc = enfrtc + (F(9) * CP) * TG(9) * DTIME
ELSE
enfrtc = enfrtc + (0!) * TG(9) * DTIME
END IF
' RETURN
' END
END SUB
DEFSNG I-N
SUB makestr2 (printv, length, decimals, printv$)
'call makestr2(printv,length,decimals,printv$)
DEFINT I-N
IF decimals = 0 THEN
printv$ = RIGHT$(" " + STR$(INT(printv)), length)
ELSEIF INSTR(STR$(printv), ".") THEN
IF LEN(STR$(printv)) - INSTR(STR$(printv), ".") > decimals THEN
printv = INT(printv * 10 ^ decimals) / (10 ^ decimals)
END IF
printv$ = STR$(printv) + " "
printv$ = LEFT$(printv$, INSTR(printv$, ".") + decimals)
' IF LEN(printv$) - INSTR(printv$, ".") > decimals THEN
' printv$ = printv$ + RIGHT$(" ", 1 + LEN(printv$) - decimals - INSTR(printv$, "."))
'END IF
printv$ = RIGHT$(" " + printv$, length)
ELSE
printv$ = RIGHT$(" " + STR$(printv), length - decimals - 1)
printv$ = printv$ + RIGHT$(" ", decimals + 1)
END IF
END SUB
DEFSNG I-N
SUB XDELP (K%, AREAIN, AREAOT, DHX, XLGTH, ROX, VISX, FLOWIN, PXIN, XCOEF, DP, M%)
' INDENTIFICATION OF TYPE OF LOSSES
' KTYPE%=1 CONTRACTION--KTYPE%=2 EXPANSION--KTYPE%=3 90-TURN --KTYPE4 I
' KTYPE%=5 OTHER TURN --KTYPE%=6 FRICTION-TUBE--KTYPE%=7 FRICTION-REGENE
' KTYPE%=8 FRICTION+MOMENTUM FOR TUBE--ISOTHERMAL
' KTYPE%=9 FRICTION+MOMENTUM FOR REGENRATOR--ISOTHERMAL
' PAGE 60 (66 IN PDF) ----------------------
DEFINT I-N
DIM COEFX(40), DPX(40)
' DATA GRAV/32.2, PSI/0.006945
GRAV = 32.2: PSI = .006945
GAMMA = DP
'IF (ABS(FLOWIN)) <= (0!) THEN RETURN
IF (ABS(FLOWIN)) <= (0!) THEN EXIT SUB
PR = 0!
POUT = 0!
DPMON = 0!
KTYPE% = K%
AIN = AREAIN / 144!
AOUT = AREAOT / 144!
DE = DHX / 12!
KLG = XLGTH / 12!
RO = ROX * 144! * 12!
VISC = VISX * 12
FLOW = ABS(FLOWIN) / 8!
PIN = PXIN
ROMEAN = RO
' IF(AIN-AOUT) THEN 30,30,31
IF AIN - AOUT <= 0 THEN
' 30 AMIN=AIN
AMIN = AIN
a = AIN / AOUT
' GOTO 32
ELSE
31 AMIN = AOUT
a = AOUT / AIN
END IF
32 ' CONTINUE
IF (a) < (.001) THEN a = 0!
VEL = FLOW / AMIN / RO
VELHD = RO * VEL ^ 2 / (2! * GRAV)
RE = RO * VEL * DE / VISC
' GOTO (1,2,3,4,5,6,7),KTYPE%
'ON KTYPE% GOTO 1, 2, 3, 4, 5, 6, 7,6,7,6,7
IF KTYPE% = 1 THEN
' 1 ' CONTINUE
' KTYPE%=1 CONTRACTION
' IF(RE-3000.0) THEN 12,12,11
IF (RE - 3000!) > 0 THEN ' 12,12,11
' 11 COEF=-0.4*A+0.5
COEF = -.4 * a + .5
' GOTO 100
COEF = 1! - a ^ 2 + COEF
ELSE
' 12 COEF=-0.4*A+1.0
COEF = -.4 * a + 1!
' GOTO 100
COEF = 1! - a ^ 2 + COEF
END IF
ELSEIF KTYPE% = 2 THEN
' 2 ' CONTINUE
' KTYPE%=2 EXPANSION
' IF(RE-3000.0) 21,21,22
IF (RE - 3000! <= 0) THEN
'21 COEF = (1! - 2.6 * A + 1.005 * A ^ 2) * (-1!)
COEF = (1! - 2.6 * a + 1.005 * a ^ 2) * (-1!)
' GOTO 100
COEF = 1! - a ^ 2 + COEF
ELSE
' 22 COEF=(1.0-2.0928*A+0.996*A^2)*(-1.0)
COEF = (1! - 2.0928 * a + .996 * a ^ 2) * (-1!)
' GOTO 100
COEF = 1! - a ^ 2 + COEF
END IF
' 3 ' CONTINUE
ELSEIF KTYPE% = 3 THEN
COEF = .22
' GOTO 102
' 4 ' CONTINUE
ELSEIF KTYPE% = 4 THEN
COEF = 3! * .22
' GOTO 102
' 5 ' CONTINUE
ELSEIF KTYPE% = 5 THEN
COEF = .11
' GOTO 102
' 6 ' CONTINUE
ELSEIF KTYPE% = 6 OR KTYPE% = 8 OR KTYPE% = 10 THEN
' KTYPE%=6 FRICTION-TUBE
IF (RE) < (1500!) THEN COEF = 16! / RE * (XLG / (DE / 4!))
IF (RE) >= (1500!) THEN COEF = .046 / RE ^ .2 * (XLG / (DE / 4!))
' GOTO 102
' 7 ' CONTINUE
ELSEIF KTYPE% = 7 OR KTYPE% = 9 OR KTYPE% = 11 THEN
' KTYPE%=7 FRICTION-REGENERATOR
' PAGE 61 (67 IN PDF) ---------------
COEF = .96 * EXP(-.019 * RE) + .54
COEF = COEF * (XLG / (DE / 4!))
' GOTO 102
END IF
' 100 COEF=1.0-A^2+COEF
102 DP = COEF * VELHD * PSI
IF (KTYPE%) < (8) THEN GOTO 101
' KTYPE%=8 FRICTION+MOMENTUM --TUBE KTYPE%=9 SAME FOR REGENERATOR
' KTYPE%=10 --TUBE MON.+FRIC. KTYPE%=11--FOR REGENERATOR ADI
IF (KTYPE%) >= (10) THEN GAMMA1 = GAMMA
KMA1 = VEL / SQR(GAMMA * PIN * 144! / RO + GRAV)
KMA2 = 1! / SQR(1! / XMA1 ^ 2 - GAMMA * COEF)
FMA21 = 1!
IF (ABS((XMA2 - XMA1) / XMA1)) <= (.001) THEN GOTO 52
51 XMA0 = (XMA1 + XMA2) * .5
FMA21 = ((GAMMA - 1!) * XMA0 ^ 2 + 2!) / ((GAMMA - 1!) * XMA1 ^ 2 + 2!)
IF (KTYPE%) < (10) THEN FMA21 = 1!
IF (KTYPE%) < (10) THEN GAMMA1 = 1!
XMA2 = 1! / SQR(1! / XMA1 ^ 2 - GAMMA * COEF + GAMMA * (GAMMA1 + 1!) / 2! / GAMMA1 * LOG((XMA1 / XMA0) ^ 2 * FMA21))
IF (ABS(XMA2 - XMA0)) > (.001) THEN GOTO 51
52 TR = 1! / FMA21
VR = XMA2 / XMA1 * SQR(1! / FMA21)
PR = XMA1 / XMA2 * SQR(1! / FMA21)
POUT = PR * PIN
DP = PIN - POUT
101 ' CONTINUE
IF (FLOWIN) <= (0!) THEN DP = (-1!) * DP
COEFX(M%) = COEF
DPX(M%) = DP
' RETURN
' END
END SUB
'==================================================================