Hi,
I am getting a run-time error when trying to compile my program and as I
am not very experienced in Fortran I am not sure what the error is or how
to fix it. The program has been lifted out of a book called "Interneal
Combustion Engines:Applied Thermosciences" and is a program to model the
heat release rate of a an engine. It requires a couplke of subroutines,
DVERK and LEQIF which are IMSL routines and I do not have. I have found
the code for DVERK but unfortunately have not found the code for LEQIF
which is why I assume I am receiving the run-time error. the error is a
zero or negative argument to a logarithm routine.
FARG - in file release2_good.for at line 664 [+1738]
main - in file release2_good.for at line 41 [+0190]
The main code is listed below and any help in resolving this matter would be greatly appreciated. Or if anybody else has a heat release program that they are willing to share that would be great too.
I can e-mail the code if the formatting is not retained in this thread to make things easier.
Thanks
Chris
C GIVEN IN BLOCK DATA SUBPROGRAM
C COMPRESSION RATIO - R
C BORE - B (CM)
C STROKE - S (CM)
C HALF-STROKE TO ROD RATIO - EPS
C ENGINE SPEED - RPM
C HEAT TRANSFER COEFFICIENT - H (J/S/M**2/K)
C BLOWBY COEFFICIENT - C (L/S)
C EQUIVALENCE RATIO - PHI
C RESIDUAL FRACTION - F
C INITIAL PRESSURE - P1 (BAR)
C INITIAL TEMPERATURE - T1 (K)
C WALL TEMPERATURE - TW (K)
C BURN ANGLE - THETAB (DEG)
C START OF HEAT RELEASE - THETAS (DEG ATDC)
C
C FIND:
C INDICATED THERMAL EFFICIENCY - ETA
C INDICATED MEAN EFFECTIVE PRESSURE - IMEP (BAR)
C
C NOTES:
C 1. COSINE BURNING LAW IS EMPLOYED, EQN. 4.61
C ` 2. FUEL - GASOLINE (C7H17)
C
EXTERNAL RATES
REAL IMEP, MNOT, M
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT, UFINAL
DIMENSION Y1(6), Y2(10), Y(6), COM(24), WM(6,9)
DATA AO/47870./, FS/.06548/, Y/1., 10000., 350., 3*0.0/,
& N/6/, THETA/-180./, THETAE/-179./, TOL/.001/, IND/1/,
& NW/6/, PI/3.141593/, Y2/10*0.0/
C
C LABEL HEADER ON OUTPUT PAGE AND PRINT INITIAL CONDITIONS
C
OMEGA = RPM*PI/30
CALL AUXLRY(THETA, V, X, EM)
Y(1) = P1
Y(3) = T1
CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, SU, Y1, CP, DLVLT,
& DLVLP)
MNOT = V/VU
UNOT = MNOT*U
WRITE(6,10)
10 FORMAT(5X,'THETA',4X,'VOLUME',3X,'BURN FRAC',2X,'PRESSURE',3X,
& 'BURN TEMP',2X,'T UNBURN',3X,'WORK',3X,'HEAT LOSS',5X,
& 'MASS',3X,'H-LEAK')
WRITE(6,15)
15 FORMAT(4X,'DEG ATDC',2X,'CM**3',6X,'--',8X,'BARS',7X,'KELVIN',
& 5X,'KELVIN',3X,'JOULES',4X,'JOULES',6X,'GRAMS',2X,
& 'JOULES')
WRITE(6,20) THETA, V, X, (Y(I),I=1,5), MNOT, Y(6)
20 FORMAT(4X,F5.0,6X,F5.0,4X,F5.3,6X,F5.2,7X,F5.0,5X,F5.0,
& 5X,F5.0,3X,F5.0,7X,F5.3,3X,F5.2)
C
C INTEGRATE THE RATE EQUATIONS 4.80 - 4.85 AND PRINT ANSWERS
C EVERY 10 DEGREES
C
DO 30 I = 1,36
DO 25 J = 1,10
CALL DVERK (N, RATES, THETA, Y, THETAE, TOL, IND, COM, NW, WM)
CALL AUXLRY(THETAE, V, X, EM)
M = EM*MNOT
THETA = THETAE
THETAE = THETA + 1.
IF(THETAS .GE. THETA .AND. THETAS .LT. THETAE)
& CALL TINITL(Y(1), T(3), PHI, F, Y(2))
IF(THETA .GT. THETAS + THETAB) Y(3) = 10000.
25 CONTINUE
WRITE(6,20) THETA, V, X, (Y(J),J=1,5), M, Y(6)
30 CONTINUE
C
C COMPUTE AND WRITE THE EAT, IMEP AND ERRORS
C
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, SB, Y2, CP, DLVLT, DLVLP,
& IER)
UFINAL = U*M
ERROR1 = 1.0 - VB*M/V
ERROR2 = 1.0 + Y(4)/(UFINAL - UNOT + Y(5) + Y(6))
ETA = Y(4)/MNOT*(1. + PHI*FS*(1. - F))/PHI/FS/(1. - F)/AO
IMEP = Y(4)/(PI/4.*B**2*S)*10.
WRITE(6,40) ETA, IMEP, ERROR1, ERROR2
40 FORMAT('ETA = ',1PE11.4,' IMEP = ',1PE11.4,'ERROR1 = ',1PE11.4,
& ' ERROR2 = ',1PE11.4)
STOP
END
C
SUBROUTINE RATES(N, THETA, Y, YPRIME)
DIMENSION Y1(6), Y2(10)
REAL Y(N), YPRIME(N), M, MNOT
COMMON /ENGINE/ R, BORE, STROKE, EPS, RPM, HEAT, CEE, THETAB,
& THETAS, PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT,
& UFINAL
DATA PI/3.141593/, Y2/10*0.0/
CALL AUXLRY(THETA, V, X, EM)
M = EM*MNOT
RAD = THETA*PI/180.
DUMB = SQRT(1. - (EPS*SIN(RAD))**2)
DV = PI/8.*BORE**2*STROKE*SIN(RAD)*(1. + EPS*COS(RAD)/DUMB)
A = (DV + V*CEE/OMEGA)/M
C1 = HEAT*(PI*BORE**2/2. + 4.*V/BORE)/OMEGA/M/10000.
C0 = SQRT(X)
C
C THREE DIFFERENT COMPUTATIONS ARE REQUIRED DEPENDING
C UPON THE SIZE OF THE MASS FRACTION BURNED
C
IF( X .GT. .999 ) GO TO 20
IF( X .GT. .001 ) GO TO 10
C
C COMPRESSION
C
CALL FARG(Y(1), Y(3), PHI, F, HL, U, VU, S, Y1, CP, DLVLT,
& DLVLP)
B = C1*VU/CP*DLVLT*(1.0 - TW/Y(3))
C = 0.
D = 0.
E = VU**2/CP/Y(3)*DLVLT**2 + 10.*VU/Y(1)*DLVLP
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = 0.
YPRIME(3) = -C1/CP*(Y(3) - TW) + VU/CP*DLVLT*YPRIME(1)/10.
GO TO 30
C
C COMBUSTION
C
10 CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, S, Y1, CPU, DLVLTU,
& DLVLP)
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, S, Y2, CPB, DDLVLTB,
& DLVLPB, IER)
B = C1*(VB/CPB*DLVLTB*C0*(1. - TW/Y(2)) + VU/CPU*DLVLTU*
& (1. - C0)*(1.- TW/Y(3)))
DX = 0.5*SIN(PI*(THETA - THETAS)/THETAB)*180./THETAB
C = -(VB - VU)*DX - VB/Y(2)*DLVLTB*(HU - HB)/CPB*(DX - (X -
& X**2)*CEE/OMEGA)
D = X*(VB**2/CPB/Y(2)*DLVLTB**2 + 10.*VB/Y(1)*DLVLPB)
E = (1. - X)*(VU**2/CPU/Y(3)*DLVLTU**2 + 10.*VU/Y(1)*DLVLPU)
HL = (1.0 - X**2)*HU + X**2*HB
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CPB/C0*(Y(2) - TW) + VB/CPB*DLVLTB*
& YPRIME(1)/10. + (HU - HB)/CPB*(DX/X - (1. - X)*
& CEE/OMEGA)
YPRIME(3) = -C1/CPU/(1. + CO)*(Y(3) - TW) + VU/CPU*DLVLTU*
& YPRIME(1)/10.
GO TO 30
C
C EXPANSION
C
20 CALL ECP(Y(1), Y(2), PHI, HL, U, VB, S, Y2, CP, DLVLT, DLVLP,
& IER)
B = C1*VB/CP*DLVLT*(1. - TW/Y(2))
C = 0.
D = VB**2/CP/Y(2)*DLVLT**2 + 10.*VB/Y(1)*DLVLP
E = 0.
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CP*(Y(2) - TW) + VB/CP*DLVLT*YPRIME(1)/10.
YPRIME(3) = 0.
C
C ALL CASES
C
30 YPRIME(4) = Y(1)*DV/10.
YPRIME(5) = C1*M*(C0*(Y(2) - TW) + (1. - C0)*(Y(3) - TW))
YPRIME(6) = CEE*M/OMEGA*HL
DO 40 I = 1,N
40 YPRIME(I) = YPRIME(I)*PI/180.
RETURN
END
C
SUBROUTINE AUXLRY(THETA, V, X, EM)
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA PI/3.141593/
RAD = THETA*PI/180.
VTDC = PI/4.*B**2*S/(R - 1.)
V = VTDC*(1. + (R - 1.)/2.*(1. - COS(RAD) + 1./EPS*(1. -
& SQRT(1. - (EPS*SIN(RAD))**2))))
X = 0.5*(1. - COS(PI*(THETA - THETAS)/THETAB))
IF(THETA .LE. THETAS) X = 0.
IF(THETA .GE. THETAS + THETAB) X = 1.
EM = EXP(-C*(RAD + PI)/OMEGA)
RETURN
END
C
BLOCK DATA
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA R/10./, B/10./, S/8./, EPS/.25/, RPM/2000./, H/500./,
& C/0.8/, THETAB/60./, THETAS/-35./, PHI/0.8/, F/.1/,
& P1/1.0/, T1/350./, TW/420./
C
END
C
SUBROUTINE TINITL(P, TU, PHI, F, TB)
DIMENSION Y1(6), Y2(10)
DATA MAXITS/50/, TOL/.001/, Y2/10*0.0/
TB = 2000.
CALL FARG(P, TU, PHI, F, HU, U, V, S, Y1, CP, DLVLT, DLVLP)
DO 10 I = 1,MAXITS
CALL ECP(P, TB, PHI, HB, U, V, S, Y2, CP, DLVLT, DLVLP, IER)
DELT = (HU - HB)/CP
TB = TB + DELT
IF(ABS(DELT/TB) .LT. TOL) RETURN
10 CONTINUE
WRITE(6,20)
20 FORMAT('CONVERGENCE FAILURE IN TINITL')
RETURN
END
SUBROUTINE ECP (P, T, PHI, H, U, V, S, Y, CP, DLVLT, DLVLP, IER)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF EQUILIBRIUM COMBUSTION
C PRODUCTS
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE )K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/G)
C S - ENTROPY (J/G/K)
C Y - A TEN DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2. 2=H2O, 3=N2, 4=O2, 5=CO, 6=H2
C 7=H, 8=O, 9=OH, 10=NO
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C IER - AN ERROR MESSAGE
C 0 - NORMAL
C 1 - CONVERGENCE FAILURE: COMPOSITION LOOP
C 2 - CALLED WITH TOO HIGH AN EQUIVALENCE RATIO; C(S)
C AND OTHER SPECIES WOULD FORM
C
REAL M, MW, K, KP, MT, MP
LOGICAL CHECK
DIMENSION M(10), K(6), C(6), Y(10), D(3), B(4), A(4,4), Y0(6),
& KP(5,6), DFDT(4), DCDT(6), DKDT(6), DCDP(6),
& DYDT(10), DYDP(10), A0(7,10), DFDP(4), CP0(10),
& H0(10), S0(10), WK(12)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./
C
C SPECIFIC HEAT DATA FROM GORDON AND MCBRDE (1971)
C
DATA A0/
1 .44608041E+01, .30981719E-02, -.12392571E-05, .22741325E-09,
1 -.15525954E-13, -.48961442E+05, -.98635982E+00,
2 .27167633E+01, .29451374E-02, -.80224374E-06, .10226682E-09,
2 -.48472145E-14, -.29905826E+05, .66305671E+01,
3 .28963194E+01, .15154866E-02, -.57235277E-06, .99807393E-10,
3 -.65223555E-14, -.90586184E+03, .61615148E+01,
4 .36219535E+01, .72518264E-03, -.19652228E-06, .36201558E-10,
4 -.28945627E-14, -.12019825E+04, .36150960E+01,
5 .29840696E+01, .14891390E-02, -.57899684E-06, .36201558E-10,
5 -.69353550E-14, -.14245228E+05, .63479156E+01,
6 .31001901E+01, .51119464E-03, .52644210E-07, -.34909973E-10,
6 .36945345E-14, -.87738042E+03, -.19629421E+01,
7 .25000000E+01, 0., 0., 0., 0., .25471627E+05, -.46011763E+00,
8 .55420596E+01, -.27550619E-04, -.31028033E-08, .45510674E-11,
8 -.43680515E-15, .29230803E+03, .49203080E+01,
9 .29106427E+01, .95931650E-03, -.19441702E-06, .13756646E-10,
9 .14224542E-15, .39353815E+04, .54423445E+01,
# .31890000E+01, .13382281E-02, -.52899318E-06, .95919332E-10,
# -.64847932E-14, .98283290E+04, .67458126E+01/
C
C EQUILIBRIUM CONSTANT DATA FROM OLIKARA AND BORMAN (1975)
C
DATA KP/
1 .432168E+00, -.112464E+05, .267269E+01, -.745744E-04,
1 .242484E-08,
2 .310805E+00, -.129540E+05, .321779E+01, -.738336E-04,
2 .344645E-08,
3 -.141784E+00, -.213308E+04, .853461E+00, .355015E-04,
3 -.310227E-08,
4 .150879E-01, -.470959E+04, .646096E+00, .272805E-05,
4 -.154444E-08,
5 -.752364E+00, .124210E+05, -.260284E+01, .259556E-03,
5 -.162687E-07,
6 -.415302E-02, .148627E+05, -.475746E+01, .124699E-03,
6 -.900227E-08/
C
C MOLECULAR WEIGHTS
C
DATA M/44.01, 18.02, 28.008, 32., 28.01, 2.018, 1.009, 16.,
& 17.009, 30.004/
C
C OTHER DATA
C
DATA MAXITS/100/, TOL/3.0E-05/, RU/8.31434/
C
C MAKE SURE SOLID CARBON WILL NOT FORM
C
IER = 2
EPS = .210/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF(PHI .GT. (.210/EPS/(0.5*ALPHA - 0.5*GAMMA)) ) RETURN
IER = 0
C
C DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C REQUIRED OR IF T<1000 IN WHICH CASE FARG WILL SUFFICE
C
C
SUM = 0
DO 10 I = 1,10
10 SUM = SUM + Y(I)
IF( T .GT. 1000. .AND. SUM .GT. 0.998) GO TO 20
CALL FARG(P,T,PHI,1.,H,U,V,S,Y0,CP,DLVLT,DLVLP)
DO 30 I = 1,10
30 Y(I) = 0
DO 40 I = 1,6
40 Y(I) = Y0(I)
IF( T .LE. 1000. ) RETURN
20 CONTINUE
C
C
C EVALUATE THE REQUISITE CONSTANTS
C
PATM = .9869233*P
DO 50 I = 1,6
50 K(I) = 10.**(KP(1,I)*ALOG(T/1000.) + KP(2,I)/T + KP(3,I) +
& KP(4,I)*T + KP(5,I)*T*T)
C(1) = K(1)/SQRT(PATM)
C(2) = K(2)/SQRT(PATM)
C(3) = K(3)
C(4) = K(4)
C(5) = K(5)*SQRT(PATM)
C(6) = K(6)*SQRT(PATM)
D(1) = BETA/ALPHA
D(2) = (GAMMA + .42/EPS/PHI)/ALPHA
D(3) = (DELTA + 1.58/EPS/PHI)/ALPHA
C
C
C SET UP ITERATION LOOP FOR NEWTON-RAPHSON ITERATION AND
C INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C AT LOW TEMPERATURES ON 24 BIT MACHINES (VAX). SET
C TRICK = 1.0 ON 48 BIT MACHINES (CDC)
C
TRICK = 10.
CHECK = ABS(PHI - 1.0) .LT. TOL
IF(CHECK) PHI = PHI*(1.0 + SIGN(TOL,PHI - 1.0))
ICHECK = 0
DO 5 J = 1,MAXITS
DO 4 JJ = 1,10
4 Y(JJ) = AMIN1(1.0,AMAX1(Y(JJ),1.0E-25))
D76 = 0.5*C(1)/SQRT(Y(6))
D84 = 0.5*C(1)/SQRT(Y(4))
D94 = 0.5*C(3)*SQRT(Y(6)/Y(4))
D96 = 0.5*C(3)*SQRT(Y(4)/Y(6))
D103 = 0.5*C(4)*SQRT(Y(4)/Y(3))
D104 = 0.5*C(4)*SQRT(Y(3)/Y(4))
D24 = 0.5*C(5)*Y(6)/SQRT(Y(4))
D26 = C(5)*SQRT(Y(4))
D14 = 0.5*C(6)*Y(5)/SQRT(Y(4))
D15 = C(6)*SQRT(Y(4))
A(1,1) = 1. + D103
A(1,2) = D14 + D24 + 1. +D84 + D104 + D96
A(1,3) = D15 + 1.0
A(1,4) = D26 + 1.0 + D76 +D96
A(2,1) = 0.0
A(2,2) = 2.0*D24 + D94 - D(1)*D14
A(2,3) = -D(1)*D15 - D(1)
A(2,4) = 2.0*D26 + 2.0 + D76 + D96
A(3,1) = D103
A(3,2) = 2.0*D14 + D24 + 2.0 + D84 + D04 + D104 - D(2)*D14
A(3,3) = 2.0*D15 + 1.0 - D(2)*D15 - D(2)
A(3,4) = D26 + D96
A(4,1) = 2.0 +D103
A(4,2) = D104 - D(3)*D14
A(4,3) = -D(3)*D15 - D(3)
A(4,4) = 0
C
C SOLVE THE MATRIX EQUATIONS 3.81 FOR COMPOSITION CORRECTIONS
C
SUM = 0
DO 55 I = 1,10
55 SUM = SUM + Y(I)
B(1) = -(SUM - 1.0)
B(2) = -(2.0*Y(2) + 2.0*Y(6) + Y(7) + Y(9) - D(1)*Y(1) -
& D(1)*Y(5))
B(3) = -(2.0*Y(1) + Y(2) + 2.0*Y(4) + Y(5) + Y(8) + Y(9)
& + Y(10) - D(2)*Y(1) - D(2)*Y(5))
B(4) = -(2.0*Y(3) + Y(10) - D(3)*Y(1) - D(3)*Y(5))
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 0, WK, IERROR)
ERROR = 0.
DO 56 L = 3,6
LL = L - 2
Y(L) = Y(L) + B(LL)/TRICK
ERROR = AMAX1(ERROR,ABS(B(LL)))
Y(L) = AMIN(1.0,AMAX(Y(L),1.0E-25))
56 CONTINUE
Y(7) = C(1)*SQRT(Y(6))
Y(8) = C(2)*SQRT(Y(4))
Y(9) = C(3)*SQRT(Y(4)*Y(6))
Y(10) = C(4)*SQRT(Y(4)*Y(3))
Y(2) = C(5)*SQRT(Y(4))*Y(6)
Y(1) = C(6)*SQRT(Y(4))*Y(5)
IF(ERROR .LT. TOL) ICHECK = ICHECK +1
IF(ERROR .LT. TOL .AND. ICHECK .GE. 2) GO TO 57
5 CONTINUE
IER = 1
C
C COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
C
57 DO 60 I = 1,6
60 DKDT(I) = 2.302585*K(I)*(KP(1,I)/T - KP(2,I)/T/T + KP(4,I)
& + 2.0*KP(5,I)*T)
DCDT(1) = DKDT(1)/SQRT(PATM)
DCDT(2) = DKDT(2)/SQRT(PATM)
DCDT(3) = DKDT(3)
DCDT(4) = DKDT(4)
DCDT(5) = DKDT(5)*SQRT(PATM)
DCDT(6) = DKDT(6)*SQRT(PATM)
DCDP(1) = -0.5*C(1)/P
DCDP(2) = -0.5*C(2)/P
DCDP(5) = 0.5*C(5)/P
DCDP(6) = 0.5*C(6)/P
C
X1 = Y(1)/C(6)
X2 = Y(2)/C(5)
X7 = Y(7)/C(1)
X8 = Y(8)/C(2)
X9 = Y(9)/C(3)
X10 = Y(10)/C(4)
DFDT(1) = DCDT(6)*X1 + DCDT(5)*X2 + DCDT(1)*X7 +DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10
DFDT(2) = 2.0*DCDT(5)*X2 + DCDT(1)*X7 + DCDT(3)*X9 -
& D(1)*DCDT(6)*X1
DFDT(3) = 2.0*DCDT(6)*X1 + DCDT(5)*X2 + DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10 - D(2)*DCDT(6)*X1
DFDT(4) = DCDT(4)*X10 - D(3)*DCDT(6)*X1
DFDP(1) = DCDP(6)*X1 + DCDP(5)*X2 + DCDP(1)*X7 + DCDP(2)*X8
DFDP(2) = 2.0*DCDP(5)*X2 + DCDP(1)*X7 - D(1)*DCDP(6)*X1
DFDP(3) = 2.0*DCDP(6)*X1 + DCDP(5)*X2 + DCDP(2)*X8 -
& D(2)*DCDP(6)*X1
DFDP(4) = -D(3)*DCDP(6)*X1
C
C SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C AND THEN DETERMINE THE DEPENDENT DERIVATICES
C
DO 70 I = 1,4
70 B(I) = -DFDT(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDT(3) = B(1)
DYDT(4) = B(2)
DYDT(5) = B(3)
DYDT(6) = B(4)
DYDT(1) = SQRT(Y(4))*Y(5)*DCDT(6) + D14*DYDT(4) + D15*DYDT(5)
DYDT(2) = SQRT(Y(4))*Y(6)*DCDT(5) + D24*DYDT(4) + D26*DYDT(6)
DYDT(7) = SQRT(Y(6))*DCDT(1) + D76*DYDT(6)
DYDT(8) = SQRT(Y(4))*DCDT(2) + D84*DYDT(4)
DYDT(9) = SQRT(Y(4)*Y(6))*DCDT(3) + D94*DYDT(4) + D96*DYDT(6)
DYDT(10) = SQRT(Y(4)*Y(3))*DCDT(4) + D104*DYDT(4) + D103*DYDT(3)
C
DO 80 I = 1,4
80 B(I) = -DFDP(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDP(3) = B(1)
DYDP(4) = B(2)
DYDP(5) = B(3)
DYDP(6) = B(4)
DYDP(1) = SQRT(Y(4))*Y(5)*DCDP(6) + D14*DYDP(4) + D15*DYDP(5)
DYDP(2) = SQRT(Y(4))*Y(6)*DCDP(5) + D24*DYDP(4) + D26*DYDP(6)
DYDP(7) = SQRT(Y(6))*DCDP(1) + D76*DYDP(6)
DYDP(8) = SQRT(Y(4))*DCDP(2) + D84*DYDP(4)
DYDP(9) = D94*DYDP(4) + D96*DYDP(6)
DYDP(10) = D104*DYDP(4) + D103*DYDP(3)
C
DO 90 I = 1,10
CP0(I) = A0(1,I) + A0(2,I)*T + A0(3,I)*T*T + A0(4,I)*T**3
& + A0(5,I)*T**4
H0(I) = A0(1,I) + A0(2,I)/2.*T + A0(3,I)/3.*T*T + A0(4,I)/4.*
& T**3 + A0(5,1)/5.*T**4 + A0(6,I)/T
90 S0(I) = A0(1,I)*ALOG(T) + A0(2,I)*T + A0(3,I)/2.*T**2
& + A0(4,I)/3.*T**3 + A0(5,I)/4.*T**4 +A0(7,I)
C
C Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
C
Y(1) = (2.*Y(3) + Y(10))/D(3) - Y(5)
Y(2) = (D(1)/D(3)*(2.*Y(3) + Y(10)) - 2.*Y(6) - Y(7) - Y(9))/2
C
MW = 0.
CP = 0.
H = 0.
S = 0.
MT = 0.
MP = 0.
DO 100 I = 1,10
IF( Y(I) .LE. 1.0E-25 ) GO TO 100
H = H + H0(I)*Y(I)
MW = MW + M(I)*Y(I)
MT = MT + M(I)*DYDT(I)
MP = MP + M(I)*DYDP(I)
CP = CP + Y(I)*CP0(I) + H0(I)*T*DYDT(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
100 CONTINUE
C
R = RU/MW
V = 10.*R*T/P
CP = R*(CP - H*T*MT/MW)
DLVLT = 1.0 + AMAX1(-T*MT/MW,0.)
DLVLP = -1.0 - AMAX1(P*MP/MW,0.)
H = H*R*T
S = R*(-ALOG(PATM) + 5)
U = H - R*T
RETURN
END
C
C
SUBROUTINE FARG (P, T, PHI, F, H, U, V, S, Y, CP, DLVLT, DLVLP)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF A FUEL, AIR, RESIDUAL
C GAS MIXTURE
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C F - RESIDUAL MASS FRACTION
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/g)
C S - ENTROPY (J/G/K)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2, 2=H2O, 3=N2, 4O2, 5=CO, 6=H2
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C
C REMARKS:
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE
C 3. ENTHALPIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. THE FUEL MOLE FRACTION IS UNITY MINUS THE SUM OF Y(I)
C
REAL M, MW, MRES, MFA, MFUEL, K, NU, N2
LOGICAL RICH, LEAN
DIMENSION A(7,6), TABLE(6), M(6), NU(6), Y(6), CP0(6), H0(6),
& S0(6)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./,
& A0/4.0652/, B0/6.0977E-02/, C0/-1.8801E-05/,
& D0/-3.5880E+04/, E0/15.45/
C
C TABLE 3.1, I = 1,6
C
DATA A/
1 .24007797E+01, .87350957E-02, -.6607878E-05, .20021861E-08,
1 .63274039E-15, -.48377527E+05, .96951457E+01,
2 .40701275E+01, -.11084499E-02, .41521180E-05, -.29637404E-08,
2 .80702103E-12, -.30279722E+05, -.32270046E+00,
3 .36748261E+01, -.12081500E-02, .23240102E-05, -.63217559E-09,
3 -.22577253E-12, -.10611588E+04, .23580424E+01,
4 .36255985E+01, -.18782184E-02, .70554544E-05, -.67635137E-08,
4 .21555993E-11, -.10475226E+04, .43052778E+01,
5 .37100928E+01, -.16190964E-02, .36923594E-05, -.20319674E-08,
5 .23953344E-12, -.14356310E+05, .29555355E+01,
6 .30574451E+01, .26765200E-02, -.58099162E-05, .55210391E-08,
6 -.18122739E-11, -.98890474E+03, -.22997056E+01/
C
C OTHER DATA
C
DATA RU/8.315/, TABLE/-1.,1.,0.,0.,1.,-1./, M/44.01, 18.02,
& 28.008, 32.000, 28.01, 2.018/
C
C COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
C
RICH = PHI .GT. 1.0
LEAN = .NOT. RICH
DLVLT = 1.0
DLVLP = -1.0
EPS = .21/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF (RICH) GO TO 10
NU(1) = ALPHA*PHI*EPS
NU(2) = BETA*PHI*EPS/2
NU(3) = 0.79 + DELTA*PHI*EPS/2
NU(4) = 0.21*(1.0 - PHI)
NU(5) = 0.
NU(6) = 0.
DCDT = 0.
GO TO 20
10 Z = 1000./T
K = EXP(2.743 + Z*(-1.761 + Z*(-1.611 + Z*.2803)))
DKDT = -K*(-1.761 + Z*(-3.222 + Z*.8409))/1000.
A1 = 1.0 - K
B = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + K*(.42*(PHI - 1.) +
& ALPHA*PHI*EPS)
C = -.42*ALPHA*PHI*EPS*(PHI - 1.0)*K
NU(5) = (-B + SQRT(B*B - 4.*A1*C))/2./A1
DCDT = DKDT*(NU(5)**2 - NU(5)*(.42*(PHI - 1.) + ALPHA*
& PHI*EPS) + .42*ALPHA*PHI*EPS*(PHI - 1.))/
& (2.*NU(5)*A1 + B)
NU(1) = ALPHA*PHI*EPS - NU(5)
NU(2) = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + NU(5)
NU(3) = .79 + DELTA*PHI*EPS/2.
NU(4) = 0.
NU(6) = .42*(PHI - 1.) - NU(5)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
C
20 TMOLES = 0.
DO 30 I = 1,6
30 TMOLES = TMOLES + NU(I)
MRES = 0.
DO 40 I = 1,6
Y(I) = NU(I)/TMOLES
40 MRES = MRES + Y(I)*M(I)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF FURLE-AIR
C
FUEL = EPS*PHI/(1. + EPS*PHI)
O2 = .21/(1. + EPS*PHI)
N2 = .79/(1. + EPS*PHI)
MFA = FUEL*(12.01*ALPHA + 1.008*BETA + 16.*GAMMA + 14.01*
& DELTA) + 32.*O2 + 28.02*N2
C
C COMPUTE MOLE FRACTIONS OF FUEL-AIR-RESIDUAL GAS
C
YRES = F/(F + MRES/MFA*(1. - F))
DO 50 I = 1,6
50 Y(I) = Y(I)*YRES
YFUEL = FUEL*(1. - YRES)
Y(3) = Y(3) + N2*(1. - YRES)
Y(4) = Y(4) + O2*(1. - YRES)
C
C COMPUTE COMPONENT PROPERTIES
C
DO 60 I = 1,6
CP0(I) = A(1,I) + A(2,I)*T + A(3,I)*T**2 + A(4,I)*T**3 +
& A(5,I)*T**4
H0(I) = A(1,I) + A(2,I)/2.*T + A(3,I)/3.*T**2 + A(4,I)/4.*
& T**3 + A(5,I)/5.*T**4 + A(6,I)/T
60 S0(I) = A(1,I)*ALOG(T) + A(2,I)*T + A(3,I)/2.*T**2 +
& A(4,I)/3.*T**3 + A(5,I)/4.*T**4 + A(7,1)
C
MFUEL = 12.01*ALPHA + 1.008*BETA + 16.000*GAMMA + 14.01*DELTA
CPFUEL = A0 + B0*T + C0*T**2
HFUEL = A0 + B0/2.*T + C0/3.*T**2 + D0/T
S0FUEL = A0*ALOG(T) + B0*T + C0/2.*T**2 + E0
C
C COMPUTE PROPERTIES OF MIXTURE
C
H = HFUEL*YFUEL
S = (S0FUEL - ALOG(YFUEL))*YFUEL
CP = CPFUEL*YFUEL
MW = MFUEL*YFUEL
DO 70 I = 1,6
H = H + H0(I)*Y(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
CP = CP + CP0(I)*Y(I) + H0(I)*T*TABLE(I)*DCDT*YRES/TMOLES
70 MW = MW + Y(I)*M(I)
C
R = RU/MW
H = R*T*H
U = H - R*T
V = 10.*R*T/P
S = R*(-ALOG(.9869*P) + S)
CP = R*CP
C
RETURN
END
C
C
SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W)
INTEGER N, IND, NW, K
DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
C Added external statement for fcn to avoid a warning message.
external fcn
C
C***********************************************************************
C *
C PURPOSE - THIS IS A RUNGE-KUTTA SUBROUTINE BASED ON VERNER'S *
C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO *
C THE SOLUTION OF A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL *
C EQUATIONS WITH INITIAL CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL *
C ERROR PROPORTIONAL TO A TOLERANCE SPECIFIED BY THE USER. (THE *
C PROPORTIONALITY DEPENDS ON THE KIND OF ERROR CONTROL THAT IS USED, *
C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) *
C *
C VARIOUS OPTIONS ARE AVAILABLE TO THE USER, INCLUDING DIFFERENT *
C KINDS OF ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS *
C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE CALCULATION (AND *
C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES. *
C *
C THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS. HOWEVER, A GOOD *
C VARIABLE-ORDER-ADAMS METHOD WILL PROBABLY BE MORE EFFICIENT IF THE *
C FUNCTION EVALUATIONS ARE VERY COSTLY. SUCH A METHOD WOULD ALSO BE *
C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE *
C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE FOR EXAMPLE *
C WITH GRAPHICAL OUTPUT. *
C *
C HULL-ENRIGHT-JACKSON 1/10/76
I am getting a run-time error when trying to compile my program and as I
am not very experienced in Fortran I am not sure what the error is or how
to fix it. The program has been lifted out of a book called "Interneal
Combustion Engines:Applied Thermosciences" and is a program to model the
heat release rate of a an engine. It requires a couplke of subroutines,
DVERK and LEQIF which are IMSL routines and I do not have. I have found
the code for DVERK but unfortunately have not found the code for LEQIF
which is why I assume I am receiving the run-time error. the error is a
zero or negative argument to a logarithm routine.
FARG - in file release2_good.for at line 664 [+1738]
main - in file release2_good.for at line 41 [+0190]
The main code is listed below and any help in resolving this matter would be greatly appreciated. Or if anybody else has a heat release program that they are willing to share that would be great too.
I can e-mail the code if the formatting is not retained in this thread to make things easier.
Thanks
Chris
C GIVEN IN BLOCK DATA SUBPROGRAM
C COMPRESSION RATIO - R
C BORE - B (CM)
C STROKE - S (CM)
C HALF-STROKE TO ROD RATIO - EPS
C ENGINE SPEED - RPM
C HEAT TRANSFER COEFFICIENT - H (J/S/M**2/K)
C BLOWBY COEFFICIENT - C (L/S)
C EQUIVALENCE RATIO - PHI
C RESIDUAL FRACTION - F
C INITIAL PRESSURE - P1 (BAR)
C INITIAL TEMPERATURE - T1 (K)
C WALL TEMPERATURE - TW (K)
C BURN ANGLE - THETAB (DEG)
C START OF HEAT RELEASE - THETAS (DEG ATDC)
C
C FIND:
C INDICATED THERMAL EFFICIENCY - ETA
C INDICATED MEAN EFFECTIVE PRESSURE - IMEP (BAR)
C
C NOTES:
C 1. COSINE BURNING LAW IS EMPLOYED, EQN. 4.61
C ` 2. FUEL - GASOLINE (C7H17)
C
EXTERNAL RATES
REAL IMEP, MNOT, M
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT, UFINAL
DIMENSION Y1(6), Y2(10), Y(6), COM(24), WM(6,9)
DATA AO/47870./, FS/.06548/, Y/1., 10000., 350., 3*0.0/,
& N/6/, THETA/-180./, THETAE/-179./, TOL/.001/, IND/1/,
& NW/6/, PI/3.141593/, Y2/10*0.0/
C
C LABEL HEADER ON OUTPUT PAGE AND PRINT INITIAL CONDITIONS
C
OMEGA = RPM*PI/30
CALL AUXLRY(THETA, V, X, EM)
Y(1) = P1
Y(3) = T1
CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, SU, Y1, CP, DLVLT,
& DLVLP)
MNOT = V/VU
UNOT = MNOT*U
WRITE(6,10)
10 FORMAT(5X,'THETA',4X,'VOLUME',3X,'BURN FRAC',2X,'PRESSURE',3X,
& 'BURN TEMP',2X,'T UNBURN',3X,'WORK',3X,'HEAT LOSS',5X,
& 'MASS',3X,'H-LEAK')
WRITE(6,15)
15 FORMAT(4X,'DEG ATDC',2X,'CM**3',6X,'--',8X,'BARS',7X,'KELVIN',
& 5X,'KELVIN',3X,'JOULES',4X,'JOULES',6X,'GRAMS',2X,
& 'JOULES')
WRITE(6,20) THETA, V, X, (Y(I),I=1,5), MNOT, Y(6)
20 FORMAT(4X,F5.0,6X,F5.0,4X,F5.3,6X,F5.2,7X,F5.0,5X,F5.0,
& 5X,F5.0,3X,F5.0,7X,F5.3,3X,F5.2)
C
C INTEGRATE THE RATE EQUATIONS 4.80 - 4.85 AND PRINT ANSWERS
C EVERY 10 DEGREES
C
DO 30 I = 1,36
DO 25 J = 1,10
CALL DVERK (N, RATES, THETA, Y, THETAE, TOL, IND, COM, NW, WM)
CALL AUXLRY(THETAE, V, X, EM)
M = EM*MNOT
THETA = THETAE
THETAE = THETA + 1.
IF(THETAS .GE. THETA .AND. THETAS .LT. THETAE)
& CALL TINITL(Y(1), T(3), PHI, F, Y(2))
IF(THETA .GT. THETAS + THETAB) Y(3) = 10000.
25 CONTINUE
WRITE(6,20) THETA, V, X, (Y(J),J=1,5), M, Y(6)
30 CONTINUE
C
C COMPUTE AND WRITE THE EAT, IMEP AND ERRORS
C
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, SB, Y2, CP, DLVLT, DLVLP,
& IER)
UFINAL = U*M
ERROR1 = 1.0 - VB*M/V
ERROR2 = 1.0 + Y(4)/(UFINAL - UNOT + Y(5) + Y(6))
ETA = Y(4)/MNOT*(1. + PHI*FS*(1. - F))/PHI/FS/(1. - F)/AO
IMEP = Y(4)/(PI/4.*B**2*S)*10.
WRITE(6,40) ETA, IMEP, ERROR1, ERROR2
40 FORMAT('ETA = ',1PE11.4,' IMEP = ',1PE11.4,'ERROR1 = ',1PE11.4,
& ' ERROR2 = ',1PE11.4)
STOP
END
C
SUBROUTINE RATES(N, THETA, Y, YPRIME)
DIMENSION Y1(6), Y2(10)
REAL Y(N), YPRIME(N), M, MNOT
COMMON /ENGINE/ R, BORE, STROKE, EPS, RPM, HEAT, CEE, THETAB,
& THETAS, PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT,
& UFINAL
DATA PI/3.141593/, Y2/10*0.0/
CALL AUXLRY(THETA, V, X, EM)
M = EM*MNOT
RAD = THETA*PI/180.
DUMB = SQRT(1. - (EPS*SIN(RAD))**2)
DV = PI/8.*BORE**2*STROKE*SIN(RAD)*(1. + EPS*COS(RAD)/DUMB)
A = (DV + V*CEE/OMEGA)/M
C1 = HEAT*(PI*BORE**2/2. + 4.*V/BORE)/OMEGA/M/10000.
C0 = SQRT(X)
C
C THREE DIFFERENT COMPUTATIONS ARE REQUIRED DEPENDING
C UPON THE SIZE OF THE MASS FRACTION BURNED
C
IF( X .GT. .999 ) GO TO 20
IF( X .GT. .001 ) GO TO 10
C
C COMPRESSION
C
CALL FARG(Y(1), Y(3), PHI, F, HL, U, VU, S, Y1, CP, DLVLT,
& DLVLP)
B = C1*VU/CP*DLVLT*(1.0 - TW/Y(3))
C = 0.
D = 0.
E = VU**2/CP/Y(3)*DLVLT**2 + 10.*VU/Y(1)*DLVLP
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = 0.
YPRIME(3) = -C1/CP*(Y(3) - TW) + VU/CP*DLVLT*YPRIME(1)/10.
GO TO 30
C
C COMBUSTION
C
10 CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, S, Y1, CPU, DLVLTU,
& DLVLP)
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, S, Y2, CPB, DDLVLTB,
& DLVLPB, IER)
B = C1*(VB/CPB*DLVLTB*C0*(1. - TW/Y(2)) + VU/CPU*DLVLTU*
& (1. - C0)*(1.- TW/Y(3)))
DX = 0.5*SIN(PI*(THETA - THETAS)/THETAB)*180./THETAB
C = -(VB - VU)*DX - VB/Y(2)*DLVLTB*(HU - HB)/CPB*(DX - (X -
& X**2)*CEE/OMEGA)
D = X*(VB**2/CPB/Y(2)*DLVLTB**2 + 10.*VB/Y(1)*DLVLPB)
E = (1. - X)*(VU**2/CPU/Y(3)*DLVLTU**2 + 10.*VU/Y(1)*DLVLPU)
HL = (1.0 - X**2)*HU + X**2*HB
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CPB/C0*(Y(2) - TW) + VB/CPB*DLVLTB*
& YPRIME(1)/10. + (HU - HB)/CPB*(DX/X - (1. - X)*
& CEE/OMEGA)
YPRIME(3) = -C1/CPU/(1. + CO)*(Y(3) - TW) + VU/CPU*DLVLTU*
& YPRIME(1)/10.
GO TO 30
C
C EXPANSION
C
20 CALL ECP(Y(1), Y(2), PHI, HL, U, VB, S, Y2, CP, DLVLT, DLVLP,
& IER)
B = C1*VB/CP*DLVLT*(1. - TW/Y(2))
C = 0.
D = VB**2/CP/Y(2)*DLVLT**2 + 10.*VB/Y(1)*DLVLP
E = 0.
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CP*(Y(2) - TW) + VB/CP*DLVLT*YPRIME(1)/10.
YPRIME(3) = 0.
C
C ALL CASES
C
30 YPRIME(4) = Y(1)*DV/10.
YPRIME(5) = C1*M*(C0*(Y(2) - TW) + (1. - C0)*(Y(3) - TW))
YPRIME(6) = CEE*M/OMEGA*HL
DO 40 I = 1,N
40 YPRIME(I) = YPRIME(I)*PI/180.
RETURN
END
C
SUBROUTINE AUXLRY(THETA, V, X, EM)
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA PI/3.141593/
RAD = THETA*PI/180.
VTDC = PI/4.*B**2*S/(R - 1.)
V = VTDC*(1. + (R - 1.)/2.*(1. - COS(RAD) + 1./EPS*(1. -
& SQRT(1. - (EPS*SIN(RAD))**2))))
X = 0.5*(1. - COS(PI*(THETA - THETAS)/THETAB))
IF(THETA .LE. THETAS) X = 0.
IF(THETA .GE. THETAS + THETAB) X = 1.
EM = EXP(-C*(RAD + PI)/OMEGA)
RETURN
END
C
BLOCK DATA
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA R/10./, B/10./, S/8./, EPS/.25/, RPM/2000./, H/500./,
& C/0.8/, THETAB/60./, THETAS/-35./, PHI/0.8/, F/.1/,
& P1/1.0/, T1/350./, TW/420./
C
END
C
SUBROUTINE TINITL(P, TU, PHI, F, TB)
DIMENSION Y1(6), Y2(10)
DATA MAXITS/50/, TOL/.001/, Y2/10*0.0/
TB = 2000.
CALL FARG(P, TU, PHI, F, HU, U, V, S, Y1, CP, DLVLT, DLVLP)
DO 10 I = 1,MAXITS
CALL ECP(P, TB, PHI, HB, U, V, S, Y2, CP, DLVLT, DLVLP, IER)
DELT = (HU - HB)/CP
TB = TB + DELT
IF(ABS(DELT/TB) .LT. TOL) RETURN
10 CONTINUE
WRITE(6,20)
20 FORMAT('CONVERGENCE FAILURE IN TINITL')
RETURN
END
SUBROUTINE ECP (P, T, PHI, H, U, V, S, Y, CP, DLVLT, DLVLP, IER)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF EQUILIBRIUM COMBUSTION
C PRODUCTS
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE )K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/G)
C S - ENTROPY (J/G/K)
C Y - A TEN DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2. 2=H2O, 3=N2, 4=O2, 5=CO, 6=H2
C 7=H, 8=O, 9=OH, 10=NO
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C IER - AN ERROR MESSAGE
C 0 - NORMAL
C 1 - CONVERGENCE FAILURE: COMPOSITION LOOP
C 2 - CALLED WITH TOO HIGH AN EQUIVALENCE RATIO; C(S)
C AND OTHER SPECIES WOULD FORM
C
REAL M, MW, K, KP, MT, MP
LOGICAL CHECK
DIMENSION M(10), K(6), C(6), Y(10), D(3), B(4), A(4,4), Y0(6),
& KP(5,6), DFDT(4), DCDT(6), DKDT(6), DCDP(6),
& DYDT(10), DYDP(10), A0(7,10), DFDP(4), CP0(10),
& H0(10), S0(10), WK(12)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./
C
C SPECIFIC HEAT DATA FROM GORDON AND MCBRDE (1971)
C
DATA A0/
1 .44608041E+01, .30981719E-02, -.12392571E-05, .22741325E-09,
1 -.15525954E-13, -.48961442E+05, -.98635982E+00,
2 .27167633E+01, .29451374E-02, -.80224374E-06, .10226682E-09,
2 -.48472145E-14, -.29905826E+05, .66305671E+01,
3 .28963194E+01, .15154866E-02, -.57235277E-06, .99807393E-10,
3 -.65223555E-14, -.90586184E+03, .61615148E+01,
4 .36219535E+01, .72518264E-03, -.19652228E-06, .36201558E-10,
4 -.28945627E-14, -.12019825E+04, .36150960E+01,
5 .29840696E+01, .14891390E-02, -.57899684E-06, .36201558E-10,
5 -.69353550E-14, -.14245228E+05, .63479156E+01,
6 .31001901E+01, .51119464E-03, .52644210E-07, -.34909973E-10,
6 .36945345E-14, -.87738042E+03, -.19629421E+01,
7 .25000000E+01, 0., 0., 0., 0., .25471627E+05, -.46011763E+00,
8 .55420596E+01, -.27550619E-04, -.31028033E-08, .45510674E-11,
8 -.43680515E-15, .29230803E+03, .49203080E+01,
9 .29106427E+01, .95931650E-03, -.19441702E-06, .13756646E-10,
9 .14224542E-15, .39353815E+04, .54423445E+01,
# .31890000E+01, .13382281E-02, -.52899318E-06, .95919332E-10,
# -.64847932E-14, .98283290E+04, .67458126E+01/
C
C EQUILIBRIUM CONSTANT DATA FROM OLIKARA AND BORMAN (1975)
C
DATA KP/
1 .432168E+00, -.112464E+05, .267269E+01, -.745744E-04,
1 .242484E-08,
2 .310805E+00, -.129540E+05, .321779E+01, -.738336E-04,
2 .344645E-08,
3 -.141784E+00, -.213308E+04, .853461E+00, .355015E-04,
3 -.310227E-08,
4 .150879E-01, -.470959E+04, .646096E+00, .272805E-05,
4 -.154444E-08,
5 -.752364E+00, .124210E+05, -.260284E+01, .259556E-03,
5 -.162687E-07,
6 -.415302E-02, .148627E+05, -.475746E+01, .124699E-03,
6 -.900227E-08/
C
C MOLECULAR WEIGHTS
C
DATA M/44.01, 18.02, 28.008, 32., 28.01, 2.018, 1.009, 16.,
& 17.009, 30.004/
C
C OTHER DATA
C
DATA MAXITS/100/, TOL/3.0E-05/, RU/8.31434/
C
C MAKE SURE SOLID CARBON WILL NOT FORM
C
IER = 2
EPS = .210/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF(PHI .GT. (.210/EPS/(0.5*ALPHA - 0.5*GAMMA)) ) RETURN
IER = 0
C
C DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C REQUIRED OR IF T<1000 IN WHICH CASE FARG WILL SUFFICE
C
C
SUM = 0
DO 10 I = 1,10
10 SUM = SUM + Y(I)
IF( T .GT. 1000. .AND. SUM .GT. 0.998) GO TO 20
CALL FARG(P,T,PHI,1.,H,U,V,S,Y0,CP,DLVLT,DLVLP)
DO 30 I = 1,10
30 Y(I) = 0
DO 40 I = 1,6
40 Y(I) = Y0(I)
IF( T .LE. 1000. ) RETURN
20 CONTINUE
C
C
C EVALUATE THE REQUISITE CONSTANTS
C
PATM = .9869233*P
DO 50 I = 1,6
50 K(I) = 10.**(KP(1,I)*ALOG(T/1000.) + KP(2,I)/T + KP(3,I) +
& KP(4,I)*T + KP(5,I)*T*T)
C(1) = K(1)/SQRT(PATM)
C(2) = K(2)/SQRT(PATM)
C(3) = K(3)
C(4) = K(4)
C(5) = K(5)*SQRT(PATM)
C(6) = K(6)*SQRT(PATM)
D(1) = BETA/ALPHA
D(2) = (GAMMA + .42/EPS/PHI)/ALPHA
D(3) = (DELTA + 1.58/EPS/PHI)/ALPHA
C
C
C SET UP ITERATION LOOP FOR NEWTON-RAPHSON ITERATION AND
C INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C AT LOW TEMPERATURES ON 24 BIT MACHINES (VAX). SET
C TRICK = 1.0 ON 48 BIT MACHINES (CDC)
C
TRICK = 10.
CHECK = ABS(PHI - 1.0) .LT. TOL
IF(CHECK) PHI = PHI*(1.0 + SIGN(TOL,PHI - 1.0))
ICHECK = 0
DO 5 J = 1,MAXITS
DO 4 JJ = 1,10
4 Y(JJ) = AMIN1(1.0,AMAX1(Y(JJ),1.0E-25))
D76 = 0.5*C(1)/SQRT(Y(6))
D84 = 0.5*C(1)/SQRT(Y(4))
D94 = 0.5*C(3)*SQRT(Y(6)/Y(4))
D96 = 0.5*C(3)*SQRT(Y(4)/Y(6))
D103 = 0.5*C(4)*SQRT(Y(4)/Y(3))
D104 = 0.5*C(4)*SQRT(Y(3)/Y(4))
D24 = 0.5*C(5)*Y(6)/SQRT(Y(4))
D26 = C(5)*SQRT(Y(4))
D14 = 0.5*C(6)*Y(5)/SQRT(Y(4))
D15 = C(6)*SQRT(Y(4))
A(1,1) = 1. + D103
A(1,2) = D14 + D24 + 1. +D84 + D104 + D96
A(1,3) = D15 + 1.0
A(1,4) = D26 + 1.0 + D76 +D96
A(2,1) = 0.0
A(2,2) = 2.0*D24 + D94 - D(1)*D14
A(2,3) = -D(1)*D15 - D(1)
A(2,4) = 2.0*D26 + 2.0 + D76 + D96
A(3,1) = D103
A(3,2) = 2.0*D14 + D24 + 2.0 + D84 + D04 + D104 - D(2)*D14
A(3,3) = 2.0*D15 + 1.0 - D(2)*D15 - D(2)
A(3,4) = D26 + D96
A(4,1) = 2.0 +D103
A(4,2) = D104 - D(3)*D14
A(4,3) = -D(3)*D15 - D(3)
A(4,4) = 0
C
C SOLVE THE MATRIX EQUATIONS 3.81 FOR COMPOSITION CORRECTIONS
C
SUM = 0
DO 55 I = 1,10
55 SUM = SUM + Y(I)
B(1) = -(SUM - 1.0)
B(2) = -(2.0*Y(2) + 2.0*Y(6) + Y(7) + Y(9) - D(1)*Y(1) -
& D(1)*Y(5))
B(3) = -(2.0*Y(1) + Y(2) + 2.0*Y(4) + Y(5) + Y(8) + Y(9)
& + Y(10) - D(2)*Y(1) - D(2)*Y(5))
B(4) = -(2.0*Y(3) + Y(10) - D(3)*Y(1) - D(3)*Y(5))
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 0, WK, IERROR)
ERROR = 0.
DO 56 L = 3,6
LL = L - 2
Y(L) = Y(L) + B(LL)/TRICK
ERROR = AMAX1(ERROR,ABS(B(LL)))
Y(L) = AMIN(1.0,AMAX(Y(L),1.0E-25))
56 CONTINUE
Y(7) = C(1)*SQRT(Y(6))
Y(8) = C(2)*SQRT(Y(4))
Y(9) = C(3)*SQRT(Y(4)*Y(6))
Y(10) = C(4)*SQRT(Y(4)*Y(3))
Y(2) = C(5)*SQRT(Y(4))*Y(6)
Y(1) = C(6)*SQRT(Y(4))*Y(5)
IF(ERROR .LT. TOL) ICHECK = ICHECK +1
IF(ERROR .LT. TOL .AND. ICHECK .GE. 2) GO TO 57
5 CONTINUE
IER = 1
C
C COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
C
57 DO 60 I = 1,6
60 DKDT(I) = 2.302585*K(I)*(KP(1,I)/T - KP(2,I)/T/T + KP(4,I)
& + 2.0*KP(5,I)*T)
DCDT(1) = DKDT(1)/SQRT(PATM)
DCDT(2) = DKDT(2)/SQRT(PATM)
DCDT(3) = DKDT(3)
DCDT(4) = DKDT(4)
DCDT(5) = DKDT(5)*SQRT(PATM)
DCDT(6) = DKDT(6)*SQRT(PATM)
DCDP(1) = -0.5*C(1)/P
DCDP(2) = -0.5*C(2)/P
DCDP(5) = 0.5*C(5)/P
DCDP(6) = 0.5*C(6)/P
C
X1 = Y(1)/C(6)
X2 = Y(2)/C(5)
X7 = Y(7)/C(1)
X8 = Y(8)/C(2)
X9 = Y(9)/C(3)
X10 = Y(10)/C(4)
DFDT(1) = DCDT(6)*X1 + DCDT(5)*X2 + DCDT(1)*X7 +DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10
DFDT(2) = 2.0*DCDT(5)*X2 + DCDT(1)*X7 + DCDT(3)*X9 -
& D(1)*DCDT(6)*X1
DFDT(3) = 2.0*DCDT(6)*X1 + DCDT(5)*X2 + DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10 - D(2)*DCDT(6)*X1
DFDT(4) = DCDT(4)*X10 - D(3)*DCDT(6)*X1
DFDP(1) = DCDP(6)*X1 + DCDP(5)*X2 + DCDP(1)*X7 + DCDP(2)*X8
DFDP(2) = 2.0*DCDP(5)*X2 + DCDP(1)*X7 - D(1)*DCDP(6)*X1
DFDP(3) = 2.0*DCDP(6)*X1 + DCDP(5)*X2 + DCDP(2)*X8 -
& D(2)*DCDP(6)*X1
DFDP(4) = -D(3)*DCDP(6)*X1
C
C SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C AND THEN DETERMINE THE DEPENDENT DERIVATICES
C
DO 70 I = 1,4
70 B(I) = -DFDT(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDT(3) = B(1)
DYDT(4) = B(2)
DYDT(5) = B(3)
DYDT(6) = B(4)
DYDT(1) = SQRT(Y(4))*Y(5)*DCDT(6) + D14*DYDT(4) + D15*DYDT(5)
DYDT(2) = SQRT(Y(4))*Y(6)*DCDT(5) + D24*DYDT(4) + D26*DYDT(6)
DYDT(7) = SQRT(Y(6))*DCDT(1) + D76*DYDT(6)
DYDT(8) = SQRT(Y(4))*DCDT(2) + D84*DYDT(4)
DYDT(9) = SQRT(Y(4)*Y(6))*DCDT(3) + D94*DYDT(4) + D96*DYDT(6)
DYDT(10) = SQRT(Y(4)*Y(3))*DCDT(4) + D104*DYDT(4) + D103*DYDT(3)
C
DO 80 I = 1,4
80 B(I) = -DFDP(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDP(3) = B(1)
DYDP(4) = B(2)
DYDP(5) = B(3)
DYDP(6) = B(4)
DYDP(1) = SQRT(Y(4))*Y(5)*DCDP(6) + D14*DYDP(4) + D15*DYDP(5)
DYDP(2) = SQRT(Y(4))*Y(6)*DCDP(5) + D24*DYDP(4) + D26*DYDP(6)
DYDP(7) = SQRT(Y(6))*DCDP(1) + D76*DYDP(6)
DYDP(8) = SQRT(Y(4))*DCDP(2) + D84*DYDP(4)
DYDP(9) = D94*DYDP(4) + D96*DYDP(6)
DYDP(10) = D104*DYDP(4) + D103*DYDP(3)
C
DO 90 I = 1,10
CP0(I) = A0(1,I) + A0(2,I)*T + A0(3,I)*T*T + A0(4,I)*T**3
& + A0(5,I)*T**4
H0(I) = A0(1,I) + A0(2,I)/2.*T + A0(3,I)/3.*T*T + A0(4,I)/4.*
& T**3 + A0(5,1)/5.*T**4 + A0(6,I)/T
90 S0(I) = A0(1,I)*ALOG(T) + A0(2,I)*T + A0(3,I)/2.*T**2
& + A0(4,I)/3.*T**3 + A0(5,I)/4.*T**4 +A0(7,I)
C
C Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
C
Y(1) = (2.*Y(3) + Y(10))/D(3) - Y(5)
Y(2) = (D(1)/D(3)*(2.*Y(3) + Y(10)) - 2.*Y(6) - Y(7) - Y(9))/2
C
MW = 0.
CP = 0.
H = 0.
S = 0.
MT = 0.
MP = 0.
DO 100 I = 1,10
IF( Y(I) .LE. 1.0E-25 ) GO TO 100
H = H + H0(I)*Y(I)
MW = MW + M(I)*Y(I)
MT = MT + M(I)*DYDT(I)
MP = MP + M(I)*DYDP(I)
CP = CP + Y(I)*CP0(I) + H0(I)*T*DYDT(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
100 CONTINUE
C
R = RU/MW
V = 10.*R*T/P
CP = R*(CP - H*T*MT/MW)
DLVLT = 1.0 + AMAX1(-T*MT/MW,0.)
DLVLP = -1.0 - AMAX1(P*MP/MW,0.)
H = H*R*T
S = R*(-ALOG(PATM) + 5)
U = H - R*T
RETURN
END
C
C
SUBROUTINE FARG (P, T, PHI, F, H, U, V, S, Y, CP, DLVLT, DLVLP)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF A FUEL, AIR, RESIDUAL
C GAS MIXTURE
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C F - RESIDUAL MASS FRACTION
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/g)
C S - ENTROPY (J/G/K)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2, 2=H2O, 3=N2, 4O2, 5=CO, 6=H2
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C
C REMARKS:
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE
C 3. ENTHALPIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. THE FUEL MOLE FRACTION IS UNITY MINUS THE SUM OF Y(I)
C
REAL M, MW, MRES, MFA, MFUEL, K, NU, N2
LOGICAL RICH, LEAN
DIMENSION A(7,6), TABLE(6), M(6), NU(6), Y(6), CP0(6), H0(6),
& S0(6)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./,
& A0/4.0652/, B0/6.0977E-02/, C0/-1.8801E-05/,
& D0/-3.5880E+04/, E0/15.45/
C
C TABLE 3.1, I = 1,6
C
DATA A/
1 .24007797E+01, .87350957E-02, -.6607878E-05, .20021861E-08,
1 .63274039E-15, -.48377527E+05, .96951457E+01,
2 .40701275E+01, -.11084499E-02, .41521180E-05, -.29637404E-08,
2 .80702103E-12, -.30279722E+05, -.32270046E+00,
3 .36748261E+01, -.12081500E-02, .23240102E-05, -.63217559E-09,
3 -.22577253E-12, -.10611588E+04, .23580424E+01,
4 .36255985E+01, -.18782184E-02, .70554544E-05, -.67635137E-08,
4 .21555993E-11, -.10475226E+04, .43052778E+01,
5 .37100928E+01, -.16190964E-02, .36923594E-05, -.20319674E-08,
5 .23953344E-12, -.14356310E+05, .29555355E+01,
6 .30574451E+01, .26765200E-02, -.58099162E-05, .55210391E-08,
6 -.18122739E-11, -.98890474E+03, -.22997056E+01/
C
C OTHER DATA
C
DATA RU/8.315/, TABLE/-1.,1.,0.,0.,1.,-1./, M/44.01, 18.02,
& 28.008, 32.000, 28.01, 2.018/
C
C COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
C
RICH = PHI .GT. 1.0
LEAN = .NOT. RICH
DLVLT = 1.0
DLVLP = -1.0
EPS = .21/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF (RICH) GO TO 10
NU(1) = ALPHA*PHI*EPS
NU(2) = BETA*PHI*EPS/2
NU(3) = 0.79 + DELTA*PHI*EPS/2
NU(4) = 0.21*(1.0 - PHI)
NU(5) = 0.
NU(6) = 0.
DCDT = 0.
GO TO 20
10 Z = 1000./T
K = EXP(2.743 + Z*(-1.761 + Z*(-1.611 + Z*.2803)))
DKDT = -K*(-1.761 + Z*(-3.222 + Z*.8409))/1000.
A1 = 1.0 - K
B = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + K*(.42*(PHI - 1.) +
& ALPHA*PHI*EPS)
C = -.42*ALPHA*PHI*EPS*(PHI - 1.0)*K
NU(5) = (-B + SQRT(B*B - 4.*A1*C))/2./A1
DCDT = DKDT*(NU(5)**2 - NU(5)*(.42*(PHI - 1.) + ALPHA*
& PHI*EPS) + .42*ALPHA*PHI*EPS*(PHI - 1.))/
& (2.*NU(5)*A1 + B)
NU(1) = ALPHA*PHI*EPS - NU(5)
NU(2) = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + NU(5)
NU(3) = .79 + DELTA*PHI*EPS/2.
NU(4) = 0.
NU(6) = .42*(PHI - 1.) - NU(5)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
C
20 TMOLES = 0.
DO 30 I = 1,6
30 TMOLES = TMOLES + NU(I)
MRES = 0.
DO 40 I = 1,6
Y(I) = NU(I)/TMOLES
40 MRES = MRES + Y(I)*M(I)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF FURLE-AIR
C
FUEL = EPS*PHI/(1. + EPS*PHI)
O2 = .21/(1. + EPS*PHI)
N2 = .79/(1. + EPS*PHI)
MFA = FUEL*(12.01*ALPHA + 1.008*BETA + 16.*GAMMA + 14.01*
& DELTA) + 32.*O2 + 28.02*N2
C
C COMPUTE MOLE FRACTIONS OF FUEL-AIR-RESIDUAL GAS
C
YRES = F/(F + MRES/MFA*(1. - F))
DO 50 I = 1,6
50 Y(I) = Y(I)*YRES
YFUEL = FUEL*(1. - YRES)
Y(3) = Y(3) + N2*(1. - YRES)
Y(4) = Y(4) + O2*(1. - YRES)
C
C COMPUTE COMPONENT PROPERTIES
C
DO 60 I = 1,6
CP0(I) = A(1,I) + A(2,I)*T + A(3,I)*T**2 + A(4,I)*T**3 +
& A(5,I)*T**4
H0(I) = A(1,I) + A(2,I)/2.*T + A(3,I)/3.*T**2 + A(4,I)/4.*
& T**3 + A(5,I)/5.*T**4 + A(6,I)/T
60 S0(I) = A(1,I)*ALOG(T) + A(2,I)*T + A(3,I)/2.*T**2 +
& A(4,I)/3.*T**3 + A(5,I)/4.*T**4 + A(7,1)
C
MFUEL = 12.01*ALPHA + 1.008*BETA + 16.000*GAMMA + 14.01*DELTA
CPFUEL = A0 + B0*T + C0*T**2
HFUEL = A0 + B0/2.*T + C0/3.*T**2 + D0/T
S0FUEL = A0*ALOG(T) + B0*T + C0/2.*T**2 + E0
C
C COMPUTE PROPERTIES OF MIXTURE
C
H = HFUEL*YFUEL
S = (S0FUEL - ALOG(YFUEL))*YFUEL
CP = CPFUEL*YFUEL
MW = MFUEL*YFUEL
DO 70 I = 1,6
H = H + H0(I)*Y(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
CP = CP + CP0(I)*Y(I) + H0(I)*T*TABLE(I)*DCDT*YRES/TMOLES
70 MW = MW + Y(I)*M(I)
C
R = RU/MW
H = R*T*H
U = H - R*T
V = 10.*R*T/P
S = R*(-ALOG(.9869*P) + S)
CP = R*CP
C
RETURN
END
C
C
SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W)
INTEGER N, IND, NW, K
DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
C Added external statement for fcn to avoid a warning message.
external fcn
C
C***********************************************************************
C *
C PURPOSE - THIS IS A RUNGE-KUTTA SUBROUTINE BASED ON VERNER'S *
C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO *
C THE SOLUTION OF A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL *
C EQUATIONS WITH INITIAL CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL *
C ERROR PROPORTIONAL TO A TOLERANCE SPECIFIED BY THE USER. (THE *
C PROPORTIONALITY DEPENDS ON THE KIND OF ERROR CONTROL THAT IS USED, *
C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) *
C *
C VARIOUS OPTIONS ARE AVAILABLE TO THE USER, INCLUDING DIFFERENT *
C KINDS OF ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS *
C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE CALCULATION (AND *
C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES. *
C *
C THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS. HOWEVER, A GOOD *
C VARIABLE-ORDER-ADAMS METHOD WILL PROBABLY BE MORE EFFICIENT IF THE *
C FUNCTION EVALUATIONS ARE VERY COSTLY. SUCH A METHOD WOULD ALSO BE *
C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE *
C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE FOR EXAMPLE *
C WITH GRAPHICAL OUTPUT. *
C *
C HULL-ENRIGHT-JACKSON 1/10/76