Home > Archive > Fortran > December 2005 > Too high values(T,P) after combustion in fortran code
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
Too high values(T,P) after combustion in fortran code
|
|
|
| Hi,
I am getting a too much high value in Temperature & Pressure after
combustion and as I am not very experienced in Fortran.
The program has been lifted out of a book called "Interneal Combustion
Engines:Applied Thermosciences" and is a program to model the otto cycle
(just compression and combustion).
The main code is listed below and any help in resolving this matter would
be greatly appreciated.
Sorry for my poor english. :)
Thanks
James,
****************************************
*************
DIMENSION Y1(6), Y2(10)
DATA AO/47870./, FS/0.06548/
C
R=10.
PHI=0.8
F=0.1
P1=1.0
T1=350.
C
CALL FARG(P1,T1,PHI,F,H1,U1,V1,S1,Y1,CP,DLVLT
,DLVLP)
V2=V1/R
S2=S1
GAMMA=CP/(CP+P1*V1/T1/10.*DLVLT**2./DLVLP)
T2=T1*(V1/V2)**(GAMMA-1.0)
P2=P1*(V1/V2)**(GAMMA)
CALL CMPRSS(V2,S2,PHI,R,F,T2,P2,U2)
V3=V2
U3=U2
P3=200.
T3=4000.
CALL COMBST(V3,U3,PHI,P3,T3,S3)
STOP
END
SUBROUTINE COMBST(V3,U3,PHI,P3,T3,S3)
DIMENSION Y2(10)
DATA MAXITS/650/, TOL/.0001/, DX3/1000./, DX4/100./
U2 = U3
V2 = V3
DO 620 I=1, MAXITS
CALL ECP(P3,T3,PHI,H3,U3,V3,S3,Y2,CP,DLVLT,DL
VLP,IER)
IF(IER.NE.O) WRITE(6,630) IER
630 FORMAT('ECP RETURNED IER=',I3)
F1=U2-U3
F2=V2-V3
DET=V3/P3*DLVLP*(CP-P3*V3/T3*DLVLT)+V3**2/T3*DLVLT*(DLVLT+DLVLP)
DX3=(F2*V3*(DLVLT+DLVLP)+F1*V3/P3*DLVLP)/DET
DX4=(-V3/T3*DLVLT*F1+F2*(CP-P3*V3/T3*DLVLT))/DET
T3=T3+DX3
P3=P3+DX4
IF(ABS(DX3)/T3.LT.TOL.AND.ABS(DX4)/P3.LT.TOL) GO TO 640
PRINT *,T3,P3,S3
PAUSE
620 CONTINUE
WRITE(6,650)
650 FORMAT('CONVERGENCE FAILURE')
640 CONTINUE
C
RETURN
END
SUBROUTINE CMPRSS(V2,S2,PHI,R,F,T2,P2,U2)
C IMPLICIT REAL (A-H,O-Z)
C IMPLICIT INTEGER (I-N)
C REAL IMEP
DIMENSION Y1(6)
DATA MAXITS/550/, TOL/.0001/, DX1/2000./, DX2/100./
S1 = S2
V1 = V2
C DO ITERATION
C
DO 520 I=1, MAXITS
CALL FARG(P2,T2,PHI,F,H2,U2,V2,S2,Y1,CP,DLVLT
,DLVLP)
IF(ABS(DX1)/T2.LT.TOL.AND.ABS(DX2)/P2.LT.TOL) GO TO 540
F1=S1-S2
F2=V1-V2
DET=(CP/P2*DLVLP+V2/10./T2*DLVLT**2)
DX1=(T2/P2*DLVLP*F1+DLVLT*F2/10.)/DET
DX2=(-DLVLT*F1+CP/V2*F2)/DET
T2=T2+DX1
P2=P2+DX2
PRINT*,T2,P2
PAUSE
520 CONTINUE
WRITE(6,550)
550 FORMAT ('CONVERGENCE FAILURE')
540 CONTINUE
RETURN
END
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 VOLUNE (CM**3/G)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTION
C 1=CO2, 2=H2O, 3=N2, 4=O2, 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 PRESSUREE AT CONSTANT TEMPERATURE
C
C REMARKS :
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE(C7H17). TO CHANGE FUELS
C ONLY THE FUEL DATA STATEMENT NEEDS MODIFICATION
C 3. ENTHALIES 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
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..... FUEL DATA
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./,
& A0/4.0652/, B0/6.0977E-02/, C0/-1.8801E-5/, D0/-3.588E+04/,
& E0/15.45/
C..... TABLE 3.1, I=1,6
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, .2955535E+01,
6 .30574451E+01, .26765200E-02, -.58099162E-05, .55210391E-08,
6 -.18122739E-11, -.98890474E+03, -.22997056E+01/
C..... OTHER DATA
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..... COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
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 + .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..... COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
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..... COMPUTE MOLE FRACTION AND MOLECULAR WEIGHT OF FUEL-AIR
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..... COMPUTE MOLE FRACTION OF FUEL-AIR-RESIDUAL GAS
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..... COMPUTE COMPONENT PROPERITIES
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,I)
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..... COMPUTE PROPERTIES OF MIXTURE
H = HFUEL*YFUEL
IF(YFUEL.EQ.0.) THEN
S = 0.
GO TO 61
ENDIF
S = (S0FUEL - ALOG(YFUEL))*YFUEL
61 CP = CPFUEL*YFUEL
MW = MFUEL*YFUEL
DO 70 I = 1,6
H = H + H0(I)*Y(I)
IF( Y(I).EQ.0. ) GO TO 65
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
65 CP = CP + CP0(I)*Y(I) + H0(I)*T*TABLE(I)*DCDT*YRES/TMOLES
70 MW = MW + Y(I)*M(I)
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
RETURN
END
SUBROUTINE ECP(P, T, PHI, H, U, V, S, Y, CP, DLVLT, DLVLP, IER)
C
C PURPOSE :
C COMPUTE THE PROPERTIES OF EQULIBRIUM COMBUSTION
C PRODUCT
C
C GIVEN :
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RETIO
C
C RETURNS :
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUNE (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 PRESSUREE 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
C REMARKS :
C 1. VALID FOR 300 < T < 4000
C 2. ASSUMES THE FUEL IS GASOLINE(C7H17). TO CHANGE FUELS
C ONLY THE FUEL DATA STATEMENT NEEDS MODIFICATION
C 3. ENTHALIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. SOUBROUTINE REQUIRED : FARG, GAUSSJ
C 5. IF ESTIMATES OF Y(I) I=1,10 ARE AVAILALE THEY SHOULD
C BE INPUT THROUGH THE ARGUMENT LIST. OTHERWISE INPUT
C ZERO'S AND THE ROUTINE WILL MAKE ITS OWN ESTIMATES
C
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)
C.... FUEL DATA
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0.0/, DELTA/0.0/
C.... SPECIFIC HEAT DATA FROM GORDON AND MCBRIDE (1971)
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, .73618264E-03, -.19652228E-06, .36201558E-10,
4 -.28945627E-14, -.12019825E+04, .36150960E+01,
5 .29840696E+01, .14891390E-02, -.57899684E-06, .10364577E-09,
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 .25420596E+01, -.27550619E-04, -.31028033E-08, .45510674E-11,
8 -.43680515E-15, .29230803E+05, .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.... EQILIBRIUM CONSANT DATA FROM OLIKARA AND BORMAN (1975)
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, -.260286E+01, .259556E-03,
5 -.162687E-07,
6 -.415302E-02, .148627E+05, -.475746E+01, .124699E-03,
6 -.900227E-08/
C.... MOLECULAR WEIGHTS
DATA M/44.01, 18.02, 28.008, 32., 28.01, 2.018, 1.009, 16.,
& 17.009, 30.004/
C.... OTHER DATA
DATA MAXITS/100/, TOL/3.0E-05/, RU/8.31434/
C.... MAKE SURE SOLID CARBON EILL NOT FORM
IER = 2
EPS = .210/(ALPHA + .25*BETA - .5*GAMMA)
IF( PHI .GT. (.210/EPS/(.5*ALPHA - .5*GAMMA)) ) RETURN
IER = 0
C.... DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C.... REQUIRED OR IF T < 1000 IN WHICH CASE FARG WILL SUFFICE
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.... EVALUATE THE REQUISITE CONSTANTS
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.... SET UP ITERATION LOOP FOR NEWTON-RAPSON ITERATION AND
C.... INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C.... AT LOW TEMPERATURE ON 24 BIT MACHINES (VAX). SET
C.... TRICK = 1.0 ON 48 BIT MACHINES (CDC)
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 + D94
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 + D94 + 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.... SOLVE THE MATRIX EQUATION 3.81 FOR COMPOSITION CORRECTIONS
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 GAUSSJ(A, 4, B)
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) = AMIN1(1.0, AMAX1(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.... COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
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) = -.5*C(1)/P
DCDP(2) = -.5*C(2)/P
DCDP(5) = .5*C(5)/P
DCDP(6) = .5*C(6)/P
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.... SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C.... AND THEN DETERMINE THE DEPENDENT DERIVATIVES
DO 70 I = 1,4
70 B(I) = -DFDT(I)
CALL GAUSSJ(A, 4, B)
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)
DO 80 I = 1,4
80 B(I) = -DFDP(I)
CALL GAUSSJ(A, 4, B)
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.... COMPUTE THE THERMODYNAMIC PROPERTIES
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,I)/5.*T**4 + A0(6,I)/T
90 S0(I) = A0(1,I)*ALOG(T) + A0(2,I)*T + A0(3,I)/2.*T*T
& + A0(4,I)/3.*T**3 + A0(5,I)/4.*T**4 + A0(7,I)
C.... Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C.... CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
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.
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)
IF(Y(I).EQ.0.) GO TO 100
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
100 CONTINUE
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) + S)
U = H - R*T
RETURN
END
C
SUBROUTINE GAUSSJ(A, N, B)
C
C PURPOSE :
C LINEAR EQUATION SOLUTION BY GAUSS-JORDAN ELIMINATION, AX = B.
C
C GIVEN :
C A - AN INPUT MATRIX OF N BY N ELEMENTS.
C B - AN INPUT MATRIX OF N.
C
C RETURN :
C ON OUTPUT ONLY B IS REPLACED BY THE CORRESPONDING SET
C OF SOLUTION VECTORS.
C
C***************************************
*****************************
PARAMETER (NMAX = 50)
DIMENSION A(N,N), AA(NMAX,NMAX), B(N),
& IPIV(NMAX), INDXR(NMAX), INDXC(NMAX)
DO 11 J=1,N
IPIV(J) = 0
DO 10 JJ=1,N
10 AA(J,JJ) = A(J,JJ)
11 CONTINUE
DO 22 I=1,N
BIG = 0.
DO 13 J=1,N
IF(IPIV(J).NE.1) THEN
DO 12 K=1,N
IF(IPIV(K).EQ.0) THEN
IF(ABS(AA(J,K)).GE.BIG) THEN
BIG = ABS(AA(J,K))
IROW = J
ICOL = K
ENDIF
ELSE IF (IPIV(K).GT.1) THEN
PAUSE ' ERROR ---> SINGULAR MATRIX'
ENDIF
12 CONTINUE
ENDIF
13 CONTINUE
IPIV(ICOL) = IPIV(ICOL) + 1
IF(IROW.NE.ICOL) THEN
DO 14 L=1,N
DUM = AA(IROW,L)
AA(IROW,L) = AA(ICOL,L)
AA(ICOL,L) = DUM
14 CONTINUE
DUM = B(IROW)
B(IROW) = B(ICOL)
B(ICOL) = DUM
ENDIF
INDXR(I) = IROW
INDXC(I) = ICOL
IF(AA(ICOL,ICOL).EQ.0) PAUSE ' ERROR ---> SINGULAR MATRIX'
PIVINV = 1./AA(ICOL,ICOL)
AA(ICOL,ICOL) = 1.
DO 16 L=1,N
AA(ICOL,L) = AA(ICOL,L) * PIVINV
16 CONTINUE
B(ICOL) = B(ICOL) * PIVINV
DO 21 LL=1,N
IF(LL.NE.ICOL) THEN
DUM = AA(LL,ICOL)
AA(LL,ICOL) = 0.
DO 18 L=1,N
AA(LL,L) = AA(LL,L) - AA(ICOL,L)*DUM
18 CONTINUE
B(LL) = B(LL) - B(ICOL)*DUM
ENDIF
21 CONTINUE
22 CONTINUE
DO 24 L=N,1,-1
IF(INDXR(L).NE.INDXC(L)) THEN
DO 23 K=1,N
DUM = AA(K,INDXR(L))
AA(K,INDXR(L)) = AA(K,INDXC(L))
AA(K,INDXC(L)) = DUM
23 CONTINUE
ENDIF
24 CONTINUE
RETURN
END
| |
| Tim Prince 2005-12-13, 7:04 pm |
| cjn99 wrote:
> Hi,
>
> I am getting a too much high value in Temperature & Pressure after
> combustion and as I am not very experienced in Fortran.
>
> The program has been lifted out of a book called "Interneal Combustion
> Engines:Applied Thermosciences" and is a program to model the otto cycle
> (just compression and combustion).
>
> IF(IER.NE.O) WRITE(6,630) IER
When working on a code of unknown quality, you should turn on the
warnings provided by your compiler, e.g.
gfortran -Wall -O -march=pentium-m james_7.f
....
james_7.f: In function ‘combst’:
james_7.f:32: warning: ‘o’ is used uninitialized in this function
This is a clue that "O" may be a typo, where "0" was meant. Although
there are options in many compilers which happen to make them
synonymous, this is bound to make it unreliable.
....
james_7.f: In function ‘gaussj’:
james_7.f:556: warning: ‘icol’ may be used uninitialized in this function
james_7.f:555: warning: ‘irow’ may be used uninitialized in this function
This tells you that the compiler will have to make its own decision how
to treat these variables, in case the data values are not as the program
author expected.
I could be nearly as critical of the excessive use of default data
typing (names I.. thru N.. integer, else real).
How are we to know what results you expect? The units are rather
unfamiliar, even to some of us who studied this stuff at one time. It
would take excessive study simply to figure out what problem is intended
to be solved.
I had difficulty with the formatting of your fixed form data. I ran it
through some old-fashioned reformatting:
C***************************************
*************
real y1(6),y2(10)
data ao/47870./,fs/0.06548/
C
r= 10.
phi= 0.8
f= 0.1
p1= 1.0
t1= 350.
C
call farg(p1,t1,phi,f,h1,u1,v1,s1,y1,cp,dlvlt
,dlvlp)
v2= v1/r
s2= s1
gamma= cp/(cp+p1*v1/t1/10.*dlvlt**2./dlvlp)
t2= t1*(v1/v2)**(gamma-1.0)
p2= p1*(v1/v2)**(gamma)
call cmprss(v2,s2,phi,r,f,t2,p2,u2)
v3= v2
u3= u2
p3= 200.
t3= 4000.
call combst(v3,u3,phi,p3,t3,s3)
stop
end
subroutine combst(v3,u3,phi,p3,t3,s3)
real y2(10)
data maxits/650/,tol/.0001/,dx3/1000./,dx4/100./
u2= u3
v2= v3
do i= 1,maxits
call ecp(p3,t3,phi,h3,u3,v3,s3,y2,cp,dlvlt,dl
vlp,ier)
if(ier.ne.o)then
write(6,10)ier
endif
f1= u2-u3
f2= v2-v3
det= v3/p3*dlvlp*(cp-p3*v3/t3*dlvlt)+v3**2/t3*dlvlt*(dlvlt+dlv
&lp)
dx3= (f2*v3*(dlvlt+dlvlp)+f1*v3/p3*dlvlp)/det
dx4= (-v3/t3*dlvlt*f1+f2*(cp-p3*v3/t3*dlvlt))/det
t3= t3+dx3
p3= p3+dx4
if(abs(dx3)/t3 < tol.and.abs(dx4)/p3 < tol)then
return
endif
print*,t3,p3,s3
pause
enddo
write(6,20)
return
10 format('ECP RETURNED IER=',i3)
20 format('CONVERGENCE FAILURE')
end
subroutine cmprss(v2,s2,phi,r,f,t2,p2,u2)
C IMPLICIT REAL (A-H,O-Z)
C IMPLICIT INTEGER (I-N)
C REAL IMEP
real y1(6)
data maxits/550/,tol/.0001/,dx1/2000./,dx2/100./
s1= s2
v1= v2
C DO ITERATION
C
do i= 1,maxits
call farg(p2,t2,phi,f,h2,u2,v2,s2,y1,cp,dlvlt
,dlvlp)
if(abs(dx1)/t2 < tol.and.abs(dx2)/p2 < tol)then
return
endif
f1= s1-s2
f2= v1-v2
det= (cp/p2*dlvlp+v2/10./t2*dlvlt**2)
dx1= (t2/p2*dlvlp*f1+dlvlt*f2/10.)/det
dx2= (-dlvlt*f1+cp/v2*f2)/det
t2= t2+dx1
p2= p2+dx2
print*,t2,p2
pause
enddo
write(6,30)
return
30 format('CONVERGENCE FAILURE')
end
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 VOLUNE (CM**3/G)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTION
C 1=CO2, 2=H2O, 3=N2, 4=O2, 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 PRESSUREE AT CONSTANT TEMPERATURE
C
C REMARKS :
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE(C7H17). TO CHANGE FUELS
C ONLY THE FUEL DATA STATEMENT NEEDS MODIFICATION
C 3. ENTHALIES 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
C***************************************
*******************************
integer m,mw,mres,mfa,mfuel,k,nu,n2
logical rich,lean
real a(7,6),table(6),y(6),cp0(6),h0(6),s0(6)
integer m(6),nu(6)
C..... FUEL DATA
data alpha/7.0/,beta/17.0/,gamma/0./,delta/0./,a0/4.0652/,b0/6.097
&7e-02/,c0/-1.8801e-5/,d0/-3.588e+04/,e0/15.45/
C..... TABLE 3.1, I=1,6
data a/.24007797e+01,.87350957e-02,-.6607878e-05,.20021861e-08,.63
&274039e-15,-.48377527e+05,.96951457e+01,.40701275e+01,-.11084499e-
&02,.41521180e-05,-.29637404e-08,.80702103e-12,-.30279722e+05,-.322
&70046e+00,.36748261e+01,-.12081500e-02,.23240102e-05,-.63217559e-0
&9,-.22577253e-12,-.10611588e+04,.23580424e+01,.36255985e+01,-.1878
&2184e-02,.70554544e-05,-.67635137e-08,.21555993e-11,-.10475226e+04
&,.43052778e+01,.37100928e+01,-.16190964e-02,.36923594e-05,-.203196
&74e-08,.23953344e-12,-.14356310e+05,.2955535e+01,.30574451e+01,.26
&765200e-02,-.58099162e-05,.55210391e-08,-.18122739e-11,-.98890474e
&+03,-.22997056e+01/
C..... OTHER DATA
data ru/8.315/,table/-1.,1.,0.,0.,1.,-1/,m/44.01,18.02,28.008,32.0
&00,28.01,2.018/
C..... COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
rich= phi > 1.0
lean= .not.rich
dlvlt= 1.0
dlvlp= -1.0
eps= .21/(alpha+0.25*beta-0.5*gamma)
if(.not.rich)then
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.
else
z= 1000./t
k= exp(2.743+z*(-1.761+z*(-1.611+.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*ep
&s)
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..... COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
endif
tmoles= 0.
do i= 1,6
tmoles= tmoles+nu(i)
enddo
mres= 0.
do i= 1,6
y(i)= nu(i)/tmoles
mres= mres+y(i)*m(i)
C..... COMPUTE MOLE FRACTION AND MOLECULAR WEIGHT OF FUEL-AIR
enddo
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..... COMPUTE MOLE FRACTION OF FUEL-AIR-RESIDUAL GAS
yres= f/(f+mres/mfa*(1.-f))
do i= 1,6
y(i)= y(i)*yres
enddo
yfuel= fuel*(1.-yres)
y(3)= y(3)+n2*(1.-yres)
y(4)= y(4)+o2*(1.-yres)
C..... COMPUTE COMPONENT PROPERITIES
do 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
s0(i)= a(1,i)*log(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,i)
enddo
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*log(t)+b0*t+c0/2.*t**2+e0
C..... COMPUTE PROPERTIES OF MIXTURE
h= hfuel*yfuel
if(yfuel == 0.)then
s= 0.
else
s= (s0fuel-log(yfuel))*yfuel
endif
cp= cpfuel*yfuel
mw= mfuel*yfuel
do i= 1,6
h= h+h0(i)*y(i)
if(y(i).ne.0.)then
s= s+y(i)*(s0(i)-log(y(i)))
endif
cp= cp+cp0(i)*y(i)+h0(i)*t*table(i)*dcdt*yre
s/tmoles
mw= mw+y(i)*m(i)
enddo
r= ru/mw
h= r*t*h
u= h-r*t
v= 10.*r*t/p
s= r*(-log(.9869*p)+s)
cp= r*cp
return
end
subroutine ecp(p,t,phi,h,u,v,s,y,cp,dlvlt,dlvlp,ier
)
C
C PURPOSE :
C COMPUTE THE PROPERTIES OF EQULIBRIUM COMBUSTION
C PRODUCT
C
C GIVEN :
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RETIO
C
C RETURNS :
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUNE (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 PRESSUREE 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
C REMARKS :
C 1. VALID FOR 300 < T < 4000
C 2. ASSUMES THE FUEL IS GASOLINE(C7H17). TO CHANGE FUELS
C ONLY THE FUEL DATA STATEMENT NEEDS MODIFICATION
C 3. ENTHALIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. SOUBROUTINE REQUIRED : FARG, GAUSSJ
C 5. IF ESTIMATES OF Y(I) I=1,10 ARE AVAILALE THEY SHOULD
C BE INPUT THROUGH THE ARGUMENT LIST. OTHERWISE INPUT
C ZERO'S AND THE ROUTINE WILL MAKE ITS OWN ESTIMATES
C
C***************************************
*******************************
real m,mw,k,kp,mt,mp
logical check
real m(10),k(6),c(6),y(10),d(3),b(4),a(4,4),y
0(6),kp(5,6),dfdt
& (4),dcdt(6),dkdt(6),dcdp(6),dydt(10),dyd
p(10),a0(7,10),dfdp(4),cp0
&(10),h0(10),s0(10)
C.... FUEL DATA
data alpha/7.0/,beta/17.0/,gamma/0.0/,delta/0.0/
C.... SPECIFIC HEAT DATA FROM GORDON AND MCBRIDE (1971)
data a0/.44608041e+01,.30981719e-02,-.12392571e-05,.22741325e-09,-
&.15525954e-13,-.48961442e+05,-.98635982e+00,.27167633e+01,.2945137
&4e-02,-.80224374e-06,.10226682e-09,-.48472145e-14,-.29905826e+05,.
&66305671e+01,.28963194e+01,.15154866e-02,-.57235277e-06,.99807393e
&-10,-.65223555e-14,-.90586184e+03,.61615148e+01,.36219535e+01,.736
&18264e-03,-.19652228e-06,.36201558e-10,-.28945627e-14,-.12019825e+
&04,.36150960e+01,.29840696e+01,.14891390e-02,-.57899684e-06,.10364
&577e-09,-.69353550e-14,-.14245228e+05,.63479156e+01,.31001901e+01,
&.51119464e-03,.52644210e-07,-.34909973e-10,.36945345e-14,-.8773804
&2e+03,-.19629421e+01,.25000000e+01,0.,0.,0.,0.,.25471627e+05,-.460
&11763e+00,.25420596e+01,-.27550619e-04,-.31028033e-08,.45510674e-1
&1,-.43680515e-15,.29230803e+05,.49203080e+01,.29106427e+01,.959316
&50e-03,-.19441702e-06,.13756646e-10,.14224542e-15,.39353815e+04,.5
&4423445e+01,.31890000e+01,.13382281e-02,-.52899318e-06,.95919332e-
&10,-.64847932e-14,.98283290e+04,.67458126e+01/
C.... EQILIBRIUM CONSANT DATA FROM OLIKARA AND BORMAN (1975)
data kp/.432168e+00,-.112464e+05,.267269e+01,-.745744e-04,.242484e
&-08,.310805e+00,-.129540e+05,.321779e+01,-.738336e-04,.344645e-08,
&-.141784e+00,-.213308e+04,.853461e+00,.355015e-04,-.310227e-08,.15
&0879e-01,-.470959e+04,.646096e+00,.272805e-05,-.154444e-08,-.75236
&4e+00,.124210e+05,-.260286e+01,.259556e-03,-.162687e-07,-.415302e-
&02,.148627e+05,-.475746e+01,.124699e-03,-.900227e-08/
C.... MOLECULAR WEIGHTS
data m/44.01,18.02,28.008,32.,28.01,2.018,1.009,16.,17.009,30.004/
C.... OTHER DATA
data maxits/100/,tol/3.0e-05/,ru/8.31434/
C.... MAKE SURE SOLID CARBON EILL NOT FORM
ier= 2
eps= .210/(alpha+.25*beta-.5*gamma)
if(phi <= .210/eps/(.5*alpha-.5*gamma))then
ier= 0
C.... DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C.... REQUIRED OR IF T < 1000 IN WHICH CASE FARG WILL SUFFICE
sum= 0.
do i= 1,10
sum= sum+y(i)
enddo
if(t <= 1000..or.sum <= 0.998)then
call farg(p,t,phi,1.,h,u,v,s,y0,cp,dlvlt,dlvlp)
do i= 1,10
y(i)= 0.
enddo
do i= 1,6
y(i)= y0(i)
enddo
if(t <= 1000.)then
return
endif
endif
patm= .9869233*p
C.... EVALUATE THE REQUISITE CONSTANTS
do i= 1,6
k(i)= 10.**(kp(1,i)*log(t/1000.)+kp(2,i)/t+kp(3,i)+kp(4,i)*
&t+kp(5,i)*t*t)
enddo
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.... SET UP ITERATION LOOP FOR NEWTON-RAPSON ITERATION AND
C.... INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C.... AT LOW TEMPERATURE ON 24 BIT MACHINES (VAX). SET
C.... TRICK = 1.0 ON 48 BIT MACHINES (CDC)
trick= 10.
check= abs(phi-1.0) < tol
if(check)then
phi= phi*(1.0+sign(tol,phi-1.0))
endif
icheck= 0
do j= 1,maxits
do jj= 1,10
y(jj)= min(1.0,max(y(jj),1.0e-25))
enddo
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+d94
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+d94+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.... SOLVE THE MATRIX EQUATION 3.81 FOR COMPOSITION CORRECTIONS
sum= 0.
do i= 1,10
sum= sum+y(i)
enddo
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 gaussj(a,4,b)
error= 0.
do l= 3,6
ll= l-2
y(l)= y(l)+b(ll)/trick
error= max(error,abs(b(ll)))
y(l)= min(1.0,max(y(l),1.0e-25))
enddo
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 < tol)then
icheck= icheck+1
endif
if(error < tol.and.icheck >= 2)then
goto40
endif
enddo
ier= 1
C.... COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
40 do i= 1,6
dkdt(i)= 2.302585*k(i)*(kp(1,i)/t-kp(2,i)/t/t+kp(4,i)+2.0*
&kp(5,i)*t)
enddo
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)= -.5*c(1)/p
dcdp(2)= -.5*c(2)/p
dcdp(5)= .5*c(5)/p
dcdp(6)= .5*c(6)/p
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)*x
&9+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.... SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C.... AND THEN DETERMINE THE DEPENDENT DERIVATIVES
do i= 1,4
b(i)= -dfdt(i)
enddo
call gaussj(a,4,b)
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)+d10
3*dydt(3)
do i= 1,4
b(i)= -dfdp(i)
enddo
call gaussj(a,4,b)
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.... COMPUTE THE THERMODYNAMIC PROPERTIES
do 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,i)/5.*t**4+a0(6,i)/t
s0(i)= a0(1,i)*log(t)+a0(2,i)*t+a0(3,i)/2.*t*t+a0(4,i)/3.
&*t**3+a0(5,i)/4.*t**4+a0(7,i)
C.... Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C.... CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
enddo
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.
mw= 0.
cp= 0.
h= 0.
s= 0.
mt= 0.
mp= 0.
do i= 1,10
if(y(i) > 1.0e-25)then
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)
if(y(i).ne.0.)then
s= s+y(i)*(s0(i)-log(y(i)))
endif
endif
enddo
r= ru/mw
v= 10.*r*t/p
cp= r*(cp-h*t*mt/mw)
dlvlt= 1.0+max(-t*mt/mw,0.)
dlvlp= -1.0-max(p*mp/mw,0.)
h= h*r*t
s= r*(-log(patm)+s)
u= h-r*t
endif
return
end
C
subroutine gaussj(a,n,b)
C
C PURPOSE :
C LINEAR EQUATION SOLUTION BY GAUSS-JORDAN ELIMINATION, AX = B.
C
C GIVEN :
C A - AN INPUT MATRIX OF N BY N ELEMENTS.
C B - AN INPUT MATRIX OF N.
C
C RETURN :
C ON OUTPUT ONLY B IS REPLACED BY THE CORRESPONDING SET
C OF SOLUTION VECTORS.
C
C***************************************
*****************************
parameter(nmax= 50)
real a(n,n),aa(nmax,nmax),b(n),ipiv(nmax),ind
xr(nmax),indxc(nm
&ax)
do j= 1,n
ipiv(j)= 0
do jj= 1,n
aa(j,jj)= a(j,jj)
enddo
enddo
do i= 1,n
big= 0.
do j= 1,n
if(ipiv(j).ne.1)then
do k= 1,n
if(ipiv(k).ne.0)then
if(ipiv(k) > 1)then
pause' ERROR ---> SINGULAR MATRIX'
endif
else
if(abs(aa(j,k)) >= big)then
big= abs(aa(j,k))
irow= j
icol= k
endif
endif
enddo
endif
enddo
ipiv(icol)= ipiv(icol)+1
if(irow.ne.icol)then
do l= 1,n
dum= aa(irow,l)
aa(irow,l)= aa(icol,l)
aa(icol,l)= dum
enddo
dum= b(irow)
b(irow)= b(icol)
b(icol)= dum
endif
indxr(i)= irow
indxc(i)= icol
if(aa(icol,icol) == 0)then
pause' ERROR ---> SINGULAR MATRIX'
endif
pivinv= 1./aa(icol,icol)
aa(icol,icol)= 1.
do l= 1,n
aa(icol,l)= aa(icol,l)*pivinv
enddo
b(icol)= b(icol)*pivinv
do ll= 1,n
if(ll.ne.icol)then
dum= aa(ll,icol)
aa(ll,icol)= 0.
do l= 1,n
aa(ll,l)= aa(ll,l)-aa(icol,l)*dum
enddo
b(ll)= b(ll)-b(icol)*dum
endif
enddo
enddo
do l= n,1,-1
if(indxr(l).ne.indxc(l))then
do k= 1,n
dum= aa(k,indxr(l))
aa(k,indxr(l))= aa(k,indxc(l))
aa(k,indxc(l))= dum
enddo
endif
enddo
return
end
| |
| Tim Prince 2005-12-13, 7:04 pm |
| Tim Prince wrote:
> cjn99 wrote:
>
>
>
>
>
>
>
> When working on a code of unknown quality, you should turn on the
> warnings provided by your compiler, e.g.
> gfortran -Wall -O -march=pentium-m james_7.f
> ...
> james_7.f: In function ‘combst’:
> james_7.f:32: warning: ‘o’ is used uninitialized in this function
> 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
> s0(i)= a(1,i)*log(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,i)
Look up Horner's rule for evaluation of polynomials.
| |
|
| >>> cp0(i)= a(1,i)+a(2,i)*t+a(3,i)*t**2+a(4,i)*t**3+
a(5,i)*t**4
>
> Look up Horner's rule for evaluation of polynomials.
Or look up "optimizing compiler" :)
--
FX
| |
| Paul Van Delst 2005-12-13, 7:04 pm |
| FX wrote:
>
>
> Or look up "optimizing compiler" :)
I don't think optimising compilers are going to solve your problems with respect to
accumulating errors using the above technique. Horner's method for evaluating polynomials
is usually the first thing that people learn when converting mathematics-on-paper to
mathematics-on-computer. In addition, most of these developed numerical methods come with
error analysis methods as well. All good stuff.
See, for example, ch.5 in N.J.Higham's book "Accuracy and stability of numerical
algorithms". I also just had a look in my Numerical Recipes book and I see that it states
(in Sec.5.3, p167):
"We assume that you know enough /never/ [emphasis theirs] to evaluate a polynomial this way:
p=c(1)+c(2)*x+c(3)*x**2+c(4)*x**3+c(5)*x
**4
Come the (computer) revolution, all persons found guilty of such criminal behaviour will
be summarily executed and their programs won't be!"
A bit tongue in ch , obviously, but the use of the technique in your quoted code is a
red flag for me that indicates poor quality code - especially when there is a very simple,
well known alternative. But (as you can probably tell) I'm anal about this sort of stuff
so keep your salt handy. :o)
cheers,
paulv
--
Paul van Delst
CIMSS @ NOAA/NCEP/EMC
| |
| Tim Prince 2005-12-13, 7:04 pm |
|
"Paul Van Delst" <Paul.vanDelst@noaa.gov> wrote in message
news:dnmt4i$l1u$1@news.nems.noaa.gov...
> FX wrote:
>
> I don't think optimising compilers are going to solve your problems with
> respect to accumulating errors using the above technique. Horner's method
> for evaluating polynomials is usually the first thing that people learn
> when converting mathematics-on-paper to mathematics-on-computer. In
> addition, most of these developed numerical methods come with error
> analysis methods as well. All good stuff.
>
> See, for example, ch.5 in N.J.Higham's book "Accuracy and stability of
> numerical algorithms". I also just had a look in my Numerical Recipes book
> and I see that it states (in Sec.5.3, p167):
>
> "We assume that you know enough /never/ [emphasis theirs] to evaluate a
> polynomial this way:
> p=c(1)+c(2)*x+c(3)*x**2+c(4)*x**3+c(5)*x
**4
> Come the (computer) revolution, all persons found guilty of such criminal
> behaviour will be summarily executed and their programs won't be!"
>
> A bit tongue in ch, obviously, but the use of the technique in your
> quoted code is a red flag for me that indicates poor quality code -
> especially when there is a very simple, well known alternative. But (as
> you can probably tell) I'm anal about this sort of stuff so keep your salt
> handy. :o)
>
Optimizing compilers should recover some of the time wasting implicit in the
quoted code, by taking advantage of common factors in the power
calculations. That is not sufficient to improve numerical properties. In
principle, the Fortran standard allows automatic translation to Horner's
method for some cases, but no other language does, and the machinery for
doing it is not widely available. In this case, with strict interpretation,
/3. would prevent full optimization. Differing evaluation of such
expressions between compilers is a common source of unexpected behavior.
| |
|
| > In this case, with strict interpretation,
> /3. would prevent full optimization.
This I don't understand, doesn't the fortran standard explicitly permit
to e.g. multiply by the reciprocal value instead of doing the division.
Joost
| |
| e p chandler 2005-12-13, 7:04 pm |
| Tim Prince wrote:
> cjn99 wrote:
>
>
>
>
> When working on a code of unknown quality, you should turn on the
> warnings provided by your compiler, e.g.
> gfortran -Wall -O -march=pentium-m james_7.f
> ...
> james_7.f: In function 'combst':
> james_7.f:32: warning: 'o' is used uninitialized in this function
> This is a clue that "O" may be a typo, where "0" was meant. Although
> there are options in many compilers which happen to make them
> synonymous, this is bound to make it unreliable.
> ...
> james_7.f: In function 'gaussj':
> james_7.f:556: warning: 'icol' may be used uninitialized in this function
> james_7.f:555: warning: 'irow' may be used uninitialized in this function
> This tells you that the compiler will have to make its own decision how
> to treat these variables, in case the data values are not as the program
> author expected.
>
In addition, using FTNCHEK, available at
http://www.dsm.fordham.edu/~ftnchek/
may help too. Running this on the OP's (reformatted to fixed form) code
brought up a number of warnings about un-used and un-set variables.
| |
| e p chandler 2005-12-13, 7:04 pm |
| Tim Prince wrote:
> Tim Prince wrote:
>
> Look up Horner's rule for evaluation of polynomials.
Sound advice. But I would do other things first before assuming that
the problem is one of numerical mis-analysis.
1. Bench check the code against known good test data.
2. Build some scaffolding code and test subroutines individually.
Excuse me for coining a corollary to Brooks' law:
"Adding optimization to buggy code makes it buggier."
-- Elliot
| |
| Mr Hrundi V Bakshi 2005-12-16, 9:58 pm |
| OP,
Look at Appendix A ("Program for Finding Coefficients of NASA Polynomials")
and Appendix C ("Table of Coefficient Sets for NASA Polynomials") in
"Combustion Chemistry", ed. W.C. Gardiner (Springer), QD516.G24 1984. If
you need either I can dig them up from the archieves and email them to you.
--
Bye,
Hrundi V.B.
|
|
|
|
|