c
c
c
c
c
SUBROUTINE GXNOX
C---------------------------------------------------------------
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C================================================================
C 1 CREK COMMON BLOCKS
C ====================
LOGICAL SLD
COMMON/CINDEX/JHCPS,JCK1,JITER,JCK2(15)
COMMON/CKINET/GCPSUM,GCK1(8)
LOGICAL LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
COMMON/CHMLGC/LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
COMMON/CSPECE/GH0(20),GCK2(20),GSMW(20),GCK3(300)
COMMON/CPARAM/GEMV,GHSUB0,JNS,GPA,GQ0,GQ1,GQ2,
1 GQ3,GQ4,GRHOP,GSM,GS1(20),GS2(20),GTK
C================================================================
SAVE GSNAME
SAVE TEN1, TEN2, GBX1, GBX2, TACT1, TACT2, GFSN,
1 L0SPR2, L0SPR3, L0SPR4, RSPR, RSOUR, ID
SAVE JXDBCK, JYDBCK, JZDBCK, JSDBCK, NS, NR
SAVE CRKCAL, CRKBUG, DBGNOX
SAVE GN2PR, GN2OX, GCO2PR, GH2OPR, GHM, GOM, GNM,
1 GNOM, GN2M, GOHM, GO2M, GCO2M, GH2OM, GASCON,
1 GCOM, GH2M, GFUM, SETMIN
SAVE JXNO, JXN, JXO, JXH, JXOH, JNOSR, JYOXID,
1 JXO2, JYO2, JYCO2, JYH2, JYH2O, JYCO, JEQUI,
1 JDEGF, JYN2, JYPR, LUREAD, IERR, JTMP
SAVE JCRKT, L0VOL, GFLOT, GFNOT, GFNOVT, GFLNT,
1 GEXC, GEXH, GEXN, GEXO
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
PARAMETER (NSM=10,NRM=3)
CHARACTER*4 GSNAME
DIMENSION GSNAME(8)
DIMENSION RSOUR(NSM),RSPR(NSM)
DIMENSION GBX1(NRM),GBX2(NRM),TACT1(NRM),TACT2(NRM)
DIMENSION TEN1(NRM),TEN2(NRM),ID(4,NRM),GFSN(NSM)
C
LOGICAL CRKCAL,QEQ,CRKBUG,DBGNOX
C
DATA GSNAME/'CO ','CO2','H ','H2 ','H2O','O ','OH ','O2'/
C Zeldovich rate constants
C GBX1 & GBX2 are in units of kg**2/m**3/kg-mol/s when TEN1
C and TEN2 are zero.
DATA GBX1/7.6E10,6.4E6,6.3E8/
DATA GBX2/1.6E10,1.5E6,2.5E9/
C TACT1 & TACT2 are in units of K ( i.e. E/Ro ) where
C is the activation energy & Ro is the gas constant.
DATA TACT1/3.8E4,3150.0,0.0/
DATA TACT2/0.0,1.95E4,2.446E4/
C TEN1 & TEN2 are exponents for the temperature
DATA TEN1/0.0,1.0,0.5/
DATA TEN2/0.0,1.0,0.5/
DATA ID/4,5,3,2,
1 7,2,3,5,
1 6,2,3,1/
C
NAMSUB='GXNOX'
IXL=IABS(IXL)
C Molecular mass in Kg/(Kg mol)
C*****************************************************************
C--- GROUP 1. Run title and other preliminaries
C*****************************************************************
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
C User may here change message transmitted to the VDU screen or
C batch-run log file.
CALL WRYT40('ground file is gxnox.f of 25.04.99 ')
CALL WRYT40('PHOENICS version number is 3.2 ')
C
C.... Read information from Q1
C integers
CALL RQ1I('NOXD','NS',NS)
CALL RQ1I('NOXD','NR',NR)
IF(LSG63) THEN
CALL RQ1I('NOXD','JXDBCK',JXDBCK)
CALL RQ1I('NOXD','JYDBCK',JYDBCK)
CALL RQ1I('NOXD','JZDBCK',JZDBCK)
CALL RQ1I('NOXD','JSDBCK',JSDBCK)
C reals
CALL RQ1R('NOXD','GN2OX',GN2OX)
CALL RQ1R('NOXD','GN2PR',GN2PR)
CALL RQ1R('NOXD','GCO2PR',GCO2PR)
CALL RQ1R('NOXD','GH2OPR',GH2OPR)
C logicals
CALL RQ1L('NOXD','CRKBUG',CRKBUG)
CALL RQ1L('NOXD','DBGNOX',DBGNOX)
ENDIF
C
CALL SUB3R(GHM,1.00797,GOM,15.9994,GNM,14.0067)
CALL SUB2R(GNOM,30.0061,GN2M,28.0134)
CALL SUB3R(GOHM,17.00737,GO2M,31.9988,GCO2M,44.00995)
CALL SUB3R(GH2OM,18.01534,GASCON,8314.0,SETMIN,1.0E-10)
C
IF(LSG62) THEN
CALL SUB2R(GCOM,28.01055,GH2M,2.01594)
ELSEIF(LSG63) THEN
GFUM=16.04303
ENDIF
C
CALL SUB2(JXNO,LBNAME('XNO'),JXN,LBNAME('XN'))
CALL SUB2(JXO,LBNAME('XO'),JXH,LBNAME('XH'))
CALL SUB2(JXOH,LBNAME('XOH'),JNOSR,LBNAME('NOSR'))
CALL SUB2(JCRKT,LBNAME('CRKT'),JXO2,LBNAME('XO2'))
IF(LSG62) THEN
CALL SUB2(JYCO2,LBNAME('YCO2'),JYH2,LBNAME('YH2'))
CALL SUB2(JYH2O,LBNAME('YH2O'),JYCO,LBNAME('YCO '))
CALL SUB2(JEQUI,LBNAME('EQUI'),JDEGF,LBNAME('DEGF'))
CALL SUB2(JYN2,LBNAME('YN2'),JYO2,LBNAME('YO2'))
JTMP=LBNAME('TMP1')
ELSEIF(LSG63) THEN
JTMP=LBNAME('TMP1')
JYPR=LBNAME('PROD')
JYOXID=LBNAME('OXID')
ENDIF
C
CRKCAL=.TRUE.
LUREAD=20
IERR=-1
CALL OPENZZ(29,IERR)
IF(IERR.NE.0) THEN
CALL WRIT40('ERROR OPENING CH4DAT IN GXNOX ')
CALL WRIT1I('ERROR NO',IERR)
CALL SET_ERR(252, 'ERROR OPENING CH4DAT IN GXNOX.',1)
C CALL EAROUT(1)
ENDIF
CALL GXCRK0( LUREAD , LUPR1 )
GEMV=0.0
LEQUIL=.TRUE.
CALL GXMAKE(L0SPR2,NX*NY,'SPR2')
CALL GXMAKE(L0SPR3,NX*NY,'SPR3')
CALL GXMAKE(L0SPR4,NX*NY,'SPR4')
C*****************************************************************
C--- GROUP 13. Boundary conditions and special sources
C*****************************************************************
ELSEIF(IGR.EQ.13.AND.ISC.EQ.1) THEN
C--coefficient = GRND
IF(NPATCH.EQ.'NOXSR') CALL FN0(CO,-L0SPR2)
ELSEIF(IGR.EQ.13.AND.ISC.EQ.12) THEN
C--value = GRND
IF(NPATCH.EQ.'NOXSR') THEN
IF(INDVAR.EQ.JXNO) THEN
CALL FN0(VAL,-L0SPR3)
C... evaluate and store local NO source per unit volume of gas
CALL SUB4(L0CO,L0SPR2,L0VAL,L0SPR3,L0SOR,L0F(JNOSR),
1 L0NO,L0F(JXNO))
DO 1311 J=1,NX*NY
F(L0SOR+J)=F(L0CO+J)*(F(L0VAL+J)-F(L0NO+J))
1311 CONTINUE
ELSEIF(INDVAR.EQ.JXN) THEN
CALL FN0(VAL,-L0SPR4)
ENDIF
ENDIF
***************************************************************
C--- GROUP 19. SECTION 4 ---- START OF ITERATION.
C***************************************************************
ELSEIF(IGR.EQ.19.AND.ISC.EQ.4) THEN
IF(CRKCAL) THEN
CALL SUB3(L0XH,L0F(JXH),L0XO,L0F(JXO),L0XOH,L0F(JXOH))
CALL SUB3(L0XO2,L0F(JXO2),L0CRKT,L0F(JCRKT),L0VOL,L0F(VOL))
IF(LSG62) THEN
CALL SUB3(L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),L0YN2,L0F(JYN2))
CALL SUB4(L0YCO2,L0F(JYCO2),L0YH2O,L0F(JYH2O),L0YH2,
1 L0F(JYH2),L0FTMP,L0F(JTMP))
ELSEIF(LSG63) THEN
CALL SUB2(L0FYPR,L0F(JYPR),L0FTMP,L0F(JTMP))
L0OXID=L0F(JYOXID)
ENDIF
IF(L0FTMP.EQ.0) THEN
WRITE(14,*) 'Temperature is not stored. Execution ended.'
CALL SET_ERR(253, 'Error. Temperature is not stored.',1)
C CALL EAROUT(1)
ENDIF
C**Equilibrium Calculation
if(dbgrnd) then
call writst
call writbl
call writ40('Crek calculation starts for current slab')
call writ1i('IZ ',IZ)
call writbl
call writst
endif
DO 1946 JX=1,NX
DO 1946 JY=1,NY
J = JY + NY*(JX-1)
C--Skip crek call for blocked cells
IF(SLD(J)) GO TO 1946
IF(L0VOL.NE.0) THEN
IF(F(L0VOL+J).LE.TINY) GO TO 1946
ENDIF
C
C--Supply all species except N2 for equilibrium calculation
C H2O = H +OH ; CO2 = CO +0.5O2 ; O2 = 2O
C H2O = 0.5H2 + OH ; H2 = 2H
C
IF(LSG62) THEN
GYPR=1.0-F(L0YN2+J)
GYN2=F(L0YN2+J)
C
ELSEIF(LSG63) THEN
GCPR = F(L0FYPR +J)
GYPR = AMAX1(TINY,GCPR - GN2PR*GCPR)
GYN2 = GN2PR*GCPR+GN2OX*(F(L0OXID+J))
ENDIF
C
C--Skip crek call if only N2 present
IF(GYN2.GT.0.98) GO TO 1946
C
C--Initialise CREK arrays S1 & S2 with mol fractions
C
DO 1941 IS=1,JNS
C CO
IF(IS.EQ.1.AND.LSG62) THEN
GS1(IS)=F(L0YCO+J)/(GCOM*GYPR)
GS2(IS)=GS1(IS)
C CO2(coal)
ELSEIF(IS.EQ.2) THEN
IF(LSG62) THEN
GS1(IS)=F(L0YCO2+J)/(GCO2M*GYPR)
C CO2(gas)
ELSEIF(LSG63) THEN
GS1(IS)=GCPR*GCO2PR/(GCO2M*GYPR)
ENDIF
GS2(IS)=GS1(IS)
C H2
ELSEIF(IS.EQ.4.AND.LSG62) THEN
GS1(IS)=F(L0YH2+J)/(GH2M*GYPR)
GS2(IS)=GS1(IS)
C H2O(coal)
ELSEIF(IS.EQ.5) THEN
IF(LSG62) THEN
GS1(IS)=F(L0YH2O+J)/(GH2OM*GYPR)
C H2O(gas)
ELSEIF(LSG63) THEN
GS1(IS)=GCPR*GH2OPR/(GH2OM*GYPR)
ENDIF
GS2(IS)=GS1(IS)
C O2
ELSEIF(IS.EQ.8.AND.LSG62) THEN
GS1(IS)=F(L0YO2+J)/(GO2M*GYPR)
GS2(IS)=GS1(IS)
C H,O & OH species
ELSE
GS1(IS)=SETMIN
GS2(IS)=GS1(IS)
ENDIF
1941 CONTINUE
C
C Temperature, pressure & other CREK input parameters
GTK = F(L0FTMP+J)
GPA = PRESS0
JHCPS = 1
C Enter CREK routine to determine enthalpy
CALL HCPS
C
GHSUM=0.0
DO 1942 IS=1,JNS
GHSUM=GHSUM+GH0(IS)*GS2(IS)
1942 CONTINUE
C
GHSUB0=GHSUM*GASCON*GTK
C
C Enter CREK to perform equilibrium dissociation calculation
CALL GXCREK
C H mass fraction from CREK
F(L0XH+J)=GS2(3)*GHM*GYPR
C O mass fraction from CREK
F(L0XO+J)=GS2(6)*GOM*GYPR
C OH mass fraction from CREK
F(L0XOH+J)=GS2(7)*GOHM*GYPR
C O2 mass fraction from CREK
F(L0XO2+J)=GS2(8)*GO2M*GYPR
C T from CREK equilibrium calculation
F(L0CRKT+J)=GTK
C--Check total mass of array GS2 = 1.0
GFRAC=0.0
DO 1944 IS=1,JNS
GFRAC=GFRAC+GS2(IS)*GSMW(IS)
1944 CONTINUE
C
IF(GFRAC.LT.0.995) THEN
WRITE(LUPR1,*)' gs1 input to crek;',
1 ' gs2 output from crek'
WRITE(LUPR1,*)' mol fractions of',
1 'equilibrium species'
WRITE(LUPR1,*)' GS1',
1 ' GS2 GH0 '
DO 1945 IS=1,JNS
WRITE(LUPR1,951) IS,GSNAME(IS),
1 GS1(IS),GS2(IS),GH0(IS)
1945 CONTINUE
CALL WRITBL
WRITE(LUPR1,*)'sum of mass fractions : ',GFRAC
WRITE(LUPR1,*)'Temperature= ',GTK,' Hsub0 = ',GHSUB0
CALL WRITBL
IF(LSG62) THEN
WRITE(LUPR1,*) ' mass fractions',' including N2'
WRITE(LUPR1,*)' input',' output '
DO 1947 IS = 1 , JNS
WRITE(LUPR1,951) IS,GSNAME(IS),
1 GS1(IS)*GSMW(IS)*GYPR,
1 GS2(IS)*GSMW(IS)*GYPR
1947 CONTINUE
WRITE(LUPR1,953) F(L0YN2+J),F(L0YN2+J)
ENDIF
951 FORMAT(1X,I2,1X,A4,1P4E11.3)
953 FORMAT(1X,' 9 N2 ',1P2E11.3)
ENDIF
1946 CONTINUE
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) 'Crek calculation ended at slab = ',IZ
CALL WRITBL
CALL WRITST
IF(IZ.EQ.NZ) CRKCAL=.FALSE.
ENDIF
C**End of equilibrium calculation*******************************
C**Calculation of NOX Zeldovich reaction sources.
CALL SUB3(L0XO,L0F(JXO),L0XN,L0F(JXN),L0XOH,L0F(JXOH))
CALL SUB3(L0XH,L0F(JXH),L0XO2,L0F(JXO2),L0XNO,L0F(JXNO))
CALL SUB3(L0FEP2,L0SPR2,L0FEP3,L0SPR3,L0FEP4,L0SPR4)
L0RHO1=L0F(DEN1)
IF(LSG62) THEN
L0YN2=L0F(JYN2)
L0VOL=L0F(VOL)
ELSEIF(LSG63) THEN
L0FYPR=L0F(JYPR)
ENDIF
C
if(dbgrnd) then
call writ40('NOX calculation starts for current slab ')
call writ2i('isweep ',isweep,'iz ',iz)
call writbl
endif
DO 1301 JX=1,NX
DO 1301 JY=1,NY
J=JY+NY*(JX-1)
c Initialize arrays
F(L0FEP2+J) = 0.0
F(L0FEP3+J) = 0.0
F(L0FEP4+J) = 0.0
c-- Store mechanism values in GFSN array
c-- Convert mass fractions to mole fractions ( kg-mol/Kg )
c H mole fraction
GFSN(1)=F(L0XH+J)/GHM
c N mole fraction
GFSN(2)=F(L0XN+J)/GNM
c NO mole fraction
GFSN(3)=F(L0XNO+J)/GNOM
c N2 mole fraction
IF(LSG62) THEN
GFSN(4)=F(L0YN2 +J)/GN2M
ELSEIF(LSG63) THEN
GFSN(4)=F(L0FYPR +J)*GN2PR/GN2M
ENDIF
c O mole fraction
GFSN(5)=F(L0XO +J)/GOM
c OH mole fraction
GFSN(6)=F(L0XOH+J)/GOHM
c O2 mole fraction
GFSN(7)=F(L0XO2+J)/GO2M
c
c-- For blocked cells, skip the following.
IF(.NOT.SLD(J)) THEN
CRJ IF(F(L0VOL+J).GT.TINY) THEN
GTK=F((L0F(JCRKT))+J)+TINY
TKINV=1.0/GTK
GRHO2=F(L0RHO1+J)*F(L0RHO1+J)
c
c-- Zeldovich mechanism - source terms
DO 1303 IS=1,NS
GFSN(IS)=AMAX1( TINY , GFSN(IS) )
RSOUR(IS) = 0.0
RSPR(IS) = 0.0
1303 CONTINUE
c
DO 1304 IR=1,NR
c-- Forward & reverse reaction rate calculations
TK1=AMAX1(AMIN1((TACT1(IR)*TKINV),100.),-100.)
TK2=AMAX1(AMIN1((TACT2(IR)*TKINV),100.),-100.)
c-- Forward & backward rates in kg-mol/s/m**3
FRATE=GBX1(IR)*EXP(-TK1)*GTK**TEN1(IR)*GRHO2
RRATE=GBX2(IR)*EXP(-TK2)*GTK**TEN2(IR)*GRHO2
c-- Species selection for reaction, ir
J1=ID(1,IR)
J2=ID(2,IR)
J3=ID(3,IR)
J4=ID(4,IR)
FF=-FRATE*GFSN(J1)*GFSN(J2)+RRATE*GFSN(J3)
1 *GFSN(J4)
DSDF=-FRATE*(GFSN(J1)+GFSN(J2))-RRATE*(GFSN(J3)
1 +GFSN(J4))
c-- Store source for N & NO in array RSOUR
RSOUR(J1)=RSOUR(J1)+FF
RSOUR(J2)=RSOUR(J2)+FF
RSOUR(J3)=RSOUR(J3)-FF
RSOUR(J4)=RSOUR(J4)-FF
RSPR(J1)=RSPR(J1)+DSDF
RSPR(J2)=RSPR(J2)+DSDF
RSPR(J3)=RSPR(J3)+DSDF
RSPR(J4)=RSPR(J4)+DSDF
1304 CONTINUE
ENDIF
c-- Coefficient for N
F(L0FEP2+J)=-RSPR(2)
c-- Value for XN ( multiply by molecular weight so that final
c reaction rate source is in units of kg/m**3/s )
F(L0FEP4+J)=-(RSOUR(2)/(RSPR(2)+TINY))*GNM+F(L0XN+J)
c-- Value for XNO ( multiply by molecular weight so that final
c reaction-rate source is in units of kg/m**3/s )
F(L0FEP3+J)=RSOUR(3)*GNOM
1301 CONTINUE
C***************************************************************
C--- GROUP 19. SECTION 6 ---- Finish of IZ slab.
C***************************************************************
ELSEIF(IGR.EQ.19.AND.ISC.EQ.6) THEN
C
IF(LSG62) THEN
C-- Equivalence ratio calculation
IF(ISWEEP.EQ.LSWEEP-1.OR.ENUFSW) THEN
IF(JEQUI.NE.0) THEN
CALL SUB4(L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),L0YH2,
1 L0F(JYH2),L0YCO2,L0F(JYCO2))
CALL SUB2(L0YH2O,L0F(JYH2O),L0EQUI,L0F(JEQUI))
C-- Valence (C-4,O-2,H-1) / molecular mass; the division by mol
C mass is to convert mass frac in the loops below into moles.
CALL SUB4R(GE1,4./28.,GE2,4./44.,GE3,2./18.,GE4,2./2.)
CALL SUB4R(GE5,2./28.,GE6,4./44.,GE7,2./18.,GE8,4./32.)
C-- Equivalence ratio = ((4*moles C)+(1*moles H)) / (2*moles O)
DO 1961 J=1,NX*NY
GNUM=GE1*F(L0YCO+J)+GE2*F(L0YCO2+J)+GE3*F(L0YH2O+J)+
1 GE4*F(L0YH2+J)
GDEN=GE5*F(L0YCO+J)+GE6*F(L0YCO2+J)+GE7*F(L0YH2O+J)+
1 GE8*F(L0YO2+J)
GEQUIV=GNUM/(GDEN+TINY)
F(L0EQUI+J)=GEQUIV
1961 CONTINUE
ENDIF
IF(JDEGF.NE.0) THEN
C-- Calculate temperatures in Fahrenheit
L0DEGF=L0F(JDEGF)
C-- Convert crkt temperatures to fahrenheit
DO 1962 J=1,NX*NY
1962 F(L0DEGF+J)=(F((L0F(JCRKT))+J)-273.)*1.8+32.
ENDIF
ENDIF
C
C-- Calculate nox outlet data
CALL GETPTC('OUTLET',GTYPE,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
IF(IZ.LT.JZF .OR. IZ.GT.JZL) RETURN
IF(IZ.EQ.JZF) THEN
CALL SUB3R (GFLOT,0.0,GFLNT,0.0,GFNOT,0.0)
CALL SUB4R(GCOM,28.01055,GO2M,31.9988,GN2M,28.0134,
1 CO2M,44.00995)
CALL SUB3R(GH2M,2.01594,GH2OM,18.01534,GNOM,30.0061)
GFNOVT=0.
ENDIF
CALL SUB4(L0P1,L0F(P1),L0R1,L0F(R1),L0RHO,L0F(DEN1),
1 L0XN,L0F(JXN))
CALL SUB4(L0XNO,L0F(JXNO),L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),
1 L0YCO2,L0F(JYCO2))
CALL SUB3(L0YN2,L0F(JYN2),L0YH2,L0F(JYH2),L0YH2O,L0F(JYH2O))
IF(QEQ(GTYPE,6.0).OR.QEQ(GTYPE,7.0)) THEN
L0AREA=L0F(AHIGH)
ELSEIF (QEQ(GTYPE,4.0).OR.QEQ(GTYPE,5.0)) THEN
L0AREA=L0F(ANORTH)
ELSEIF (QEQ(GTYPE,2.0).OR.QEQ(GTYPE,3.0)) THEN
L0AREA=L0F(AEAST)
ELSE
WRITE (LUPR1,*) 'Outlet boundary is not an AREA type.'
WRITE (LUPR1,*) 'Outlet summary coding skipped.'
RETURN
ENDIF
CALL GETCOV('OUTLET',P1,GCO,GVAL)
IF(IZ.EQ.JZF) CALL SUB4R(GEXC,0.,GEXH,0.,GEXO,0.,GEXN,0.)
DO 1963 JX=JXF,JXL
DO 1963 JY=JYF,JYL
J=JY+(JX-1)*NY
GFLOW=-GCO*F(L0R1+J)*(GVAL-F(L0P1+J))*F(L0AREA+J)
GFLONO=GFLOW*F(L0XNO+J)
RMIXP=(F(L0YCO+J)/GCOM+F(L0YO2+J)/GO2M+F(L0YN2+J)/GN2M
1 +F(L0YCO2+J)/GCO2M+F(L0YH2+J)/GH2M
1 +F(L0YH2O+J)/GH2OM)*GASCON
GFLNOV=GFLONO* (GASCON/(GNOM*RMIXP))
GFLOT=GFLOT+GFLOW
GFNOT= GFNOT + GFLONO
GFNOVT= GFNOVT + GFLNOV
GEXC=GEXC+(F(L0YCO+J)*12./28.+F(L0YCO2+J)*12./44.)*GFLOW
GEXH=GEXH+(F(L0YH2+J)+F(L0YH2O+J)*2./18.)*GFLOW
GEXO=GEXO+(F(L0YCO+J)*16./28.+F(L0YCO2+J)*32./44.
1 +F(L0YH2O+J)*16./18.+F(L0YO2+J))*GFLOW
GEXN=GEXN+(F(L0YN2+J))*GFLOW
1963 CONTINUE
IF(IZ.LT.JZL) RETURN
GNOEX=1.E6*GFNOT/GFLOT
GNOVX=1.E6*GFNOVT/GFLOT
IF(ISWEEP.LT.LSWEEP.AND..NOT.ENUFSW) RETURN
WRITE(LUPR1,*)
WRITE(LUPR1,*)' *******************************'
WRITE(LUPR1,*)' * NOX FORMATION : OUTLET DATA *'
WRITE(LUPR1,*)' *******************************'
WRITE(LUPR1,1969) GFLOT,GFNOT,GNOEX,GNOVX,GEXC,GEXH,
1 GEXO,GEXN
CALL WRITBL
1969 FORMAT(/6X,' Exit mass flow rate = ',1PE12.3,' kg/s',
1 /6X,' ===================',
1 /6X,' Exit NO mass flow rate = ',1PE12.3,' kg/s'
1 /6X,' ======================',
1 /6X,' Exit NO mass fraction = ',1PE12.3,' ppm by mass'
1 /6X,' ===================== ',
1 /6X,' Exit NO vol fraction = ',1PE12.3,' ppm by vol',
1 /6X,' ==================== ',
1 /6X,' Exit C mass flow rate = ',1PE12.3,' kg/s'
1 /6X,' =====================',
1 /6X,' Exit H mass flow rate = ',1PE12.3,' kg/s'
1 /6X,' =====================',
1 /6X,' Exit O mass flow rate = ',1PE12.3,' kg/s'
1 /6X,' =====================',
1 /6X,' Exit N mass flow rate = ',1PE12.3,' kg/s'
1 /6X,' =====================',//)
ELSEIF(LSG63) THEN
C
IF(IZ.EQ.NZ) THEN
IF(ISWEEP.EQ.LSWEEP.OR.ENUFSW) THEN
C
CALL SUB2(JFU,LBNAME('FUEL'),JMF,LBNAME('MIXF'))
CALL SUB2(L0P1,L0F(P1),L0XN,L0F(JXN))
CALL SUB2(L0XNO,L0F(JXNO),L0FU,L0F(JFU))
L0MF=L0F(JMF)
C
GCO = FIXP
GVAL= 0.0
CALL SUB3R(GFLOT,0.0,GFLOF,0.0,GFLOFU,0.0)
CALL SUB2R(GFLON,0.0,GFLONO,0.0)
DO 1965 J = 1 , NX*NY
GFLOW = -GCO*(GVAL-F(L0P1+J))
GFLOT = GFLOT + GFLOW
GFLOF = GFLOF + GFLOW*F(L0MF+J)
GFLOFU = GFLOFU + GFLOW*F(L0FU+J)
GFLON = GFLON + GFLOW*F(L0XN+J)
GFLONO = GFLONO + GFLOW*F(L0XNO+J)
1965 CONTINUE
GFEX = (GFLOF /(GFLOT+1.E-10))
GFUEX = (GFLOFU/(GFLOT+1.E-10))
GNEX = (GFLON /(GFLOT+1.E-10))*1.E6
GNOEX = (GFLONO/(GFLOT+1E-10))*1.E6
CALL WRITST
WRITE(LUPR1,*)'* NOX FORMATION : OUTLET DATA *'
WRITE(LUPR1,*)'*******************************'
WRITE(LUPR1,*) ' Isweep = ',ISWEEP
WRITE(LUPR1,1966) GFLOT,GFEX,GFUEX,GNEX,GNOEX
1966 FORMAT(/6X,' Exit mass flow rate = ',
1 1PE12.4,' kg/s',
1 /6X,' ===================',/,
1 /6X,' Exit mixture fraction = ',1PE12.4,
1 /6X,' =====================',/,
1 /6X,' Exit fuel mass fraction = ',1PE12.4,
1 /6X,' =======================',/,
1 /6X,' Exit N mass fraction = ',1PE12.4,' ppm',
1 /6X,' =======================',/,
1 /6X,' Exit NO mass fraction = ',1PE12.4,' ppm',
1 /6X,' ======================',//)
C
ENDIF
ENDIF
ENDIF
ENDIF
NAMSUB='gxnox'
END
C----------------------------------------------------------------------
c
SUBROUTINE GXCREK
C***********************************************************************
C.... This subroutine is called only from GXNOX in PHOENICS-3.2; but it
C may be called from other subroutines if the user provides them.
C The PHOENICS Encyclopaedia provides information about the function
C of the CREK package generally.
C***********************************************************************
C
COMMON/CEQUIL/AL(7,20),ATOM(2,7),B0(7),PI(7)
COMMON/NCEQUL/CATOM(7)
CHARACTER*1 CATOM
COMMON/CINDEX/IHCPS,IMAT,ITER,JJ,N1,N2,N3,NA,NGLOB,NLM,NSM,NQ,
1IDCO,IDCO2,IDH2,IDH2O,IDN2,IDO2
COMMON/CKINET/CPSUM,ER,PPLN,FQ,RGAS,RGASIN,SMINV,TKINV,TLN
COMMON/CPARAM/EMV,HSUB0,NS,PA,Q0,Q1,Q2,Q3,Q4,RHOP,SM,
1S1(20),S2(20),TK
LOGICAL LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
COMMON/CHMLGC/LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
COMMON/CASUB/ASUB(20,3)
CHARACTER*4 ASUB
COMMON/CSPECE/H0(20),S0(20),SMW(20),SSAVE(20),Z(7,2,20)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER*1 CARB,HYDR
LOGICAL EQZ,QEQ
DATA FACTOR/5.0/,PATM/101325.0/,SMALL/1.0E-6/
DATA CARB/'C'/,HYDR/'H'/
C
C NORMAL SOLUTION
C
C.... Determine equivalence ratio; if outside interval (.1,10), assume
C no reaction and 230694 adiabatic non-reacted mixture properties;
C save given estimates or program generated estimates if tk is small.
C if solution is successful, 230694 to calling program. otherwise,
C enter problem cell logic below
C
c CALL WRITST
c CALL WRITBL
c CALL WRIT40('Entering subroutine GXCREK, from GXNOX ')
c CALL WRITBL
c CALL WRITST
NAMSUB='GXCREK'
CALL ERATIO
IF(.NOT.LREACT) GO TO 400
C.... Allow combustion for all mixtures....
EMVSAV=EMV
PPLN=ALOG(PA/PATM)
TSAVE=TK
DO 10 I=1,NS
10 SSAVE(I)=S2(I)
IF(TK.LT.SMALL) GO TO 100
C
CALL SPECE
C
IF(LCONVG) GO TO 900
C
C PROBLEM CELL
C
C.... Solution logic is different for four types of problems as follows:
C MODE 1 ... LEQUIL = T, LADIAB = T
C MODE 2 ... LEQUIL = T, LADIAB = F
C MODE 3 ... LEQUIL = F, LADIAB = T
C MODE 4 ... LEQUIL = F, LADIAB = F
C
C.... Always try rt=0.0 (lnrg=f) after solution failure when lequil=f.
C logic to find solution is controlled in chapters 1 and 2 below
C where in each section, new estimates are determined either by
C saved given ones, new assigned ones, or solution found not at
C required conditions. The variables nextok and nextng are assigned
C the statement numbers of where to go if the solution attempt is
C successful or not, respectively.
C
ASSIGN 900 TO NEXTOK
ASSIGN 100 TO NEXTNG
GO TO 520
C
100 MODE=4
IF(LEQUIL) MODE=MODE-2
IF(LADIAB) MODE=MODE-1
C
C* * * * * * * * * * * * * * * * * * * * * * * * CHAPTER 1 * * * * * * *
C
C EQUILIBRIUM
C.... This chapter makes equilibrium estimates amd initiates strategy for
C cases in which convergence was not achieved on first call to spece.
C
LEQUIL=.TRUE.
LADIAB=.TRUE.
IF(MODE.LT.3) GO TO 130
C
C.... First use given estimates for equil soln in mode 3 and 4 problems
C
TK=TSAVE
DO 120 I=1,NS
120 S2(I)=SSAVE(I)
ASSIGN 200 TO NEXTOK
ASSIGN 130 TO NEXTNG
GO TO 500
C
C.... Complete combustion estimate
C
130 IF(IDCO2.EQ.0) GO TO 170
C.... If no carbon in fuel, skip to garbage estimates
X=0.0
Y=0.0
DO 132 I=1,NS
S2(I)=SMALL
DO 131 L=1,NLM
IF(EQZ(AL(L,I))) GO TO 131
IF(CATOM(L).EQ.CARB) X=X+AL(L,I)*S1(I)
IF(CATOM(L).EQ.HYDR) Y=Y+AL(L,I)*S1(I)
131 CONTINUE
132 CONTINUE
S2(IDN2)=S1(IDN2)
C
ER1=(4.0*X+Y)/(2.0*X+Y)
ER2=2.0+Y/(2.0*X+1.e-20)
IF(ER.LE.1.0) GO TO 133
IF(ER.LT.ER1) GO TO 134
IF(ER.LT.ER2) GO TO 135
C.... ER.gt.ER2 ---> all c and h in co, h2 and unburned cxhy
S2(IDCO)=2.0*(X+Y/4.0)
S2(IDH2)=Y*(1.0+Y/(4.0*X))
C.... Mole number for cxhy should be (er-2*(1+y/(4*x)))
GO TO 136
C.... Fuel-lean combustion ---> only co2 and h2o formed
133 S2(IDCO2)=X*ER
S2(IDH2O)=Y*ER/2.0
S2(IDO2)=(1.0-ER)*(X+Y/4.0)
GO TO 136
C.... Sightly fuel-rich ---> co,co2, and h2o present
134 S2(IDCO)=2.0*(X+Y/4.0)*(ER-1.0)
S2(IDCO2)=Y*(1.0-ER)/2.0+X*(2.0-ER)
S2(IDH2O)=Y*ER/2.0
GO TO 136
C.... Fuel-rich ---> co,h2, and h2o present
135 S2(IDCO)=ER*X
S2(IDH2)=Y*(ER-1.0)/2.0+X*(ER-2.0)
S2(IDH2O)=Y/2.0+X*(2.0-ER)
136 SM=0.0
DO 137 I=1,NS
137 SM=SM+S2(I)
TK=1500.0
TKINV=6.6666666E-4
IHCPS=1
XLO=TK
DO 139 K=1,10
CALL HCPS
HSUM=0.0
DO 138 I=1,NS
138 HSUM=HSUM+H0(I)*S2(I)
TK=TK*(1.0+(HSUB0*RGASIN*TKINV-HSUM)/CPSUM)
TKINV=1.0/TK
XHI=ABS(TK-XLO)
XLO=TK
IF(XHI.LT.1.0) GO TO 141
139 CONTINUE
WRITE(14,140) K,XHI
140 FORMAT(/10X,'MIXTURE TEMPERATURE FAILED TO CONVERGE IN ',
1 I2,' ITERATIONS TEMP DIFFERENCES =',1PE12.3/)
141 CONTINUE
C
IF(LDEBUG) WRITE(14,142) (S2(K),K=1,NS),SM,TK
142 FORMAT(/5X,'COMPLETE COMBUSTION ESTIMATE',/1P11E10.2/1P11E10.2/)
IF(MODE.EQ.1) ASSIGN 900 TO NEXTOK
IF(MODE.EQ.2) ASSIGN 300 TO NEXTOK
IF(MODE.GE.3) ASSIGN 200 TO NEXTOK
ASSIGN 170 TO NEXTNG
GO TO 500
C
C.... Garbage estimates (gordon and mcbride)
C
170 TK=3800.0
SM=0.1/FLOAT(NS)
DO 171 I=1,NS
171 S2(I)=SM
SM=0.1
IF(MODE.EQ.1) ASSIGN 900 TO NEXTOK
IF(MODE.EQ.2) ASSIGN 300 TO NEXTOK
IF(MODE.GE.3) ASSIGN 200 TO NEXTOK
ASSIGN 600 TO NEXTNG
GO TO 500
C
C** ** ** ** ** ** ** ** ** ** ** ** CHAPTER 2 ** ** **
C
C KINETIC
C.... Section for kinetic solution from adiabatic equilibrium estimates
C (mode 3 and 4 only)
C
C.... Near-equilibrium solution (kinetic with emv=1.0e-3)
C
200 LEQUIL=.FALSE.
IX=0
EMV=1.0E-3
XLO=EMV
C.... Increase minor species from equilibrium estimates
DO 201 I=1,NS
IF(S2(I).LT.SMALL) S2(I)=SMALL
201 CONTINUE
ASSIGN 230 TO NEXTOK
ASSIGN 210 TO NEXTNG
GO TO 500
C.... Failure on near-equil with emv=xlo, decrease emv by an order of
C magnitude and attempt again, iterating this way up to 12 times.
210 EMV=EMV*0.1
XLO=EMV
IX=IX+1
IF(IX.EQ.12) GO TO 610
TK=TSAVE
DO 211 I=1,NS
211 S2(I)=SSAVE(I)
ASSIGN 230 TO NEXTOK
ASSIGN 210 TO NEXTNG
GO TO 500
C.... Have near-equil solution, so first try directly to obtain
C required solution at given emv
230 EMV=EMVSAV
IF(MODE.EQ.3) ASSIGN 900 TO NEXTOK
IF(MODE.EQ.4) ASSIGN 300 TO NEXTOK
ASSIGN 250 TO NEXTNG
GO TO 500
C
C UPPER BRANCH MARCHING
C.... Have a kinetic solution but at emv .lt. emvsv. start at known
C solution and increase emv by factor to move towards a soln there,
C if successful, repeat until emvsav is reached; if not successful
C start half interval searching described below
C
250 XLO=EMV
EMV=XLO*FACTOR
IF(EMV.GT.EMVSAV) EMV=EMVSAV
XHI=EMV
IX=0
TK=TSAVE
DO 251 I=1,NS
251 S2(I)=SSAVE(I)
ASSIGN 250 TO NEXTOK
ASSIGN 270 TO NEXTNG
GO TO 500
C
C HALF-INTERVAL SEARCHING
C.... Have solution at xlo but not at xhi, hence start interval
C searching by setting emv to the logarithmic average;
C if iterating more than ten times, terminate.
C
270 IX=IX+1
TK=TSAVE
DO 271 I=1,NS
271 S2(I)=SSAVE(I)
IF(IX.GT.10) GO TO 620
EMV=SQRT(XLO*XHI)
XHI=EMV
ASSIGN 250 TO NEXTOK
ASSIGN 270 TO NEXTNG
GO TO 500
C
C*** *** *** *** *** *** *** *** CHAPTER 3 *** ***
C
C NON-ADIABATIC
C.... Section for non-adiabatic solutions from adiabatic estimates
C (mode 2 and 4 only)
C Try directly to obtain non-adiabatic solution; if not successful,
C start half-interval scaling from the adiabatic solution by
C defining a scaling factor fq (0.0-1.0) to multiply the non-adiabatic
C term (q) in the energy equation in spece
C
300 LADIAB=.FALSE.
XLO=0.0
XHI=1.0
FQ=1.0
IX=0
310 ASSIGN 320 TO NEXTOK
ASSIGN 330 TO NEXTNG
GO TO 500
C
320 IF(QEQ(FQ,1.0)) GO TO 900
FQ=1.0
XLO=FQ
IX=0
GO TO 310
C
330 IX=IX+1
IF(IX.GT.10) GO TO 340
TK=TSAVE
DO 331 I=1,NS
331 S2(I)=SSAVE(I)
XHI=FQ
FQ=0.5*(XLO+XHI)
GO TO 310
C
340 FQ=1.0
GO TO 630
C
C**** **** **** **** **** **** CHAPTER 4 **** ***
C
C FAILURE EXITS
C.... Failed equil or kinetic soln or equiv ratio outside (0.1,10)
C return adiabatic, non-reacted mixture properties
C
400 SM=0.0
DO 401 I=1,NS
S2(I)=S1(I)
401 SM=SM+S2(I)
TK=1000.0
XLO=TK
IHCPS=1
TKINV=1.0E-3
DO 403 I=1,10
CALL HCPS
HSUM=0.0
DO 402 K=1,NS
HSUM=HSUM+H0(K)*S2(K)
402 CONTINUE
TK=TK*(1.0+(HSUB0*RGASIN*TKINV-HSUM)/CPSUM)
TKINV=1.0/TK
XHI=ABS(TK-XLO)
XLO=TK
IF(XHI.LT.1.0) GO TO 404
403 CONTINUE
WRITE(14,140) I,XHI
404 CONTINUE
RHOP=PA/(RGAS*TK*SM)
GO TO 900
C
C***** ***** ***** ***** ***** CHAPTER 5 ***** *
C
C PROBLEM CELL CALL TO SPECE
C.... Take the estimates generated in chapters 1,2 and attempt a solution
C with full equations. If successful, update the save answers with the
C solution and return to statement number nextok. If not, the action
C depends on whether an equilbm or kinetic soln is sought. Failed
C equil soln, return to statement number nextng, while failure in a
C kinetic soln will be followed by an attempt with lnrg=f ---> rt=0.0
C and same estimates. Setting rt=0.0 implies that a change in temp
C field has no effect on species distribution for that particular
C iteration, but does allow the species changes to influence the temp
C change ---> partial decoupling of the energy equation. If successful
C with lnrg=f, repeat with lnrg=t (full equations); if still no good,
C return to statement number nextng.
C
500 CALL SPECE
IF(LCONVG) GO TO 540
C SOLUTION FAILED TRY RT=0.0
IF(LNRG) GO TO 520
LNRG=.TRUE.
510 GO TO NEXTNG, (100,130,170,210,250,270,330,600)
520 IF(LEQUIL) GO TO 510
LNRG=.FALSE.
C**MPD 530 LNRG=.FALSE. >>>> UNREFERENCED STATEMENT LABEL <<<< ??
GO TO 500
C HAVE SOLN BUT WITHOUT RT TERMS, SEND THIS
C SOLN AS ESTIMATES WITH RT CALCULATION
540 IF(LNRG) GO TO 550
LNRG=.TRUE.
GO TO 500
C
C.... Solution is successful, update save answers and continue.
C
550 TSAVE=TK
DO 560 I=1,NS
560 SSAVE(I)=S2(I)
C
GO TO NEXTOK, (200,230,300,320,900)
C
C****** ****** ****** ****** CHAPTER 6 ******
C
C ERROR MESSAGES
C
600 CONTINUE
C WRITE(14,601)
C 601 FORMAT(/10X,3(4H****),' FAILURE TO FIND EQUIL SOLN...',
C 1 'AVG INLET PROPS RETURNED'/)
C RETURN INITIAL GUESS....NOT THE UNREACTED MIXTURE PROPERTIES
GO TO 670
610 WRITE(14,611)
611 FORMAT(/10X,3(4H****),' FAILURE TO FIND NEAR-EQUIL SOLN...',
1 'AVG INLET PROPS RETURNED'/)
GO TO 650
620 WRITE(14,621)
621 FORMAT(/10X,3(4H****),' FAILURE TO OBTAIN KINETIC SOLN AFTER',
1 ' TEN INTERVAL HALVING...AVG INLET PROPS RETURNED'/)
GO TO 650
630 WRITE(14,631)
631 FORMAT(/10X,3(4H****),' NON-ADIABATIC SOLN FAILED...',
1 'ADIAB SOLN RETURNED'/)
GO TO 670
C
C.... Restore failed problem mode prior to return
C
650 IF(MODE.EQ.2) LADIAB=.FALSE.
IF(MODE.EQ.3) LEQUIL=.FALSE.
IF(MODE.EQ.4) LADIAB=.FALSE.
GO TO 400
C
C.... Failed non-adiabatic solution...return adiabatic
c equil or kinetic solution
C
670 TK=TSAVE
DO 671 I=1,NS
671 S2(I)=SSAVE(I)
C
900 NAMSUB='gxcrek'
END
c