Stirling engine, heatpump and refrigerator quickbasic program

The engine of the future project

The following is a translation of the ADIABATIC SECOND ORDER DESIGN PROGRAM written in Fortran by P.A. Rios.

I found it on:

http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19830022057_1983022057.pdf

and translated it into Quick Basic.

In the future I will make 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 quickbasic or qbasic and load the program. The input is in the DATA lines on the end, and can be changed there.

Now I'm working on a better version.

This page is ment for progammers.

For the html-javascript version of the Isothermal second order calculation Program (Written by Wiliam R. Martini in Fortran) Click here

The program is in the source code.

A more advanced version of this program will be published in a Quick Basic version on this site soon.

Click here for HOME

'==================================================================== DECLARE FUNCTION xlarge! (X!(), NDIV%) DECLARE FUNCTION small! (X!(), NDIV%) DECLARE SUB makestr (printv!, length, decimals, printv$) DECLARE SUB grafiek (opdracht$, x1!, y1!, x2!, y2!, kleur!) DECLARE SUB closefilenr (filecount%, FILENUMMER%) DECLARE SUB givefilenr (filecount%, FILENUMMER%) DECLARE SUB pdint (X!, DMW!(), DMC!(), RVT!, DVC!(), NDIV%, DM!, PR!(), XINT!, DPR!(), XI1!, XI2!, XNHT!) DECLARE SUB VOLW (W!(), WI!(), DW!(), DWI!(), CFI!, SFI!, ZZW!, NDIV%, SIFI!(), COFI!(), SALF!(), CALF!(), DALF!) DECLARE SUB VOLC (DALF!, NF%, C!(), CI!(), DC!(), DCI!(), ZZC!, NDIV%, SIFI!(), COFI!(), SALF!(), CALF!(), PDR!, ALF!) DECLARE FUNCTION plot! (X(), H!) 'www.scribd.com/doc/20191898/Stirling-Engine-Design-Manual ' Rios 41.23 job (),kosdruk-notify=wrmjc*region=120k 'exec fortgclg 'fort.sysin dd * 'Expects number of ducts= number of regerators! MS% = 1' mode switch for pumping loss Martini (off in Richard's exemple) CLS 'AFC = Cooler free flow area, cm^2 'AFH = Heater free flow area, cm^2 'AFR = Regenerator free flow area, cm^2 'ALF = 4.7123889 (270 degrees) 'ARG = sin(PV angle) 'BDR = (input) Regenerator diameter, cm 'BEC = (input) Piston end clearence, cm 'BPD = (input) Piston diameter, cm 'BPL = (input) Hot cap length, cm 'BRC = (input) Piston gap, cm (of brg) 'BRL = (input) Regenerator length, cm 'BRO = (input) Regenerator density factor 'BST = (input) Piston stroke, cm 'BTC = (input) Effective cold temperature, K 'BTC1 = Cold metal temperature, K 'BTR = Regenerator temperature, K 'BTW = (input) Hot effective gas temperature, K 'BTW1 = Hot metal temperature, K 'BWD = (input) Effective regenerator wire diameter, cm 'C() = Cold space as fraction of the stroke amplitude at mid-increment ' C() varies from 0 to 2 and back. 'CALF() = Sin change per radian increment 'CALFP = CALF() 'CFI = Cos of phase angle 'CI() = Same as C() for beginning of increment 'CMMAX = Largest cold dimensionless mass 'CMU = Cold hydrogen viscosity 'CNTU = Number of heat transfer units in cold space 'COFI() = Cos values for cold space 'CON = Conduction loss, watts 'COR = PMAX ^ (XNHT - 2) * DALF ^ (XNHT - 1) REGEL 340, 345, 346 'CP1 = Hydrogen heat capacity 'CRC = SQR(ZZC^2-CALF()^2) 'CRW = CRC in hot space 'CTD = (input) Cooler tube inside diameter, cm 'CTLL = (input) Total cooler tube length, cm 'CTLS = (input) Cooler cool tube length, cm 'CV1 = Hydrogen heat capacity 'DALF = 2phi/ndiv 'DC() = Angle derivative of C() 'DCI() = Angle derivative of CI() 'DDD = (input) Cooler duct diameter, cm 'DLL = (input) Cooler duct length, cm 'DM = Sum of changes in mass (DMRE) 'DMC = Cold dimension mass change XDMC() 'DMRE = Sum of changes in mass (DM) 'DMW = Hot dimensionless mass change XDMW() 'DMX = Dimensionless change in mass relating to X, the fraction from ' the cold end 'DP = Change in pressure 'DPR() = DP array 'DTC = Cooler metal temperature - effective temperature 'DTW = HOT metal temperature - effective temperature 'DTH = Delta Th 'DV = Dead volume, cm^3 'DVC = DC() 'DVCI = DCI() 'DVW = DW() 'DVWI = DWI() 'DW()= Angle derivative of W() 'DWI() = Angle derivative of WI() 'DX = 1/xnds 'EX1 = 1-XNHT 'EX2 = 2-XNHT 'FC = Cold friction factor 'FFF = Friction flow credit, watts 'FH = Hot friction factor 'FR() = Regenerator friction factors (3 pts.) 'FI = Phase angle, rad. 'FI1 = Phase angle in deg. 'FIPV = pv angle (output) arcsin (ARG) 'FR() = Regenerator friction factor 'G1 = Y value subplot 'G2 = Y value subplot 'GDMS() = Calculated mass flow values 'GGV = Dead volume at side of hot cap, cm^3 'GINT() = Flow loss variable 'GI2() = Pressure drop value 'GI3() = Pressure drop value 'GLH = Heater pressure drop integral 'GLR = Regenerator pressure drop integral 'GLS = Cooler pressure drop integral 'H(1) = Fraction of total reduced dead volume from cold end to midway ' in cooler. 'H(2) = Fraction of total reduced dead volume from the cold end through ' the cooler 'H(3) = Fraction of total reduced dead volume from the cold end through ' half the regenerator 'H(4) = Fraction of total reduced dead volume through de regenerator 'H(5) = Fraction of total reduced dead volume through the middle of the gas ' heater (1-H(5) includes the rest of the heater and clearence on the ' end and sides of the hot cap) 'HAC = Cold active volume amplitude, cm^3 'HAV = Hot active volume amplitude, cm^3 'HCV = Reduced cooler and cold ducting dead volume, dimensionless 'HEC = Reduced cold end clearence dead volume, dimensionless 'HGV = reduced hot cap gap dead volume, dimensionless 'HHC = reduced hot clearance dead volume, dimensionless 'HHV = reduced heater dead volume, dimensionless 'HMU= Hot hydrogen viscosity 'HNTU = .1044 * HTLS / (HTD * REH ^ .2)' regel 419,420 hydrogen ?? 'HRV = Reduced regenerator dead volume, dimensionless 'HT = Basic heat input, watts 'HTD = (input) Heater tube inside diameter, cm 'HTE = Heat to engine, watts 'HTLL = (input) Total heater tube length, cm 'HTLS = (input) Heated heater tube length, cm 'HTW = Hot end heat transfer integral, dimensionless 'IND() = Array that shows if mass change is positive, or negative in ' warm and cold sides 'J% = V [Temporary angle variable, radians 'K% = 1 if warm mass change is positive, 2 if negative 'L% = 1 if cold mass change is positive, 2 if negative 'LUP% = Iterational counter 'M% = X value for plot calculation 'MBR% = (input) Number of regenerators 'MCT% = (input) Number of cooler tubes per cylinder 'MHT% = (input) Number of heater tubes for cylinder 'MW% = Dimensionless mass in hot space = (mass, grams)(R)(BTW)/(PMX1(HAV)) 'N% = NDIV% or X value for plot subroutine 'NN% = 1 up to phase angle, 2 after 'NDIV% = (input) Number of divisions per crank rotation (must be a multiple ' of 4 so that the phase angle at 90 degrees can be an even number of ' divisions) ' (Programm must be revised if NDIV% is not 360) 'NDIV1% = NDIV% + 1 'NDS% = (input) Number of divisions in dead space 'NE% = ndiv/4 + 1 'NET% = (input) Regenerator filler option (+5 = metnet, -5 = screen) 'NF% = ndiv/4 'NFF% = NF%+1 'NFI% = (input) (phase angle)(NDIV%)/360 'NFIN% = Main loop final counter, for first part = phase angle, for second ' part = end of cycle 'NIN% = NDS%+1 'NITE% = Cycle counter (counts to 15) 'NL% = (ndiv/2)+1 'NLOP% = Option counter limits changes in option to 7 (removed in final ' version) '0; = IND(K%,L%)-1,2,3 or 4, starts as 1 'NOC% = (input) number of cylinders 'NS% = (ndiv/4)+2 'NST% = Main loop initial counter, for the first part =1, for the second ' part = phase angle 'NT% = (ndiv/4)+2 'NWR% = (input) Governs printout, zero for overall results only, different ' from zero added pv data 'P = pressure, dimensionless 'PALF = Thermal diffusitivity of piston 'PDR = (input) Piston rod diameter, cm 'PI4 = phi/4 = .78539816 'PIE = 3.1415926# 'PAVG = Dimensionless average pressure 'PMAX = Maximum pressure, dimensionless 'PMIN = Minimum pressure, dimensionless 'PMX = Maximum pressure (MPa) 'PMX1 = (input) Avg. pressure MPa 'PR() = Pressure, dimensionless, fraction of maximum pressure 'PO = Basic power, watts 'POF = Net power, watts 'PS = Dimensiomless pressure from end of previous cycle 'PW = Pressure at halfway point for increment 'QB = Beta for shuttle heat loss calculation 'QCP = Cooler windage, watts 'QDK = Reheat factor 'QFS = Pumping loss factor 'QHC = Shuttle loss, watts 'QHG = Pumping loss, watts 'QHP = Heater windage, watts 'QHR = Reheat loss, watts 'QLM = Reheat factor 'QL1 = Shuttle factor 'QNPH = Reheat pressurization effect 'QNTU = Regenerator transfer units, dimensionless 'QP = Windage factor 'QR() = Regenerator windage loss values, watts 'QR1, QR2, QR3 line 409,410,411,412 'QRP = Average regenerator windage, watts 'R = Gas constant, joules/(gm)(K%) 'R2 = Constant = R(gc)2 'RE() = Regenerator Reynolds number in cold, middle, and hot part 'REC = Cold Reynolds number 'REH = Hot Reynolds number 'RER = Regenerator Reynolds number 'RMU = Regenerator hydrogen viscosity, g/cm sec. 'RNTU = Regenerator heat transfer units 'RP = Maximum pressure/minimum pressure 'RVT = Displacer mass ratio 'S = Pressure at halfway point, dimensionless 'SALF() = Sin values for cold space 'SALFP = Average sin values for cold space 'SFI = Sin of phase angle 'SHR = (input) Specific heat ratio for working gas 'SIFI() = SALF() 'SIFIP = SALF(1) 'SMC = Cold mass + 1/2 change in mass 'SMW = Hot mass + 1/2 change in mass 'SPD = (input) Engine speed, rad/sec 'TEC() = Dimensionless cold gas temperature 'TEST = Ensures that difference in dimensionless mass <.001 'TEST1 = Ensures that difference in dimensionless pressure <.005 'TEW() = Dimensionless hot gas temperature 'TMPC = Average TEC() 'TMPW = Average TEW() 'TCDM = Dimensionless average cold temperature for entire cycle 'TWDM = Dimensionless average warm temperature for entire cycle 'UDM() = Critical mass flow values from subplot 'UIN() = Critical pressure drop integral values from subplot 'UI23, 24, 33, 34 = Critical pressure drop values 'UPA = Power piston area, cm^2 'UTR = Temperature ratio = (hot metal temp, K%)/(cold metal temp, K%) 'VC = C() 'VCC = Cold volume cm^3 'VCD = Cold dead volume, cm^3 'VCI = CI() 'VD = Reduced dead volume, dimensionless 'VH = Hot volume, cm^3 'VHD = Hot dead volume, cm^3 'VRD = Regenerator dead volume, cm^3 'VT = Total volume, cm^3 'VW = W() 'VWI = WI() 'W() = Hot space as fraction of the stroke amplitude, calculated at mid ' increment 'WC = Dimensionless cold work 'WI() = Same as W() for beginning of increment 'WMMAX = Largest hot dimensionless mass 'WW = Hot work, dimensionless 'X = short term variable 'XDMC() = Change in cold mass, grams 'XDMW() = Change in hot mass, grams 'XI1 = Pressure drop integral - account for the relationship between the ' shapes of mass and pressure fluctuations 'XI2 = Influence of mass flow time variation on the heat transfer 'XI3 = xi1/xi2 'XINT = Basic pressure drop integral - for windage 'XMC = Cold gas mass, relative to total inventory 'XMCX() = Cold gas mass, grams 'XMT() = Total mass, grams 'XMW = Hot gas mass, relative to total inventory 'XMWS = Hot dimensionless gas mass from previous cycle 'XMWX() = Hot gas mass, grams 'XND = NDIV% 'XNDS = NDS% 'XNHT = (input) Value for exponent in heat transfer relation of regenerator ' matrix 'XX = Short term variable 'Y =|DMX| 'ZEF = Indicated efficiency, % 'ZZC = (input) Connecting rod length/.5 stroke for cold piston 'ZZW = (input) Connecting rod length/.5 stroke for hot piston '########### I N I T ############### 'D.3 Fortran Listing 'If memory is a problem all dimensions to 720 can be changed to 361 DIM XDMC(720), XDMW(720), IND(2, 2), PR(720), XMCX(720), XMWX(720) DIM TEC(720), TEW(720), XMT(720) DIM W(720), C(720), WI(720), CI(720), DW(720), DC(720), DWI(720) DIM DCI(720), DPR(720) DIM SIFI(720), COFI(720), SALF(720), CALF(720), H(5), UIN(5) DIM UDM(5), GDMS(20), GINT(20), GI2(20), GI3(20), RE(3), FR(3), QR(3) 'constants are now setup PIE = 4 * ATN(1)'PIE = 3.1415926# PI4 = ATN(1)'PI4 = .78539816# GOSUB 9900 'read variable text strings '##### MAIN ########### 'output basic file CALL givefilenr(filecount%, OUTFILE%) bestand$ = "STIRLOUT.TXT" OPEN bestand$ FOR OUTPUT AS #OUTFILE% PRINT #OUTFILE%, "Stirling program output: "; bestand$ GOSUB inputs 'load defaulds IF (NWR%) <> 0 THEN GOSUB printout CALL closefilenr(filecount%, OUTFILE%) INPUT "End of programm"; I$ END inputs: 'Specified engine dimensions and operating conditions are read 'in from data cards RESTORE 9999 READ name$ IF name$ = "end" OR name$ = "END" OR name$ = "" THEN RETURN READ ZZC 'stops if next data is end, END or blank IF (ZZC) = 0 THEN END ' end of program END IF 'read(5.5) ZZC,ZZW,XNHT,PMX1,SPD,NET%,NOC% 'read (5+1) SHR,NFI%,NDIV%,NWR%,NDS%,DDD,DLL 'read (5+1710) BTC,BTW,BPD,BST,BPL,BRC,BEC 'read (5+1720) PDR,BRO,BDR,BWD,BRL,MBR%,MCT% 'read (5+1730) CTD,CTLL,CTLS,HTD,HTLL,HTLS,MHT% READ ZZW, XNHT, PMX1, SPD, NET%, NOC% READ SHR, NFI%, NDIV%, NWR%, NDS%, DDD, DLL READ BTC, BTW, BPD, BST, BPL, BRC, BEC READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT% READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT% READ CON PRINT ZZC, ZZC$ PRINT ZZW, ZZW$ PRINT XNHT, XNHT$ PRINT PMX1, PMX1$ PRINT SPD, SPD$ PRINT NET%, NET$ PRINT NOC%, NOC$ PRINT PRINT SHR, SHR$ PRINT NFI%, NFI$ PRINT NDIV%, NDIV$ PRINT NWR%, NWR$ PRINT NDS%, NDS$ PRINT DDD, DDD$ PRINT DLL, DLL$ INPUT I$ PRINT BTC, BTC$ PRINT BTW, BTW$ PRINT BPD, BPD$ PRINT BST, BST$ PRINT BPL, BPL$ PRINT BRC, BRC$ PRINT BEC, BEC$ PRINT 'READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT% PRINT PDR, PDR$ PRINT BRO, BRO$ PRINT BDR, BDR$ PRINT BWD, BWD$ PRINT BRL, BRL$ PRINT MBR%, MBR$ PRINT MCT%, MCT$ INPUT I$ 'READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT% PRINT CTD, CTD$ PRINT CTLL, CTLL$ PRINT CTLS, CTLS$ PRINT HTD, HTD$ PRINT HTLL, HTLL$ PRINT HTLS, HTLS$ PRINT MHT%, MHT$ PRINT CON, CON$ INPUT I$ PRINT #OUTFILE%, "name= "; name$ PRINT #OUTFILE%, ZZC, ZZC$ PRINT #OUTFILE%, ZZW, ZZW$ PRINT #OUTFILE%, XNHT, XNHT$ PRINT #OUTFILE%, PMX1, PMX1$ PRINT #OUTFILE%, SPD, SPD$ PRINT #OUTFILE%, NET%, NET$ PRINT #OUTFILE%, NOC%, NOC$ PRINT #OUTFILE%, PRINT #OUTFILE%, SHR, SHR$ PRINT #OUTFILE%, NFI%, NFI$ PRINT #OUTFILE%, NDIV%, NDIV$ PRINT #OUTFILE%, NWR%, NWR$ PRINT #OUTFILE%, NDS%, NDS$ PRINT #OUTFILE%, DDD, DDD$ PRINT #OUTFILE%, DLL, DLL$ PRINT #OUTFILE%, PRINT #OUTFILE%, BTC, BTC$ PRINT #OUTFILE%, BTW, BTW$ PRINT #OUTFILE%, BPD, BPD$ PRINT #OUTFILE%, BST, BST$ PRINT #OUTFILE%, BPL, BPL$ PRINT #OUTFILE%, BRC, BRC$ PRINT #OUTFILE%, BEC, BEC$ PRINT #OUTFILE%, 'READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT% PRINT #OUTFILE%, PDR, PDR$ PRINT #OUTFILE%, BRO, BRO$ PRINT #OUTFILE%, BDR, BDR$ PRINT #OUTFILE%, BWD, BWD$ PRINT #OUTFILE%, BRL, BRL$ PRINT #OUTFILE%, MBR%, MBR$ PRINT #OUTFILE%, MCT%, MCT$ PRINT #OUTFILE%, 'READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT% PRINT #OUTFILE%, CTD, CTD$ PRINT #OUTFILE%, CTLL, CTLL$ PRINT #OUTFILE%, CTLS, CTLS$ PRINT #OUTFILE%, HTD, HTD$ PRINT #OUTFILE%, HTLL, HTLL$ PRINT #OUTFILE%, HTLS, HTLS$ PRINT #OUTFILE%, MHT%, MHT$ PRINT #OUTFILE%, 'Hot and cold metal temperature are set at given hot and cold 'gas temperature: BTW1 = BTW BTC1 = BTC 'Power piston area UPA = BPD * BPD * PI4 PRINT "UPA = power piston area ="; UPA; " CM^2" 'temperarure ratio, inverse of normal UTR = BTW / BTC 'hot cylinder volume amplitude is reference volume for all 'dimensionless volumes HAV = UPA * BST / 2 PRINT "HAV = hot active volume amplitude ="; HAV; " CM^3" GGV = PIE * BPD * BPL * BRC PRINT "GGV = dead volume at side of hot cap ="; GGV; " CM^3" 'initialize loop counter set for 3 iterations 'LUP% = 1 '339 **** Return point after 1 cycle **** FOR lup = 1 TO 3 'THE FOLLOWING ARE REDUCED DIMENSIONLESS DEAD VOLUMES RELATING TO 'MAXIMUM AND MINIMUM EFFECTIVE GAS TEMP, AND HAV HGV = GGV * 2 * BTW / (HAV * (BTW + BTC)) HEC = (UPA - PDR * PDR * PI4) * BEC * UTR / HAV AFC = PI4 * CTD * CTD * MCT%'cooler_free_flow_area HCV = (AFC * CTLL + DDD * DDD * PI4 * DLL * MBR%) * UTR / HAV'cooler+cold_ducting vol. AFR = PI4 * BDR * BDR * MBR% * BRO 'regenerator flow aeria HRV = AFR * BRL * 2 * BTW / (HAV * (BTW + BTC))' reduced regenerator dead volume AFH = PI4 * HTD * HTD * MHT% 'heater free flow earia HHV = AFH * HTLL / HAV 'reduced heater dead volume HHC = UPA * BEC / HAV 'reduced hot clearence volume VD = HGV + HEC + HCV + HRV + HHV + HHC '52 The following are hot, regenerator and cold dead volumes in cm^2 VHD = AFH * HTLL + UPA * BEC 'HOT_DEAD_VOLUME VRD = GGV + AFR * BRL 'REGENERATOR_DEAD_VOLUME VRD VCD = (UPA - PDR * PDR * PI4) * BEC + HCV * HAV / UTR 'COLD_DEAD_VOLUME IF lup = 1 THEN PRINT "HGV = reduced hot gap dead volume ="; HGV PRINT "HEC = cold end clearence dead volume dimensionless ="; HEC PRINT "AFC = cooler free flow area ="; AFC; " CM^2" PRINT "HCV = cooler + cold ducting volume ="; HCV PRINT "AFR = regenerator flow aeria ="; AFR; " CM^2" PRINT "HRV = reduced regenerator dead volume ="; HRV PRINT "AFH = heater free flow earia ="; AFH; " CM^2" PRINT "HHV = reduced heater dead volume ="; HHV PRINT "HHC = reduced_hot_clearance_dead_volume ="; HHC PRINT "VD = reduced dead volume="; VD INPUT I$ HAC = HAV - PI4 * PDR * PDR * BST / 2 CALL makestr(UPA * BEC, 12, 4, printv$): PRINT printv$; PRINT " = piston end clearance dead volume, cm^3" CALL makestr(AFH * HTLL, 12, 4, printv$): PRINT printv$; PRINT " = heater tube dead volume, cm^3" PRINT "+ -----------" CALL makestr(VHD, 12, 4, printv$): PRINT printv$; PRINT " = hot dead volume, cm^3" CALL makestr(2 * HAV, 12, 4, printv$): PRINT printv$; PRINT " = hot active volume, cm^3" PRINT "+ -----------" CALL makestr(VHD + 2 * HAV, 12, 4, printv$): PRINT printv$; PRINT " = total max hot volume, cm^3" PRINT CALL makestr((UPA - PDR * PDR * PI4) * BEC, 12, 4, printv$): PRINT printv$; PRINT " = piston end clearance dead volume, cm^3" CALL makestr(AFC * CTLL, 12, 4, printv$): PRINT printv$; PRINT " = cooler tube dead volume, cm^3" CALL makestr(HCV * HAV / UTR - AFC * CTLL, 12, 4, printv$): PRINT printv$; PRINT " = connecting ducts dead volume, cm^3" PRINT "+ -----------" CALL makestr(VCD, 12, 4, printv$): PRINT printv$; PRINT " = cold dead volume, cm^3" CALL makestr(2 * HAC, 12, 4, printv$): PRINT printv$; PRINT " = cold active volume, cm^3" PRINT "+ -----------" CALL makestr(VCD + 2 * HAC, 12, 4, printv$): PRINT printv$; PRINT " = total max cold volume, cm^3" PRINT CALL makestr(AFR * BRL, 12, 4, printv$): PRINT printv$; PRINT " = regenerator dead volume, cm^3" CALL makestr(GGV, 12, 4, printv$): PRINT printv$; PRINT " = hot gap dead volume, cm^3" PRINT "+ -----------" CALL makestr(VRD, 12, 4, printv$): PRINT printv$; PRINT " = regenerator dead volume, cm^3" PRINT CALL makestr(VRD + VCD + VHD, 12, 4, printv$): PRINT printv$; PRINT " = total dead volume, cm^3", INPUT I$ END IF 'on the first iteration only the input values are written out '58 IF(LUP%-1)343,343,346 IF lup = 1 THEN '343 HAC = HAV - PI4 * PDR * PDR * BST / 2 PRINT 'write(6,1805) PRINT "Cold metal temp (k) ="; BTC, , PRINT "Hot metal temp (k) ="; BTW PRINT "Piston diameter(cm) ="; BPD, PRINT "Piston stroke (cm) ="; BST PRINT "Hot cap length (cm) ="; BPL, , PRINT "Piston gap (cm) ="; BRC PRINT "Piston end clearence(cm)="; BEC, PRINT "Piston rod diameter (cm)="; PDR PRINT "Regenerator porosity ="; BRO, PRINT "Regenerator diameter(cm)="; BDR PRINT "Regenerator wire dia(cm)="; BWD, PRINT "Regenerator length (cm)="; BRL PRINT "Number of regenerators ="; MBR%, PRINT "Number of cooler tubes ="; MCT% PRINT "Cool tube diameter (cm) ="; CTD, PRINT "Total cool tube length(cm)="; CTLL ' INPUT I$ PRINT "Cool ct length (cm)="; CTLS, PRINT "Heater tube diameter (cm)="; HTD PRINT "Total heater tube length(cm)="; HTLL, PRINT "Heat heater tube length (cm)="; HTLS PRINT "Number of heater tubes ="; MHT%, PRINT "Heat transfer exponent ="; XNHT PRINT "Average pressure MPA ="; PMX1, PRINT "Engine speed (rad/sec)="; SPD PRINT "NETNET option ="; NET%, , PRINT "Number of cylinders ="; NOC% PRINT "Cold crank ratio ="; ZZC, , PRINT "Hot crank ratio ="; ZZW GOSUB printinput 'INPUT I$ END IF '68 The reduced dead volume is divided into fractions, re-evaluated 'each iteration because of change in tempertatures '346 H(1) = (HEC + HCV / 2) / VD H(2) = H(1) + HCV / (2 * VD) H(3) = H(2) + HRV / (2 * VD) H(4) = H(2) + HRV / VD H(5) = H(4) + HHV / (2 * VD) '75 Calculations from 349 to 295 are on first iteration only. Engine 'volumes and volume derivatives are calculated and decision matrix 'is defined. 'IF (LUP% - 1) 349,349,295 IF lup = 1 THEN '349 XND = NDIV% NDIV1% = NDIV% + 1 'DALF=0.0174 radians/incriment DALF = 2 * PIE / XND' 6.2831853# / XND NT% = NDIV% / 4 + 2 NE% = NT% - 1 'NOT USED ????????????????? '85 Calculation starts at 270 deg=4.71 radians in radians if the 'front side of the cold piston is used and at 90 deg=1.57 radians 'in radians if the back side is used and room must be allowed for 'a piston drive rod. This gives the proper curve shape. 'Calculation always starts with zero cold live volume. 'IF (PDR) 4060,4060,4070 IF PDR > 0 THEN ALF = PIE / 2' 1.5707963# ELSE ALF = PIE + PI4' 4.7123889# END IF NF% = NDIV% / 4 '95 Call subroutine to calculate duplacable space above or under cold 'piston at the midpoint and at the beginning of each angle 'increment as a fraction of the piston stroke amplitude. 'Subroutine also calculates derivatives to be used later. CALL VOLC(DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(), COFI(), SALF(), CALF(), PDR, ALF) 'subroutine VOLC(DALF,NF%,C(),CI(),DC(),DCI(),ZZC,NDIV%,SIFI(),COFI(),SALF(),CALF(),PDR,ALF) '101 The phase angle must be 90 deg. for a heat engine and 270 deg for 'a heat pump FI = DALF * NFI% FI1 = FI * 180 / PIE 'WRITE(6,11) SHR, NDIV%, FI1, NDS%, DDD, DLL 'CALL makestr(SHR, printv$) PRINT "Specific heat ratio ="; SHR, ', printv$, 'CALL makestr(NDIV%, printv$) PRINT , "Div. per cycle ="; NDIV%' , printv$ ' CALL makestr(FI1, printv$) PRINT "Phase angle (deg.) ="; FI1, , ', printv$, ' CALL makestr(NDS%, printv$) PRINT "Incr. in dead space int.="; NDS%', printv$ ' CALL makestr(DDD, printv$) PRINT "Duct diameter (cm) ="; DDD, ', printv$, 'CALL makestr(DLL, printv$) PRINT , "Duct length (cm)="; DLL', printv$ INPUT I$ SFI = SIN(FI) CFI = COS(FI) 'Now the hot space fractions axd, derivatives are calculated for 'each increment CALL VOLW(W(), WI(), DW(), DWI(), CFI, SFI, ZZW, NDIV%, SIFI(), COFI(), SALF(), CALF(), DALF) 'choice matrix is defined- see note #2 IND(1, 1) = 1 IND(1, 2) = 3 IND(2, 1) = 4 IND(2, 2) = 2 END IF 'First guess at initial quantities 'Counter NN% is 1 at start OF cycle when Give cold volume is dsad 'It becomes 2 after live hot volume becomes zero, phase angle is '270 deg. for a heat pump and 90 deg. for a heat engine. '295 NN% = 1 'None of the gas mass is assumed to be in the cold space to start, 'volume is zero XMC = 0 'Relative pressure is assumed to be maximum, approximation. P = 1 'At start all gas mass is assumed to be in the hot space, ignores 'dead volume XMW = 1 - CFI 'mass from pevious cycle is assumed at start to be zero, to assure 'at least 2 cycles to convergence. XMWS = 0 'previous cycle pressure PS = 1 '135 initial assumption for value of IND(k,l) is correct for heat 'engine wrong for a heat pump. NO = 1 'initialize dimensionless works WW = 0 WC = 0 NITE% = 1 NST% = 1 NFIN% = NFI% 'displaced mass ratio RVT = HAC * UTR / HAV CI(NDIV1%) = CI(1) WI(NDIV1%) = WI(1) '148 ****************************************** 'start of main do loop, return point after each increment RELOOP = 1 WHILE RELOOP = 1 RELOOP = 0 '404 'do 102 I=NST%,NFIN% FOR I = NST% TO NFIN% 'transfer volumes and derivtives from storage VW = W(I) VC = C(I) VWI = WI(I) VCI = CI(I) DVW = DW(I) DVC = DC(I) DVWI = DWI(I) DVCI = DCI(I) 'splits to 4 options 'go to (201,202,203,204),NO 'IF NO = 1 THEN gOTO 201 'IF NO = 2 THEN gOTO 202 'IF NO = 3 THEN gOTO 203 'IF NO = 4 THEN gOTO 204 IF NO = 1 THEN '201 integration program for mass increasing in both hot and cold spaces. 'NO=1 (see note #3 'Compute pressure change based upon initial conditions DP = -SHR * P * (RVT * DVCI + DVWI) / (RVT * VCI + VWI + SHR * VD) * DALF 'change in pressure dispMr coldV hotV heatR rdeadV 'finds pressure at mid increment S = P + DP / 2 'calculates final pressure change basic upon mid point values DP = -SHR * S * (RVT * DVC + DVW) / (RVT * VC + VW + SHR * VD) * DALF 'calculates mass changes DMW = S * DVW * DALF + VW * DP / SHR DMC = -(DMW + VD * DP) / RVT 'determins choice matrix IF DMW >= 0 THEN K% = 1 ELSE K% = 2 END IF IF DMC >= 0 THEN L% = 1 ELSE L% = 2 END IF NO = IND(K%, L%) 'if choice is changed next iteration will be through a different 'option 'gOTO 400 ELSEIF NO = 2 THEN '202 integration program for mass decreasing in both hot and cold spaces. 'NO=2 (see option 1 for detailed explanation) IF XMC < 0 THEN XMC = 0! END IF IF XMW < 0 THEN XMW = 0 END IF DP = -SHR * (XMC * RVT * DVCI / VCI + XMW * DVWI / VWI) / (XMC * RVT / P + XMW / P + SHR * VD) * DALF DMC = XMC * (DVCI * DALF / VCI + DP / SHR / P) DMW = -RVT * DMC - VD * DP S = P + DP / 2 SMC = XMC + DMC / 2 SMW = XMW + DMW / 2 DP = -SHR * (SMC * RVT * DVC / VC + SMW * DVW / VW) / (SMC * RVT / S + SMW / S + SHR * VD) * DALF DMC = SMC * (DVC * DALF / VC + DP / SHR / S) DMW = -RVT * DMC - VD * DP IF DMW <= 0 THEN K% = 2 ELSE K% = 1 END IF IF DMC <= 0 THEN L% = 2 ELSE L% = 1 END IF NO = IND(K%, L%) 'gOTO 400 ELSEIF NO = 3 THEN '203 integration program for mass decreasing in cold space and increasing 'in hot space, NO=3 (see option 1 for detailed explanation) ' if(XMC) 704,703,703 IF XMC < 0 THEN XMC = 0 END IF DP = -SHR * (P * DVWI + XMC * RVT * DVCI / VCI) / (VWI + XMC * RVT / P + SHR * VD) * DALF DMC = XMC * (DVCI * DALF / VCI + DP / SHR / P) DMW = -RVT * DMC - VD * DP S = P + DP / 2 SMC = XMC + DMC / 2 SMW = XMW + DMW / 2 DP = -SHR * (S * DVW + SMC * RVT * DVC / VC) / (VW + SMC * RVT / S + SHR * VD) * DALF DMC = SMC * (DVC * DALF / VC + DP / SHR / S) DMW = -RVT * DMC - VD * DP IF DMW >= 0 THEN K% = 1 ELSE K% = 2 END IF IF DMC <= 0 THEN L% = 2 ELSE L% = 1 END IF NO = IND(K%, L%) 'gOTO 400 ELSEIF NO = 4 THEN '204 Integration program for mass decreasing in cold space and decreasing 'in hot space, NO=4 (see option 1 for detailed explanation) IF XMW < 0 THEN XMW = 0 END IF DP = -SHR * (P * RVT * DVCI + XMW * DVWI / VWI) / (RVT * VCI + XMW / P + SHR * VD) * DALF DMW = XMW * (DVWI * DALF / VWI + DP / SHR / P) DMC = -(DMW + VD * DP) / RVT S = P + DP / 2 SMC = XMC + DMC / 2 SMW = XMW + DMW / 2 DP = -SHR * (S * RVT * DVC + SMW * DVW / VW) / (RVT * VC + SMW / S + SHR * VD) * DALF DMW = SMW * (DVW * DALF / VW + DP / SHR / S) DMC = -(DMW + VD * DP) / RVT IF DMW <= 0 THEN K% = 2 ELSE K% = 1 END IF IF DMC >= 0 THEN L% = 1 ELSE L% = 2 END IF NO = IND(K%, L%) 'gOTO 400 ELSE INPUT "NO IS NOT 1,2,3 OR 4"; I$ END IF '400 increment pressure and mass P = P + DP XMC = XMC + DMC XMW = XMW + DMW 'calculates works PW = P - DP / 2 WC = WC + PW * DVC * DALF WW = WW + PW * DVW * DALF 'records results into arrays PR(I) = P DPR(I) = DP XMCX(I) = XMC XMWX(I) = XMW XDMC(I) = DMC XDMW(I) = DMW ' ******** end of main loop ********** NEXT I IF NN% = 1 THEN 'reset main loop for last part of cycle NST% = NFI% + 1 NFIN% = NDIV% NN% = 2 RELOOP = 1 ' gOTO 404 ELSE '205 Tests for convergence at end of cycle. The change in the fraction 'of mass in the hot space from one cycle to the next must be less 'than 0.1%, and the change in the pressure from one cycle to the next 'must be less than 0.5%. However, NO more than 15 cycles are 'allowed TEST = SQR((XMWS - XMW) ^ 2) TEST1 = SQR((PS - P) ^ 2) IF (NITE% <= 15) THEN IF (TEST > .001) OR (TEST1 > .005) THEN 'Reinitialise for next cycle NN% = 1 XMC = 0 PS = P XMWS = XMW WW = 0 WC = 0 NST% = 1 NFIN% = NFI% NITE% = NITE% + 1 NO = 4 RELOOP = 1 ' gOTO 404 END IF END IF END IF WEND'RELOOP '307 The dimendionless pressures and works have been calcolated for one 'cycle. Now the additional heat and power losses will be calculated. 'Calculate average dimensionless pressure. PAVG = 0 '311 FOR I = 1 TO NDIV% PAVG = PAVG + PR(I) NEXT I PAVG = PAVG / NDIV% 'determine maximum and minimum dimensionless pressure PMAX = xlarge(PR(), NDIV%) PMIN = small(PR(), NDIV%) 'Adjust dimensionless works to relate to newly determend maximum 'pressure WC = WC / PMAX WW = WW / PMAX 'pressure ratio RP = PMAX / PMIN 'Find maximum masses and adjust them to maximum pessure. CMMAX = xlarge(XMCX(), NDIV%) WMMAX = xlarge(XMWX(), NDIV%) CMMAX = CMMAX / PMAX WMMAX = WMMAX / PMAX 'calc. max. pressure, mpa PMX = PMAX * PMX1 / PAVG 'Calculates angle between pressure wave and volume wave for a heat 'engine. ARG = 2 * RP / (RP - 1) * WW / PIE' 3.1416 'if (1-ARG**2) 1607,1608,1608 IF 1 - ARG ^ 2 < 0 THEN RETURN'gOTO 1607 ' FIPV = arsin(ARG) ' FIPU = ATN(SIN(ARG) / SQR(1 - SIN(ARG) ^ 2)) FIPV = ATN(ARG / (SQR(1 - ARG ^ 2))) XNDS = NDS% 'Calculates values used in flow loss calculations and flow integrals X = 0 DX = 1 / XNDS NIN% = NDS% + 1 COR = PMAX ^ (XNHT - 2) * DALF ^ (XNHT - 1) 'do 854 I=1,NIN% FOR I = 1 TO NIN% 'GOSUB PDINT CALL pdint(X, XDMW(), XDMC(), RVT, DC(), NDIV%, DMRE, PR(), XINT, DPR(), XI1, XI2, XNHT) XINT = XINT / DALF / PMAX DMRE = DMRE / PMAX / (2 * PIE)' 6.2832 XI1 = XI1 * COR / (1.5708 * DMRE) ^ (1 - XNHT) XI2 = XI2 * COR / (1.5708 * DMRE) ^ (2 - XNHT) XI3 = XI1 / XI2 GDMS(I) = DMRE GINT(I) = XINT GI2(I) = XI2 GI3(I) = XI3 X = X + DX NEXT I '854 continue 'Interpolates flow integrals FOR I = 1 TO 5 UIN(I) = plot(GINT(), H(I)) UDM(I) = plot(GDMS(), H(I)) NEXT I UI23 = plot(GI2(), H(2)) UI24 = plot(GI2(), H(4)) UI33 = plot(GI3(), H(2)) UI34 = plot(GI3(), H(4)) ' **** Calculation of constants ***** 'specific for hydrogen gas HMU = 8.873E-05 + .0000002 * (BTW - 293) CMU = 8.873E-05 + .0000002 * (BTC - 293) BTR = (BTW - BTC) / LOG(BTW / BTC) RMU = 8.873E-05 + .0000002 * (BTR - 293) CP1 = 14.6 CV1 = 10.46 R2 = 8.23168E+07 R = 4.116 '373 ******** Cold exchanger pressure drop ********** REC = UDM(1) * PMX * SPD * HAC * CTD / (BTC * AFC * CMU * R) IF (REC <= 2000) THEN FC = 16 / REC ELSE FC = EXP(-1.34 - .2 * LOG(REC)) END IF GLS = CTLL * SPD * SPD * HAC * HAC * FC * UIN(1) / (CTD * AFC * AFC * BTC * R2) QP = NOC% * SPD * PMX * HAC / (2 * PIE) QCP = QP * GLS '********** Hot exchanger pressure drop *********** REH = UDM(5) * PMX * SPD * HAC * HTD / (BTW * AFH * HMU * R) '384 if(REH-2000) 1988,1988,1989 IF (REH <= 2000) THEN FH = 16 / REH ELSE FH = EXP(-1.34 - .2 * LOG(REH)) END IF GLH = HTLL * SPD ^ 2 * HAC ^ 2 * BTW * FH * UIN(5) / (HTD * AFH ^ 2 * BTC ^ 2 * R2) QHP = QP * GLH '**** Screen---netnet option ********* RER = PMX * HAC * SPD * BWD / (AFR * R) RE(1) = RER * UDM(2) / (BTC * CMU) RE(2) = RER * UDM(3) / (BTR * RMU) RE(3) = RER * UDM(4) / (BTW * HMU) FOR I = 1 TO 3 IF (NET%) <= 0 THEN IF (RE(I) <= 60) THEN FR(I) = EXP(1.73 - .93 * LOG(RE(I))) ELSE IF (RE(I) <= 1000) THEN FR(I) = EXP(.714 - .365 * LOG(RE(I))) ELSE FR(I) = EXP(.015 - .125 * LOG(RE(I))) END IF END IF ELSE FR(I) = 2.73 * (1 + 10.397 / RE(I)) END IF NEXT I 'continue '407 ******* regenerator pressure drop ***** GLR = BRL * SPD * SPD * HAC * HAC / (BWD * AFR * AFR * R2 * BTC) QR1 = QP * GLR * UIN(2) * FR(1) QR2 = QP * GLR * UIN(3) * FR(2) * BTR / BTC QR3 = QP * GLR * UIN(4) * FR(3) * UTR QRP = (QR1 + QR3 + 4 * QR2) / 6 'Calculates effective hot and cold gas temperatures based upon the 'number of transfer units in the heat exchangers, specific for 'hydrogen CNTU = .112 * CTLS / (CTD * REC ^ .2) DTC = WC * (SHR - 1) / (2 * UDM(1) * SHR * (EXP(2 * CNTU) - 1)) BTC = BTC1 * (1 - DTC) HNTU = .1044 * HTLS / (HTD * REH ^ .2) DTH = WW * (SHR - 1) / (2 * UDM(5) * SHR * (EXP(2 * HNTU) - 1)) BTW = BTW1 * (1 - DTH) 'Note, temperature ratio is redefined for next iteration UTR = BTW / BTC '424 *** Reheat loss **** RNTU = BRL * 4.37 / (BWD * SQR(PI4 * 2 * RE(1))) QNTU = BRL * 4.031 / (BWD * SQR(PI4 * 2 * RE(3))) QNPH = AFR * BRL * .195 / (PI4 * HAC * UDM(2) * (UTR - 1)) QDK = QNPH * (UI33 + UI34 * UDM(2) / UDM(4)) / 2 QLM = (1 + QDK) / (RNTU / UI23 + QNTU * UDM(2) / (UDM(4) * UI24)) QHR = UDM(2) * CP1 * (BTW - BTC) * SPD * PMX * HAC * QLM * NOC% / (R * BTC * 2) '**** shuttle loss ***** QL1 = 231.2 * SQR(SPD * BRC * BRC) QB = (2 * QL1 * QL1 - QL1) / (2 * QL1 * QL1 - 1) QHC = .00146 * BST * (BTW - BTC) * PI4 * BPD * BST * QB * NOC% / (BRC * BPL) '**** Pumping loss ******* IF (MS% AND 1) = 0 THEN QFS = (RP / (BTW / (BTW - 2 * BTC) - BST / BPL)) + (1 / (BTW / ((BTW - 2 * BTC) + BST / BPL))) ELSE ' Martini writes in the manual the calculation above is wrong !!!!!!!!!!!!!! ' has to be X = (BTW + BTC1) / (BTW - BTC1) Y = BST / BPL QFS = -RP / (X + Y) + 1 / (X - Y) END IF QHG = ABS(SPD * PMX * GGV * BST * SHR * QFS * ARG * NOC% / ((SHR - 1) * BPL * RP * 8)) '**** Basic powers ***** 'PO = (WW * HAV + WC * HAC) * (4.5) * PMX * SPD * NOC% / PIE PO = (WW * HAV + WC * HAC) * (.5) * PMX * SPD * NOC% / PIE '441 ****** Net power ******** POF = PO - QCP - QHP - QRP 'Get ready to report on one iteration and prepare for the next. 'Reset hot end dimensionless heat transfer integral HTW = 0 'The program tries to keep PMAX=1. This adjustment of the pressure 'and mass does this '448 do 509 I=1,NDIV% FOR I = 1 TO NDIV% PR(I) = PR(I) / PMAX XMCX(I) = XMCX(I) / PMAX XMWX(I) = XMWX(I) / PMAX NEXT I '452 Dimensionless hot and cold gas temperatures for esch increment. 'if they are less than zero correct to zero. WI(NDIV1%) = WI(1) CI(NDIV1%) = CI(1) FOR I = 1 TO NDIV% IF (XMCX(I)) > 0 THEN TEC(I) = PR(I) * CI(I + 1) / XMCX(I) ELSE TEC(I) = 0 END IF IF (XMWX(I)) > 0 THEN TEW(I) = PR(I) * WI(I + 1) / XMWX(I) ELSE TEW(I) = 0 END IF NEXT I '466 Dimensionless average hot and cold gas temperatures for ful cycle. TEW(NDIV1%) = TEW(1) TEC(NDIV1%) = TEC(1) PR(NDIV1%) = PR(1) XMCX(NDIV1%) = XMCX(1) XMWX(NDIV1%) = XMWX(1) TWDM = 0 TCDM = 0 '474 do 573 I=1,NDIV% FOR I = 1 TO NDIV% DMW = XMWX(I + 1) - XMWX(I) IF (DMW) < 0 THEN TMPW = (TEW(I) + TEW(I + 1)) / 2 TWDM = TWDM + (TMPW - 1) * DMW END IF DMC = XMCX(I + 1) - XMCX(I) IF (DMC) < 0 THEN TMPC = (TEC(I) + TEC(I + 1)) / 2 TCDM = TCDM + (TMPC - 1) * DMC END IF NEXT I TWDM = TWDM * SHR / (SHR - 1) TCDM = TCDM * SHR / (SHR - 1) 'Hot end heat transfer integral for full cycle and total gas mass at 'each point in the cycle. Total mass should not change. '408 do 1021 I=1,NDIV% FOR I = 1 TO NDIV% HTW = HTW + (WI(I + 1) - WI(I)) * (PR(I) + PR(I + 1)) / 2 1021 XMT(I) = XMCX(I) * RVT + XMWX(I) + PR(I) * VD NEXT I 'Basic heat input, watts HT = HTW * SPD * PMX * HAV * NOC% / (2 * PIE) 'Specific static condution heat loss for the 4l23 engine CON = 9680 'put as input 'flow friction credit, watts FFF = (QHP + .5 * QRP) * (-1) 'Heat to engine, watts HTE = HT + QHR + QHC + QHG + CON + FFF 'indicated efficiency, % ZEF = 100 * POF / HTE 'print out results of one iteration 'write(6,12) LUP% 'write(6,3010)PO,HT 'write(6,3020) QHP,QHR 'write(6,1925) QRP,QHC,QCP,QHG,POF,CON,ZEF,FFF,HTE 'write(6,1921) BTW,BTC,RVT,VD PRINT "Iteration "; lup PRINT "Basic power (watts)="; PO, PRINT "Basic heat (watts)="; HT 'PRINT PRINT "Heater windage (watts)="; QHP, PRINT "Reheat loss (watts)="; QHR PRINT "Regen. windage (watts)="; QRP, PRINT "Shuttle loss (watts) ="; QHC PRINT "Cooler windage (watts)="; QCP, PRINT "Appendix loss (watts)="; QHG PRINT "Net power (watts) ="; POF, PRINT "Conduction (watts) ="; CON PRINT "Indicated efficiency(%)"; ZEF, PRINT "Flow friction (watts) ="; FFF PRINT , , , "Heat to engine (watts)="; HTE PRINT "Effective hot temperature (k)="; BTW, PRINT "Effective cold temperature(k)="; INT(BTC * 10) / 10 PRINT "Displaced mass ratio ="; RVT, PRINT "Reduced dead volume ="; VD INPUT I$ GOSUB printresult 'After all losses are taken into account LUP% is indexed. The program 'does 3 iterations with printouts before going into a summary 'LUP% = LUP% + 1 '510 if(LUP%-3) 339,339,1607 NEXT lup 'if input volume NWR% is other than zero the following summary 'Information is printed at the end of the computation RETURN printout: '1607 'IF (NWR%) = 0 THEN gOTO 242 '1613 write(6,51)TWDM,TCDM PRINT "Dimensionless average bas temp. hot end ="; TWDM PRINT "Dimensionless average bas temp. cold end="; TCDM INPUT I$ 'Print out each 10 degrees, angle, hot volume, hot gas temperature. 'cold volume, cold gas temperature, total volume, pressure '1149 write(6,20) '518 do 3001 I=10,NDIV%,10 xmaximum = 0 vtmaximum = 0 CLS PRINT PRINT " Pressure,MPA-volume,CN3 data for one cylinder" PRINT " ==============================================================================" PRINT "angle de. cold vol.cm^3 pressure MPa "; PRINT "angle de. cold vol.cm^3 pressure MPa" PRINT " hot vol.cm^3 total vol.cm^3"; PRINT " hot vol.cm^3 total vol.cm^3" FOR I = 10 TO NDIV% STEP 10 IF I <= NDIV% / 2 THEN hpos = 1 vpos = 5 + I / 10 ELSE hpos = 41 vpos = 5 + (I - NDIV% / 2) / 10 END IF X = PR(I) * PMX VH = VHD + HAV * WI(I) VC = VCD + CI(I) * HAC VT = VH + VRD + VC IF VT > vtmaximum THEN vtmaximum = VT IF X > xmaximum THEN xmaximum = X LOCATE vpos, hpos IF I < 100 THEN PRINT " "; PRINT I; CALL makestr(VH, 8, 2, printv$): PRINT printv$; CALL makestr(VC, 8, 2, printv$): PRINT printv$; CALL makestr(VT, 8, 2, printv$): PRINT printv$; CALL makestr(X, 8, 2, printv$): PRINT printv$; ' IF I = 140 THEN INPUT I$ NEXT I PRINT #OUTFILE%, PRINT #OUTFILE%, "Pressure,MPA-volume,CN3 data for one cylinder" PRINT #OUTFILE%, "angle de. hot vol.cm^3", PRINT #OUTFILE%, "cold vol.cm^3", PRINT #OUTFILE%, "total vol.", PRINT #OUTFILE%, "pressure MPa" FOR I = 10 TO NDIV% STEP 10 X = PR(I) * PMX VH = VHD + HAV * WI(I) VC = VCD + CI(I) * HAC VT = VH + VRD + VC IF I < 100 THEN PRINT #OUTFILE%, " "; PRINT #OUTFILE%, I; CALL makestr(VH, 13, 7, printv$): PRINT #OUTFILE%, printv$, CALL makestr(VC, 13, 7, printv$): PRINT #OUTFILE%, printv$, CALL makestr(VT, 13, 7, printv$): PRINT #OUTFILE%, printv$, CALL makestr(X, 13, 9, printv$): PRINT #OUTFILE%, printv$ NEXT I PRINT '524 Starts over with next DATA SET RETURN printinput: PRINT #OUTFILE%, PRINT #OUTFILE%, "Cold metal temp (k) ="; BTC PRINT #OUTFILE%, "Hot metal temp (k) ="; BTW PRINT #OUTFILE%, "Piston diameter(cm) ="; BPD PRINT #OUTFILE%, "Piston stroke (cm) ="; BST PRINT #OUTFILE%, PRINT #OUTFILE%, "Hot cap length (cm) ="; BPL PRINT #OUTFILE%, "Piston gap (cm) ="; BRC PRINT #OUTFILE%, "Piston end clearence(cm)="; BEC PRINT #OUTFILE%, "Piston rod diameter (cm)="; PDR PRINT #OUTFILE%, PRINT #OUTFILE%, "Regenerator porosity ="; BRO PRINT #OUTFILE%, "Regenerator diameter(cm)="; BDR PRINT #OUTFILE%, "Regenerator wire dia(cm)="; BWD PRINT #OUTFILE%, "Regenerator length (cm)="; BRL PRINT #OUTFILE%, PRINT #OUTFILE%, "Number of regenerators ="; MBR% PRINT #OUTFILE%, "Number of cooler tubes ="; MCT% PRINT #OUTFILE%, "Cool tube diameter (cm) ="; CTD PRINT #OUTFILE%, "Total cool tube length(cm)="; CTLL PRINT #OUTFILE%, PRINT #OUTFILE%, "Cool ct length (cm)="; CTLS PRINT #OUTFILE%, "Heater tube diameter (cm)="; HTD PRINT #OUTFILE%, "Total heater tube length(cm)="; HTLL PRINT #OUTFILE%, "Heat heater tube length (cm)="; HTLS PRINT #OUTFILE%, "Number of heater tubes ="; MHT% PRINT #OUTFILE%, "Heat transfer exponent ="; XNHT PRINT #OUTFILE%, PRINT #OUTFILE%, "Average pressure MPA ="; PMX1 PRINT #OUTFILE%, "Engine speed (rad/sec)="; SPD PRINT #OUTFILE%, "NETNET option ="; NET% PRINT #OUTFILE%, "Number of cylinders ="; NOC% PRINT #OUTFILE%, PRINT #OUTFILE%, "Cold crank ratio ="; ZZC PRINT #OUTFILE%, "Hot crank ratio ="; ZZW RETURN printresult: PRINT #OUTFILE%, PRINT #OUTFILE%, "Iteration "; lup PRINT #OUTFILE%, "Basic power (watts)="; PO, PRINT #OUTFILE%, "Basic heat (watts)="; HT 'PRINT PRINT #OUTFILE%, "Heater windage (watts)="; QHP, PRINT #OUTFILE%, "Reheat loss (watts)="; QHR PRINT #OUTFILE%, "Regen. windage (watts)="; QRP, PRINT #OUTFILE%, "Shuttle loss (watts) ="; QHC PRINT #OUTFILE%, "Cooler windage (watts)="; QCP, PRINT #OUTFILE%, "Appendix loss (watts)="; QHG PRINT #OUTFILE%, "Net power (watts) ="; POF, PRINT #OUTFILE%, "Conduction (watts) ="; CON PRINT #OUTFILE%, "Indicated efficiency(%)"; ZEF, PRINT #OUTFILE%, "Flow friction (watts) ="; FFF PRINT #OUTFILE%, , , , "Heat to engine (watts)="; HTE PRINT #OUTFILE%, "Effective hot temperature (k)="; BTW, PRINT #OUTFILE%, "Effective cold temperature(k)="; INT(BTC * 10) / 10 PRINT #OUTFILE%, "Displaced mass ratio ="; RVT, PRINT #OUTFILE%, "Reduced dead volume ="; VD RETURN 9900 'read variable text strings RESTORE 9902'variable text READ ZZC$, ZZW$, XNHT$, PMX1$, SPD$, NET$, NOC$ READ SHR$, NFI$, NDIV$, NWR$, NDS$, DDD$, DLL$ READ BTC$, BTW$, BPD$, BST$, BPL$, BRC$, BEC$ READ PDR$, BRO$, BDR$, BWD$, BRL$, MBR$, MCT$ READ CTD$, CTLL$, CTLS$, HTD$, HTLL$, HTLS$, MHT$ READ CON$ RETURN 9902 'variable text DATA "ZZC = (input) Connecting rod length/.5 stroke for cold piston" DATA "ZZW = (input) Connecting rod length/.5 stroke for hot piston" DATA "XNHT = (input) Value for exponent in heat transfer relation of regenerator matrix" DATA "PMX1 = (input) Avg. pressure MPa" DATA "SPD = (input) Engine speed, rad/sec" DATA "NET% = (input) Regenerator filler option (+5 = metnet, -5 = screen)" DATA "NOC% = (input) number of cylinders" DATA "SHR = (input) Specific heat ratio for working gas" DATA "NFI% = (input) (phase angle)(NDIV%)/360" DATA "NDIV% = (input) Number of divisions per crank rotation (must be a multiple of 4)" DATA "NWR% = (input) Governs printout, zero for overall results only, different from zero added pv data" DATA "NDS% = (input) Number of divisions in dead space" DATA "DDD = (input) Cooler duct diameter, cm" DATA "DLL = (input) Cooler duct length, cm" DATA "BTC = (input) Effective cold temperature, K" DATA "BTW = (input) Hot effective gas temperature, K" DATA "BPD = (input) Piston diameter, cm" DATA "BST = (input) Piston stroke, cm" DATA "BPL = (input) Hot cap length, cm" DATA "BRC = (input) Piston gap, cm" DATA "BEC = (input) Piston end clearence, cm" DATA "PDR = (input) Piston rod diameter, cm" DATA "BRO = (input) Regenerator density factor" DATA "BDR = (input) Regenerator diameter, cm" DATA "BWD = (input) Effective regenerator wire diameter, cm" DATA "BRL = (input) Regenerator length, cm" DATA "MBR% = (input) Number of regenerators" DATA "MCT% = (input) Number of cooler tubes per cylinder" DATA "CTD = (input) Cooler tube inside diameter, cm" DATA "CTLL = (input) Total cooler tube length, cm" DATA "CTLS = (input) Cooler cool tube length, cm" DATA "HTD = (input) Heater tube inside diameter, cm" DATA "HTLL = (input) Total heater tube length, cm" DATA "HTLS = (input) Heated heater tube length, cm" DATA "MHT% = (input) Number of heater tubes for cylinder" DATA "CON = Conduction loss, watts" 9999 'input DATA "GM 4L23" DATA 5.874, 5.874, .204 , 17.9263, 157.080, 5 , 4 'DATA 5.874, 5.874, .204 , 9.6516, 209.44 , 5 , 4 DATA 1.39 , 90 , 360 , 1 , 10 , .76 , 71 DATA 330 , 1033.15, 10.16, 4.65 , 6.4 , .0406, .0406 'DATA 330 , 1033 , 10.16, 4.65 , 6.4 , .0406, .0406 DATA 4.060, .8 , 3.5 , .0043 , 2.5 , 6 , 312 DATA .115 , 12.9 , 12.02, .472 , 41.8 , 25.58, 36, 0.0 'CON=specific static conduction heat loss 9680 for the 4L23 engine DATA 9680 DATA end 'integer data must be right justified. There are 7 data fields per 'line of 7 columns each 'data field layout is as follows 'ZZC, ZZW, XNHT, PMX1, SPD, NET%, NOC% 'SHR, NFI%, NDIV%, NWR%, NDS%, DDD, DLL 'BTC, BTW, BPD, BST, BPL, BRC, BEC 'PDR, BRO, BDR, BWD, BRL, MBR%, MCT% 'CTD, CTLL, CTLS, HTD, HTLL,HTLS,MHT% 'Next data set repeats last 'change according to data printout. 'Data must be within control cards(//), 'read(5.5) ZZC,ZZW,XNHT,PMX1,SPD,NET%,NOC% 'read (5+1) SHR,NFI%,NDIV%,NWR%,NDS%,DDD,DLL 'read (5+1710) BTC,BTW,BPD,BST,BPL,BRC,BEC 'read (5+1720) PDR,BRO,BDR,BWD,BRL,MBR%,MCT% 'read (5+1730) CTD,CTLL,CTLS,HTD,HTLL,HTLS,MHT% 2000 'input variable names DATA BTC, BTW,BPD,BST,BPL, BRC,BEC,PDR DATA BRO,BDR,BWD,BRL,MBR%,MCT%,CTD,CTLL DATA CTLS,HTD,HTLL,HTLS,MHT%,XNHT DATA PMX1,SPD,NET%, NOC%,ZZC,ZZW DATA SHR, NDIV%, FI1, NDS%, DDD, DLL DATA CON 2002 'inpus abreviations DATA "Cold metal temp (k) =","Hot metal temp (k) =" DATA "Piston diameter(cm) =","Piston stroke (cm) =" DATA "Hot cap length (cm) =","Piston gap (cm) =" DATA "Piston end clearence(cm)=","Piston rod diameter (cm)=" DATA "Regenerator porosity =","Regenerator diameter(cm)=" DATA "Regenerator wire dia(cm)=","Regenerator length (cm)=" DATA "Number of regenerators =","Number of cooler tubes =" DATA "Cool tube diameter (cm) =","Total cool tube length(cm)=" DATA "Cool ct length (cm)=","Heater tube diameter (cm)=" DATA "Total heater tube length(cm)=","Heat heater tube length (cm)=" DATA "Number of heater tubes =","Heat transfer exponent =" DATA "Average pressure MPA =","Engine speed (rad/sec)=" DATA "NETNET option =","Number of cylinders =" DATA "Cold crank ratio =","Hot crank ratio =" DATA "Specific heat ratio =","Div. per cycle =","Phase angle (deg.) =" DATA "Incr. in dead space int.=","Duct diameter (cm) =","Duct length (cm)=" DATA "SPECIFIED CONDUCTION=" 2006 'init load inputs abreviationsni$() RESTORE 2000 FOR var = 1 TO nivar READ ni$(var) NEXT var RETURN SUB closefilenr (filecount%, filenr%) 'CALL closefilenr (filecount%,filenr%) 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 SUB givefilenr (filecount%, filenr%) 'call givefilenr (filecount%, filenr%) IF filecount% = 0 THEN filecount% = 1 filenr% = 0 WHILE (filecount% OR 2 ^ filenr%) = filecount% filenr% = filenr% + 1 WEND filecount% = filecount% + 2 ^ filenr% END SUB SUB makestr (printv, length, decimals, printv$) 'call makestr(printv,length,decimals,printv$) 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 SUB pdint (X, DMW(), DMC(), RVT, DVC(), NDIV%, DM, PR(), XINT, DPR(), XI1, XI2, XNHT) ' CALL pdint(X, XDMW(), XDMC(), RVT, DC(), NDIV%, DMRE, PR(), XINT, DPR(), XI1, XI2, XNHT) 'dimension DMW(720), DMC(720), DVC(720), PR(720), DPR(720) '649 Subroutine for pressure drop integral calculation DM = 0 XINT = 0 XI1 = 0 EX1 = 1 - XNHT XI2 = 0 EX2 = 2 - XNHT FOR I = 1 TO NDIV% DMX = DMC(I) - X * (DMW(I) / RVT + DMC(I)) Y = ABS(DMX) DM = DM + Y A = DPR(I) * Y ^ EX1 IF (DMX) < 0 THEN A = -A XI1 = XI1 + A XI2 = XI2 + Y ^ EX2 101 XINT = XINT + Y * DMX / PR(I) * DVC(I) NEXT I 'xndiv = NDIV% 'XNDIV KOMT VERDER NIET VOOR END SUB FUNCTION plot (X(), H) 'function TO INTERPOLATE '599 subroutine to ................. 'dimension X(20) 'dim X(20) M% = 10 * H Z = H - M% / 10 M% = M% + 1 G1 = X(M%) N% = M% + 1 G2 = X(N%) plot = Z * G2 + (1 - Z) * G1 END FUNCTION FUNCTION small (X(), NDIV%) 'subroutine to find smallest of a list smallnow = X(1) FOR I = 2 TO NDIV% IF smallnow > X(I) THEN smallnow = X(I) END IF NEXT I small = smallnow END FUNCTION SUB VOLC (DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(), COFI(), SALF(), CALF(), PDR, ALF) 'subroutine to list cold volumes and derivatives 'CALL VOLC(DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(), 'COFI(), SALF(), CALF(), PDR, ALF) NDIV1% = NDIV% + 1 FOR I = 1 TO NDIV1% COFI(I) = COS(ALF) SALF(I) = SIN(ALF) ALF = ALF + DALF NEXT I FOR I = 1 TO NDIV% CALF(I) = (SALF(I + 1) - SALF(I)) / DALF NEXT I CALF(NDIV1%) = CALF(1)'Richard has NDIVI here 'CALF(NDIVI) = CALF(1)'Richard has NDIVI here FOR I = 1 TO NDIV% SIFI(I) = SALF(I) SALF(I) = (SALF(I) + SALF(I + 1)) / 2 NEXT I COFI(NDIV1%) = COFI(1) SIFI(NDIV1%) = SIFI(1) N% = NF% * 4 FOR I = 1 TO N% CRC = SQR(ZZC ^ 2 - CALF(I) ^ 2) ' see note #1 IF (PDR) > 0 THEN C(I) = 1 - SALF(I) + CRC - ZZC CI(I) = 1 - SIFI(I) + CRC - ZZC DC(I) = -CALF(I) * (1 - SALF(I) / CRC) DCI(I) = -CALF(I) * (1 - SIFI(I) / CRC) ELSE C(I) = 1 + SALF(I) - CRC + ZZC CI(I) = 1 + SIFI(I) - CRC + ZZC DC(I) = CALF(I) * (1 - SALF(I) / CRC) DCI(I) = CALF(I) * (1 - SIFI(I) / CRC) END IF NEXT I END SUB SUB VOLW (W(), WI(), DW(), DWI(), CFI, SFI, ZZW, NDIV%, SIFI(), COFI(), SALF(), CALF(), DALF) 'subroutine for pressure drop integral calculation 'subroutine to list hot volumes and derivates 'CALL VOLW(W(), WI(), DW(),DWI(),CFI,SFI,ZZW,NDIV%,SIFI(),COFI(),SALF(),CALF(),DALF) SIFIP = SIFI(1) * CFI - COFI(1) * SFI FOR I = 1 TO NDIV% SALF1 = SIFI(I + 1) * CFI - COFI(I + 1) * SFI SALFP = (SIFIP + SALF1) / 2 CALFP = (SALF1 - SIFIP) / DALF CRW = SQR(ZZW ^ 2 - CALFP ^ 2) W(I) = 1 + SALFP - CRW + ZZW WI(I) = 1 + SIFIP - CRW + ZZW DW(I) = CALFP * (1 - SALFP / CRW) DWI(I) = CALFP * (1 - SIFIP / CRW) SIFIP = SALF1 NEXT I END SUB FUNCTION xlarge (X(), NDIV%) 'subroutine to find largest of a list xlargenow = X(1) FOR I = 2 TO NDIV% IF xlargenow < X(I) THEN xlargenow = X(I) END IF NEXT I xlarge = xlargenow END FUNCTION '==================================================================