×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Contact US

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Run-time error in heat release rate program

Run-time error in heat release rate program

Run-time error in heat release rate program

(OP)
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

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members! Already a Member? Login


Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close