c
C.... File name ................. GXDENS.HTM .................... 021110
C
FUNCTION DENSITY(I)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
1 /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
COMMON/HBASE/IH01,IH02,KH01,KH01H,KH01L,KH02,KH02H,KH02L,L0H012
common/dbe/dbear
logical dbear
C
c
IF(IGRND.EQ.-1) THEN ! default when RHO1 or 2 is greater than zero
DENSITY=PRPRTY
ELSEIF(IGRND.LE.2) THEN ! when RHO1 or 2=GRND1 or GRND2
c
C.... Density is a linear function or the reciprocal of a linear function
C of either a scalar (enthalpy or other scalar); or two scalars; or
C three:
c
IF(L0SCAL.GT.0) DENSITY=RHOAG + RHOBG*F(L0SCAL+I)
IF(L0CA.GT.0) DENSITY=DENSITY + RHOCG*F(L0CA+I)
IF(L0CB.GT.0) DENSITY=DENSITY + RHODG*F(L0CB+I)
c
IF(IGRND.EQ.2) DENSITY=1./DENSITY
c
ELSEIF(IGRND.EQ.3) THEN ! when RHO1 or 2=GRND1
c
C.... Density for the compressible flow and also shallow-water theory:
cccc call writ3r('rhoag ',rhoag,'rhobg ',rhobg,'rhocg ',rhocg)
DENSITY=RHOAG * (AMAX1(0.,ABSPRES(I)))**RHOBG + RHOCG
c
C.... Multiplication by C3 is useful when two compressible gases are
C present in the same domain. C3 has of course to be computed so as
C to "travel" with the gas which it marks:
c
IF(L0SCAL.NE.0) DENSITY=DENSITY*F(L0SCAL+I)
c
C.... For shallow-water theory with NZ > 1, and compressibility in the
C upper (lower-IZ) cell only, the density is set to -RHO2 except in
C topmost cells (IZ=1):
c
IF(SHALWT) DENSITY=- RHO2
c
ELSEIF(IGRND.EQ.4) THEN
c
C.... Density is a linear function of temperature:
c
DENSITY=RHOAG + RHOBG*(F(L0TEM+I)+TEMP0)
c
ELSEIF(IGRND.EQ.5) THEN
c
C.... Density according to the ideal-gas law:
ABST=ABSTEMP(I) ! get absolute temperature
IF(L0SCAL.GT.0.AND.RHOAG.GT.0.0) THEN ! for a mixture of 2 gases
SPEVOL=RHOAG + (RHOBG-RHOAG) * F(L0SCAL+I) ! the scalar is C1 or
DENSITY=ABSPRES(I) / (SPEVOL * ABST) ! C2, according to phase
ELSE ! RHOAG and RHOBG are
DENSITY=RHOBG * ABSPRES(I) / ABST ! the ideal-gas constant
ENDIF ! of the gases
c
ELSEIF(IGRND.EQ.6) THEN
c
C.... Ideal-gas law for mixture of the three chemical 'species' fuel,
C oxidant and product involved in the simple chemical reaction
C model (Parameters RHOA,RHOB and RHOC are the molecular weights
C of fuel, oxidant and products respectively):
c
ABST=ABSTEMP(I)
IF(ABST.GT.0.0) THEN
RMIX=(F(L0FUEL+I)/RHOAG + F(L0OXID+I)/RHOBG +
1 F(L0PROD+I)/RHOCG) * GASCON
DENSITY=ABSPRES(I)/(RMIX * ABST)
ELSE
call writ40('absolute temperature less than zero')
call writ40('density set to 1.0 ')
density=1.0
ENDIF
c
ELSEIF(IGRND.EQ.7) THEN
c
C.... Ideal-gas law for a mixture of the 7 gases involved in a simple
C model of the combustion of hydrocarbons in air.
Click here for explanation
C.... Mass fraction of anything but air is FNOTAIR.
C!!!! Note that prior to the generalization of March 2002, it has been
c that any material injected which did not have a finite MIXF value
c was air, with .232 O2 content and .768 N2 content.
c This presumption is preserved here.
c
C.... FO, FC and FH are the mass fractions of the elements O, C and H,
C regardless of which compound they are contained in.
FC=0.0 ; FH=0.0 ; FO=0.0 ; FN=0.0 ! initialise to zero
FNOTAIR=0.0
IF(L0MIXF.NE.0) THEN
FMXF0=F(L0MIXF+I) ! mass fraction of material from inlet 0
FH=FH + RHOBG * FMXF0 ! stream * hydrogen content (rhobg)
FC=FC + RHOAG * FMXF0 ! carbon content (rhoag)
FN=FN + RHOCG * FMXF0 ! nitrogen content (rhocg)
! inlet 0 contains no oxygen
FNOTAIR=FNOTAIR + FMXF0
ENDIF
IF(L0MXF1.NE.0) THEN
FMXF1=F(L0MXF1+I) ! mass fraction of material from inlet 1
FH=FH + FHMXF1 * FMXF1 ! stream * hydrogen content (FH1)
FC=FC + FCMXF1 * FMXF1 ! carbon content (FC1)
FO=FO + FOMXF1 * FMXF1 ! oxygen content (FO1)
FN=FN + FNMXF1 * FMXF1 ! nitrogen content (FN1)
FNOTAIR=FNOTAIR + FMXF1
ENDIF
IF(L0MXF2.NE.0) THEN
FMXF2=F(L0MXF2+I) ! mass fraction of material from inlet 2
FH=FH + FHMXF2 * FMXF2 ! stream * hydrogen content (FH2)
FC=FC + FCMXF2 * FMXF2 ! carbon content (FC2)
FO=FO + FOMXF2 * FMXF2 ! oxygen content (FO2)
FN=FN + FNMXF2 * FMXF2 ! nitrogen content (FN2)
FNOTAIR=FNOTAIR + FMXF2
ENDIF
IF(L0MXF3.NE.0) THEN
FMXF3=F(L0MXF3+I) ! mass fraction of material from inlet 3
FH=FH + FHMXF3 * FMXF3 ! stream * hydrogen content (FH3)
FC=FC + FCMXF3 * FMXF3 ! carbon content (FC3)
FO=FO + FOMXF3 * FMXF3 ! oxygen content (FO3)
FN=FN + FNMXF3 * FMXF3 ! nitrogen content (FN3)
FNOTAIR=FNOTAIR + FMXF3
ENDIF
IF(L0MXF4.NE.0) THEN
FMXF4=F(L0MXF4+I) ! mass fraction of material from inlet 4
FH=FH + FHMXF4 * FMXF4 ! stream * hydrogen content (FH4)
FC=FC + FCMXF4 * FMXF4 ! carbon content (FC4)
FO=FO + FOMXF4 * FMXF4 ! oxygen content (FO4)
FN=FN + FNMXF4 * FMXF4 ! nitrogen content (FN4)
FNOTAIR=FNOTAIR + FMXF4
ENDIF
IF(L0MXF5.NE.0) THEN
FMXF5=F(L0MXF5+I) ! mass fraction of material from inlet 5
FH=FH + FHMXF5 * FMXF5 ! stream * hydrogen content (FH5)
FC=FC + FCMXF5 * FMXF5 ! carbon content (FC5)
FO=FO + FOMXF5 * FMXF5 ! oxygen content (FO5)
FN=FN + FNMXF5 * FMXF5 ! nitrogen content (FN5)
FNOTAIR=FNOTAIR + FMXF5
ENDIF
IF(L0MXF6.NE.0) THEN
FMXF6=F(L0MXF6+I) ! mass fraction of material from inlet 6
FH=FH + FHMXF6 * FMXF6 ! stream * hydrogen content (FH6)
FC=FC + FCMXF6 * FMXF6 ! carbon content (FC6)
FO=FO + FOMXF6 * FMXF6 ! oxygen content (FO6)
FN=FN + FNMXF6 * FMXF6 ! nitrogen content (FN6)
FNOTAIR=FNOTAIR + FMXF6
ENDIF
IF(L0MXF7.NE.0) THEN
FMXF7=F(L0MXF7+I) ! mass fraction of material from inlet 7
FH=FH + FHMXF7 * FMXF7 ! stream * hydrogen content (FH7)
FC=FC + FCMXF7 * FMXF7 ! carbon content (FC7)
FO=FO + FOMXF7 * FMXF7 ! oxygen content (FO7)
FN=FN + FNMXF7 * FMXF7 ! nitrogen content (FN7)
FNOTAIR=FNOTAIR + FMXF7
ENDIF
IF(L0FO.NE.0) F(L0FO+I)=FO
IF(L0FC.NE.0) F(L0FC+I)=FC
IF(L0FH.NE.0) F(L0FH+I)=FH
IF(L0FN.NE.0) F(L0FN+I)=FN
IF(WOOD) THEN ! -------------------'wood' section
C.... Mass fraction of wood derivatives, FWD:
FWD=F(L0FWD+I) ! FWD is treated as a source-free variable
FH=FH + HINWD*FWD ! its contribution to h-content of 'gas'
FC=FC + CINWD*FWD ! and to c
FO=FO + OINWD*FWD ! and to o
FNOTAIR=FNOTAIR + FWD
F(L0CHAR+I)=AMIN1(F(L0CHAR+I), FWD*CHINWD) ! char concentration is
YCHAR=F(L0CHAR+I) ! computed elsewhere
YWOOD=F(L0WOOD+I) ! wood also
ELSE ! when 'wood is not present
C.... Coal- (or oil-) burning in gaseous phase (i.e. FGAS=1):
FGAS=1.0
YWOOD=0.0
YCHAR=0.0
FWD=0.0
ENDIF ! end of 'wood'section
C.... Compute the species concentrations (N2 can be in coal & air):
FAI=1.0 - FNOTAIR
FO=FO + 0.232 * FAI
FN=FN + 0.768 * FAI
YN2=FN ! initial guess
IF(WOOD) THEN
C.... Mass fraction of phase 1 in the gaseous state:
C Note that it is necessary to recognise that wood and char may be
C found anywhere because their burning is kinetically controlled;
FGAS=1.- YWOOD - YCHAR ! proportion of phase 1 which is true gas
RFGAS=1./FGAS ! its reciprocal
C.... Mass fractions of elements in the gas
GO=(FO - OINWD*YWOOD)*RFGAS
GC=(FC - CINWD*YWOOD - YCHAR)*RFGAS
GH=(FH - HINWD*YWOOD)*RFGAS
ELSE
GO =FO
GC =FC
GH =FH
ENDIF
C.... Determine composition of gas on the presumption that mixed is
C burned:
ZCO2=3.666667*GC ! 3.66667=44 / 12=mlwt CO2 / mlwt C
ZH2O=9.0*GH ! 9.0 =18 / 2=mlwt H2O / mlwt H
ZO2=GO - 2.666667*GC - 8.0*GH ! 2.6667=32 / 12; 8.0=16 / 2
YVOL=0.0 ! OK ?????
IF(ZO2.GE.0.0) THEN ! region 1: O2, H2O, CO2, N2are present
ZH2=0.0
YH2=0.0
ZCO=0.0
YCO=0.0
ZVOL=0.0
YO2=ZO2*FGAS
YCO2=ZCO2*FGAS
YH2O=ZH2O*FGAS
IF(WOOD) THEN
YVOL=YVOL + VLINWD*YWOOD ! ??? how can volatiles be finite on R 1
YH2O=YH2O + MOINWD*YWOOD
ENDIF
YN2=1.0 - YO2 - YCO2 - YH2O ! ?????
yn2=amax1(yn2,0.0)
RMIX=GASCON*(YO2*0.03125 + YH2O*0.055556 + YCO2*0.022727
1 + YN2*0.035714 + TINY)
HSUB=0.0
ELSE
C.... O2 is not present, i.e. it is region 2 or 3
ZO2=0.0
YO2=0.0
C.... Suppose free CO, CO2, H2, H2O exist at the full-oxidation limit,
C i.e. GOfull=((16/2)*GH + (32/12)*GC)*(1. - GOfull)/(1. - GO)
GOFULL=(8.*GH + 2.666667*GC)/(1.- GO + 8.*GH + 2.666667*GC)
FACTFU=(1.- GOFULL)/(1.- GO)
ZCO2FU=FACTFU*3.666667*GC
ZH2OFU=FACTFU*9.0*GH
C.... At the part-oxidation limit only CO and H2 exist:
GOPART=1.333333*GC/(1.- GO + 1.333333*GC)
FACTPA=(1.- GOPART)/(1.- GO)
ZCOPA=FACTPA*2.333333*GC
ZH2PA=FACTPA*GH
C.... For the actual mixture:
FRAC=(GO - GOPART)/(GOFULL - GOPART)
IF(FRAC.GE.0.0) THEN
C.... The composition of region 2 (i.e. CO2,H2O,CO,H2 and N2):
ZVOL=0.0
YVOL=0.0
ZCO2=ZCO2FU*FRAC
ZH2O=ZH2OFU*FRAC
ZCO=ZCOPA *(1.- FRAC)
ZH2=ZH2PA *(1.- FRAC)
YCO2=ZCO2*FGAS
YH2O=ZH2O*FGAS
YCO=ZCO*FGAS
YH2=ZH2*FGAS
IF(WOOD) THEN
YVOL=YVOL + VLINWD*YWOOD
YH2O=YH2O + MOINWD*YWOOD
ENDIF
YN2=1.0 - YCO2 - YCO - YH2O - YH2 ! ?????
yn2=amax1(yn2,0.0)
RMIX=GASCON*(YCO2*0.022727 + YH2O*0.055556 + YCO*0.035714
1 + YH2*0.5 + YN2*0.035714 + TINY)
HSUB=YCO*HCOCO2 + YH2*HHH2O
ELSE
C.... The composition of region 3 (i.e. CO,H2,N2 and volatiles):
ZCO2=0.0
YCO2=0.0
ZH2O=0.0
YH2O=0.0
ZCO=1.75*GO
IF(WOOD) THEN
CONSTA=CINVL
CONSTB=HINVL
ELSE
CONSTA=RHOAG
CONSTB=RHOBG
ENDIF
ZVOL=(GC - ZCO*0.428571)/(CONSTA+tiny)
ZH2=GH - ZVOL*CONSTB
YCO=ZCO*FGAS
YH2=ZH2*FGAS
YVOL=ZVOL*FGAS
IF(WOOD) THEN
YVOL=YVOL + VLINWD*YWOOD
YH2O=YH2O + MOINWD*YWOOD
ENDIF
YN2=1.0 - YCO - YH2O - YH2 - YVOL ! ?????
yn2=amax1(yn2,0.0)
RMIX=GASCON*(YCO*0.035714 + YH2*0.5 + YN2*0.035714
1 + YVOL*RMWTVL + TINY)
HSUB=YCO*HCOCO2 + YH2*HHH2O + YVOL*HVOL
ENDIF
ENDIF
C.... Calculate the absolute temperatures and the first-phase density:
IF(WOOD) HSUB=HSUB + YCHAR*HCHAR + YWOOD*HWOOD
C... Save HSUB into H0 store
F(L0H012+I)=HSUB
TEM=(F(L0H1+I) - HSUB)/F(L0CP1+I)
TEM=AMIN1(VARMAX(TEMP1), AMAX1(100., TEM, VARMIN(TEMP1)))
DENS=ABSPREs(I)/(RMIX*TEM)
C.... Wood and char are supposed to contribute their mass to the mixture,
C without significantly increasing its volume.
IF(WOOD) DENS=DENS + YWOOD + YCHAR
DENSITY=AMIN1(VARMAX(DEN1), AMAX1(0., DENS, VARMIN(DEN1)))
IF(.NOT.ONEPHS) THEN
TEM2=(F(L0H2+I) - HCHX)/F(L0CP2+I)
C... save HCHX into H02 store
F(KH02+I)=HCHX
TEM2=AMIN1(VARMAX(TEMP2), AMAX1(100., TEM2, VARMIN(TEMP2)))
IF(L0TMP2.NE.0) F(L0TMP2+I)=TEM2
C.... Estimated interface temperature:
CINT1=CINT(H1)
CINT2=CINT(H2)
FACTEM=(TEM*CINT1 + TEM2*CINT2)/(CINT1 + CINT2)
PHINT1=F(L0CP1+I)*FACTEM + HSUB
PHINT2=F(L0CP2+I)*FACTEM + HCHX
C.... Interface values used for conductive transport between phases
F(L0H1IF+I)=PHINT1
F(L0H2IF+I)=PHINT2
ENDIF
IF(LBYCO. NE.0) F(L0YCO+I)=YCO
IF(LBYCO2.NE.0) F(L0YCO2+I)=YCO2
IF(LBYH2. NE.0) F(L0YH2+I)=YH2
IF(LBYH2O.NE.0) F(L0YH2O+I)=YH2O
IF(LBYO2. NE.0) F(L0YO2+I)=YO2
IF(LBRMIX.NE.0) F(L0RMIX+I)=RMIX
IF(LBYN2. NE.0) F(L0YN2+I)=YN2
IF(LBYVOL.NE.0) F(L0YVOL+I)=YVOL
IF(l0tem. NE.0) F(L0TEM+I)=TEM
IF(l0tmp2. NE.0) F(L0TMP2+I)=TEM2
c
IF(L0HSUB.NE.0) F(L0HSUB+I)=HSUB
c-------------------------------------------- end of 7-gases option
ELSEIF(IGRND.EQ.8) THEN ! Density for cvd special-purpose program
CALL CVDPRP(DENSITY,IFILEP,I,2)
c
ELSEIF(IGRND.EQ.9) THEN ! Density from CHEMKIN-database
CALL GXCHKI(I,DENSITY)
c
ELSEIF(IGRND.EQ.10) THEN ! Density=1/linear function of temperature
IF(L0TEM.GT.0) THEN
DENSITY=1.0/(RHOAG + RHOBG*(F(L0TEM+I)+TEMP0))
ELSE
call writ40('temperature must be solved for or stored')
call writ1i('igrnd ',igrnd)
CALL SET_ERR(215,'Error. Temperature must be solved or stored.'
* ,1)
C call earout(1)
endif
ENDIF
IF(IGRND.NE.-1) THEN
LBDEN=ITWO(DEN1,DEN2,IPHASE.EQ.1)
IF(LBDEN.GT.0) THEN
DENSITY=AMAX1(VARMIN(LBDEN),AMIN1(VARMAX(LBDEN),DENSITY))
IF(DENSITY.LE.0.0) THEN
call writ1r('density ',density)
call writ2i('den1 ',den1,'den2 ',den2)
call writ1r('density ',density)
call writ1r('varmin ',varmin(lbden))
call writ1r('varmax ',varmax(lbden))
call writ40('density.le.0.0, in gxdens; run aborted')
call writ3i('igrnd ',igrnd,'cell no.',i,'izstep ',izstep)
call writ3i('isweep ',isweep,'istep ',istep,
1 'indvar ',indvar)
CALL SET_ERR(216,'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ENDIF
ENDIF
END
c-----------------------------------------------------------------------
C
SUBROUTINE SLBDEN(IPILOPT,USEPIL,dbgloc)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(21),IPRL,
1 IBTAU,IGF4(16),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
1 IPRPS,IRADX,IRADY,IRADZ,IVFOL
1 /F0/ IF01(29),L0XYDX,L0XYDY,IF02(3),L0XYRV,L0XYXG,IF03,
1 L0XYYG,IF04,L0XYDXG,L0XYDYG,IF05(106),L0AHZ,IF06(17),
1 L0XYDZ,L0XYDZG,IF07(137)
COMMON/NAMFN/NAMFUN,NAMSUB
cccc common/denonly/igrnden
LOGICAL USEPIL,DBGLOC,SLD
common/dbe/dbear
logical dbear
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB='SLBDEN'
if(dbgloc) call banner(1,namsub,220306)
IGR=9
ISC=IPROP
if(dbgloc) call writ1i('ipilopt ',ipilopt)
C.... Call GROUND for the user-set property:
IF(IPILOPT.EQ.0) THEN
IF(USEGRD) THEN
if(dbgloc) then
call writ40('GROUND is called to set a property... ')
call writ2i('Group= ',igr,'Section=',isc)
endif
CALL GROUND
ENDIF
GO TO 800
C.... For a constant property go to the slab loop
ELSEIF(IPILOPT.EQ.-1.AND.USEPIL) THEN
if(dbgloc) call writ40('property is uniform')
GO TO 700
ENDIF
c
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
C.... Density of the first or second phase:
IF(IPHASE.EQ.1) THEN
RHOAG=RHO1A
RHOBG=RHO1B
RHOCG=RHO1C
ELSE
RHOAG=RHO2A
RHOBG=RHO2B
RHOCG=RHO2C
ENDIF
C.... Density using Extended Simple Chemically-Reacting System:
IF(IPILOPT.EQ.9 .AND. GRNDNO(9,RHOAG)) THEN
CALL GXSCRI
GO TO 800
ENDIF
c.... make necessary initial settings
CALL SETDEN
c
C.... Loop over slab to get and set cell properties:
700 IF(IPRPS.EQ.0) THEN ! PRPS is not stored; one material only, no blockages
IGRND=IPILOPT
if(dbrho) then
write(14,*) 'setting phase-',iphase,
1 ' density via PIL with igrnd=',igrnd
endif
DO 60 I=1,NXNYST
60 F(KPROP+I)=DENSITY(I) ! set slab values via PIL values
ELSEIF(USEPIL) THEN
c.... One material only, with blockages
IGRND=IPILOPT
DO 70 I=1,NXNYST
IF(SLD(I)) THEN
F(KPROP+I)=TINY
ELSE
F(KPROP+I)=DENSITY(I)
ENDIF
70 CONTINUE
ELSEIF(SURF .OR. HOL) THEN
C.... Properties are set for either scalar equation method (SURF=T), or
C height-of-liquid method (HOL=T):
CALL SURHOL(KPROP,IPROP,IFILEP,IPHASE.EQ.1)
ELSE
if(dbgloc) call writ40('prfile called; use imat for each cell')
C.... Properties are set using material flag for an individual cell.
CALL PRFILE(IPILOPT,dbgloc)
ENDIF
C----------------------------------------------------------------------
C.... Corrections, debug print-out, and other property adjustments:
if(dbgloc .and. IPILOPT.eq.7) then
if(izstep.ge.izdb1 .and.izstep.le.izdb2.and.
1 isweep.ge.iswdb1.and.isweep.le.iswdb2 ) then
if(l0mdot.ne.0) call prn('mdt*',-l0mdot)
if(l0fo.ne.0) call prn('fo ',-l0fo)
if(l0fc.ne.0) call prn('fc ',-l0fc)
if(l0fh.ne.0) call prn('fh ',-l0fh)
if(l0hsub.ne.0) call prn('hsub',-l0hsub)
if(wood) then
call prn('wood',-l0wood)
call prn('char',-l0char)
call prn('vola',-l0vola)
endif
if(lbyco. ne.0) call prn('yco ',lbyco)
if(lbyco2.ne.0) call prn('yco2',lbyco2)
if(lbyh2. ne.0) call prn('yh2 ',lbyh2)
if(lbyh2o.ne.0) call prn('yh2o',lbyh2o)
if(lbyo2. ne.0) call prn('yo2 ',lbyo2)
if(lbrmix.ne.0) call prn('rmix',lbrmix)
if(lbyn2. ne.0) call prn('yn2 ',lbyn2)
if(lbyvol.ne.0) call prn('yvol',lbyvol)
if(temp1. ne.0) call prn('tem ',temp1)
call prn('h1 ',-l0h1)
if(.not.onephs) then
call prn('tem2',temp2)
call prn('h1if',-l0h1if)
call prn('h2 ',-l0h2)
call prn('h2if',-l0h2if)
endif
endif
ENDIF
C.... Call GREX to correct a property set above:
800 IF(USEGRX) CALL GREX3
IF(USEALT) CALL ALTPRP
IF(USEGRD) THEN
IF(IPILOPT.GT.0) CALL GROUND
ENDIF
if(dbgloc) call banner(2,'slbden',0)
NAMSUB='slbden'
END
c-----------------------------------------------------------------------
C
SUBROUTINE SETDEN
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,NWHOLE,IGSH,
1 IGF3(37),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
COMMON /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
LOGICAL QEQ,LTZ
c-----------------------------------------------------------------------
IF(SOLVE(1)) L0P=L0F(1)
c???? can the next be simplified?
IF(IPHASE.EQ.1) THEN
IF(ITEM1.NE.0) THEN
L0TEM=L0F(ITEM1)
ELSEIF(TEMP1.GT.0) THEN
L0TEM=L0F(TEMP1)
ELSEIF(TMP1.LE.GRND.OR.QEQ(TMP1,1.11111111)) THEN
L0TEM=L0F(TEMP1)
ENDIF
L0DEN=L0F(DEN1)
ELSE
IF(ITEM2.NE.0) THEN
L0TEM=L0F(ITEM2)
ELSEIF(TMP2.LE.GRND) THEN
L0TEM=L0F(TEMP2)
ENDIF
L0DEN=L0F(DEN2)
IF(GEN2.NE.0) L0GEN=L0F(GEN2)
ENDIF
L0SCAL=0 ; L0CA=0 ; L0CB=0
IF(IGRND.EQ.1 .OR. IGRND.EQ.2) THEN
C.... Density is a function of either one variable (enthalpy or other
C scalar); or two, or three scalars:
IF(IBUOYB.GT.0) THEN
IPH=IPHASE - 1
L0SCAL=L0F(IBUOYB+IPH)
IF(IBUOYC.GT.0) L0CA=L0F(IBUOYC+IPH)
IF(IBUOYA.GT.0) THEN
L0CB=L0F(IBUOYA+IPH)
RHODG=RHO2A
ELSE
L0CB=L0CA
RHODG=0.0
ENDIF
ELSE
LBH =ITWO(H1,H2,IPHASE.EQ.1)
LBSCAL=ITWO(LBH , NINT(RHOCG) , QEQ(RHOCG,0.0))
L0SCAL=L0F(LBSCAL)
ENDIF
ELSEIF(IGRND.EQ.3) THEN
C.... Density for the compressible flow and also shallow-water theory:
SHALWT=LTZ(RHO2) .AND. ONEPHS .AND. IZSTEP.GT.1
IF(STORE(C3)) L0SCAL=L0F(C3)
ELSEIF(IGRND.EQ.5) THEN
IF(STORE(C1)) L0SCAL=L0F(C1 + ITWO(0,1,IPHASE.EQ.1))
ELSEIF(IGRND.EQ.6) THEN
C.... Ideal-gas law for mixture of three chemical 'species':
L0FUEL=L0F(LBFUEL) ; L0OXID=L0F(LBOXID) ; L0PROD=L0F(LBPROD)
ELSEIF(IGRND.EQ.7) THEN
C.... Ideal-gas law for a mixture of the 7 gases involved in a simple
C model of the combustion in air of one or more hydrocarbons
C containing C, H, N and an inert ash (coal, oil or wood). The gases
C are: HCx, O2, N2, CO2, CO, H2O, and H2.
c
c allowance is made for material entering from up to 8 supply points,
c the local mass fractions of which being measures by the variables
c MIXF, MXF1, etc up to MXF7
c
c Each stream can have its own elemental mass fraction FC, FH, FO, FN
c which must be read from the Q1 file (except that the values for the
c MIXF stream are supplied by way of ????)
IF(LBMIXF.GT.0) L0MIXF=L0F(LBMIXF)
c
IF(LBMXF1.GT.0) L0MXF1=L0F(LBMXF1)
IF(LBMXF2.GT.0) L0MXF2=L0F(LBMXF2)
IF(LBMXF3.GT.0) L0MXF3=L0F(LBMXF3)
IF(LBMXF4.GT.0) L0MXF4=L0F(LBMXF4)
IF(LBMXF5.GT.0) L0MXF5=L0F(LBMXF5)
IF(LBMXF6.GT.0) L0MXF6=L0F(LBMXF6)
IF(LBMXF7.GT.0) L0MXF7=L0F(LBMXF7)
c
IF(LBFC.GT.0) L0FC= L0F(LBFC)
IF(LBFH.GT.0) L0FH= L0F(LBFH)
IF(LBFO.GT.0) L0FO= L0F(LBFO)
IF(LBFN.GT.0) L0FN= L0F(LBFN)
c
L0H1 =L0F(H1)
L0CP1=L0F(ISPH1)
IF(TEMP1.NE.0) L0TEM=L0F(TEMP1)
IF(.NOT.ONEPHS) THEN
L0MDOT=L0F(INTMDT)
L0H2=L0F(H2)
L0CP2=L0F(ISPH2)
IF(TEMP2.NE.0) L0TMP2=L0F(TEMP2)
ENDIF
IF(LBYCO.NE.0) L0YCO=L0F(LBYCO)
IF(LBYCO2.NE.0) L0YCO2=L0F(LBYCO2)
IF(LBYH2.NE.0) L0YH2=L0F(LBYH2)
IF(LBYH2O.NE.0) L0YH2O=L0F(LBYH2O)
IF(LBYO2.NE.0) L0YO2=L0F(LBYO2)
IF(LBYN2.NE.0) L0YN2=L0F(LBYN2)
IF(LBRMIX.NE.0) L0RMIX=L0F(LBRMIX)
IF(LBYVOL.NE.0) L0YVOL=L0F(LBYVOL)
IF(WOOD) THEN
L0WOOD=L0F(LBWOOD)
c moved to start LBFWD=LBNAME('FWD ')
L0FWD=L0F(LBFWD)
L0CHAR=L0F(LBCHAR)
ENDIF
ENDIF
cccc if(flag) call banner(2,'setden',220306)
END
c-----------------------------------------------------------------------
C
c.... INIDEN is called from INIPRP at the beginning of a run to make
c once for all settings related to density
c
SUBROUTINE INIDEN
INCLUDE 'satear'
INCLUDE 'satgrd'
INCLUDE 'prpcmn'
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
1 /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
LOGICAL GRNDNO,EQZ,GRPROP(9),GRN,STORE
DIMENSION NAMPRP(30)
COMMON/HBASE/IH01,IH02,KH01,KH01H,KH01L,KH02,KH02H,KH02L,L0H012
EQUIVALENCE (NAMPRP(1),DENST1)
C
if(dbrho) call banner(1,'INIDEN',220306)
C.... Check what options will be used for first-phase density, and store
c the result in the grprop array:
c only, If that something is needed, I will find a less
c obscure way of doing it. I shall de-acriate the sybroutine.
c
cc CALL PRPGRN(RHO1,1,GRPROP) ! fill grprop array. Only item 7 is used below
IF(.NOT.STORE(1)) THEN
IF(GRNDNO(3,RHO1).OR.GRNDNO(5,RHO1).OR.GRNDNO(6,RHO1).OR.
1 GRNDNO(7,RHO1).OR.GRNDNO(3,RHO2).OR.GRNDNO(5,RHO2).OR.
1 GRNDNO(6,RHO2).OR.GRNDNO(7,RHO2) ) THEN
CALL WRIT40('RHOx=GRND3,5,6,7 requires pressure... ')
CALL SET_ERR(217,'Error. RHOx=GRND3,5,6,7 requires pressure...'
* ,2)
endif
endif
C.... Default settings for compressible flow.
IF(GRNDNO(3,RHO1)) THEN
DRH1DP=RHO1
ELSEIF(GRNDNO(3,RHO2)) THEN
DRH2DP=RHO2
ELSEIF(GRNDNO(4,RHO1).OR.GRNDNO(5,RHO1).OR.GRNDNO(6,RHO1)) THEN
IF(ITEM1.EQ.0 .AND. TEMP1.LT.0 .AND. .NOT.GRN(TMP1)) THEN
CALL WRIT40('RHO1=GRND4,5,6 requires temperature... ')
CALL SET_ERR(218,'Error. RHO1=GRND4,5,6 requires temperature...'
* ,2)
C CALL EAROUT(2)
ENDIF
IF(GRNDNO(5,RHO1)) DRH1DP=RHO1
ELSEIF(GRNDNO(4,RHO2).OR.GRNDNO(5,RHO2).OR.GRNDNO(6,RHO2)) THEN
IF(ITEM2.EQ.0 .AND. TEMP2.LT.0 .AND. .NOT.GRN(TMP2)) THEN
CALL WRYT40('RHO2=GRND4,5,6 requires temperature... ')
CALL SET_ERR(219,'Error. RHO2=GRND4,5,6 requires temperature...'
* ,2)
ENDIF
IF(GRNDNO(5,RHO2)) DRH2DP=RHO2
c 7gases
c ELSEIF(GRPROP(7) .OR. GRNDNO(7,RHO2)) THEN !----------- 7 gases model
ELSEIF(GRNDNO(7,RHO1).OR. GRNDNO(7,RHO2)) THEN !----------- 7 GASES
IF(LBMXF1.NE.0) THEN
CALL RQ1R('7GAS','FH1',FHMXF1)
CALL RQ1R('7GAS','FC1',FCMXF1)
CALL RQ1R('7GAS','FO1',FOMXF1)
CALL RQ1R('7GAS','FN1',FNMXF1)
CALL RQ1R('7GAS','HF1',HMXF1)
ENDIF
IF(LBMXF2.NE.0) THEN
CALL RQ1R('7GAS','FH2',FHMXF2)
CALL RQ1R('7GAS','FC2',FCMXF2)
CALL RQ1R('7GAS','FO2',FOMXF2)
CALL RQ1R('7GAS','FN2',FNMXF2)
CALL RQ1R('7GAS','HF1',HMXF2)
ENDIF
IF(LBMXF3.NE.0) THEN
CALL RQ1R('7GAS','FH3',FHMXF3)
CALL RQ1R('7GAS','FC3',FCMXF3)
CALL RQ1R('7GAS','FO3',FOMXF3)
CALL RQ1R('7GAS','FN3',FNMXF3)
CALL RQ1R('7GAS','HF1',HMXF3)
ENDIF
IF(LBMXF4.NE.0) THEN
CALL RQ1R('7GAS','FH4',FHMXF4)
CALL RQ1R('7GAS','FC4',FCMXF4)
CALL RQ1R('7GAS','FO4',FOMXF4)
CALL RQ1R('7GAS','FN4',FNMXF4)
CALL RQ1R('7GAS','HF1',HMXF4)
ENDIF
IF(LBMXF5.NE.0) THEN
CALL RQ1R('7GAS','FH5',FHMXF5)
CALL RQ1R('7GAS','FC5',FCMXF5)
CALL RQ1R('7GAS','FO5',FOMXF5)
CALL RQ1R('7GAS','FN5',FNMXF5)
call writ2r('fh5 ',FHMXF5,'fcf ',fcmxf5)
call writ2r('fo5 ',FoMXF5,'fnf ',fnmxf5)
CALL RQ1R('7GAS','HF1',HMXF5)
ENDIF
IF(LBMXF6.NE.0) THEN
CALL RQ1R('7GAS','FH6',FHMXF6)
CALL RQ1R('7GAS','FC6',FCMXF6)
CALL RQ1R('7GAS','FO6',FOMXF6)
CALL RQ1R('7GAS','FN6',FNMXF6)
CALL RQ1R('7GAS','HF1',HMXF6)
ENDIF
IF(LBMXF7.NE.0) THEN
CALL RQ1R('7GAS','FH7',FHMXF7)
CALL RQ1R('7GAS','FC7',FCMXF7)
CALL RQ1R('7GAS','FO7',FOMXF7)
CALL RQ1R('7GAS','FN7',FNMXF7)
CALL RQ1R('7GAS','HF1',HMXF7)
ENDIF
IF(GRNDNO(7,RHO1)) THEN
RHOAG=RHO1A ; RHOBG=RHO1B ; RHOCG=RHO1C
ELSE
RHOAG=RHO2A ; RHOBG=RHO2B ; RHOCG=RHO2C
ENDIF
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRITST
CALL WRIT40('Initial print-out for RHO1=7GASES model:')
CALL WRIT3R('CINCL ',RHOAG,'HINCL ',RHOBG,
1 'NINCL ',RHOCG)
ENDIF
C.... The heat of combustion per unit mass of CO:
HCCO2=32.792E6 ; HCCO=9.208E6 ; HHH2O=120.9E6
HCOCO2=(12.0/28.0) * (HCCO2 - HCCO)
HCHX=RHOAG*HCCO2 + RHOBG*HHH2O
HVOL=HCHX
C.... Assume that the volatile gas is propane:
AMWTVL=44.
IF(.NOT.NULLPR)
1 CALL WRIT3R('HCOCO2 ',HCOCO2,'HHH2O ',HHH2O,'NCHX ',HCHX)
IF(.NOT.ONEPHS) THEN
CALL GXMAKE(L0H1IF,NXNYST,'H1IF')
CALL GXMAKE(L0H2IF,NXNYST,'H2IF')
ENDIF
C.... Memory allocation for the debug print-out:
if(dbrho) then
call gxmake(l0hsub,nxnyst,'hsub')
call gxmake(l0vola,nxnyst,'vola')
endif
C.... Model for the 'wood' burning, i.e. any solid material which shares
c the gas-phase velocities:
cccc LBFWD=LBNAME('FWD ') ! storage of FWD is the signal
WOOD=LBFWD.NE.0
IF(.NOT.NULLPR) CALL WRIT1L('WOOD ',WOOD)
IF(WOOD) THEN
cccc LBWOOD=LBNAME('WOOD') ! but there may be 2 derivatives:
cccc LBCHAR=LBNAME('CHAR') ! 'wood' and 'char'
C.... Get data for the wood from SPEDAT:
C.... MOINWD, the moisture fraction of the wood, is defined as the
C amount of water in unit mass of MOIST wood:
CALL GETSDR('WOODBURN','MOINWD ',MOINWD) ! read in moisture in wood
C.... CINWD, HINWD and VLINWD denote the mass fractions of
c carbon, hydrogen and volatiles in DRY wood:
CALL GETSDR('WOODBURN','VLINWD',VLINWD) ! read in volatile
CALL GETSDR('WOODBURN','HINVL ',HINVL) ! and its hydrogen
CINVL=1.- HINVL ! but its carbon
HINWD=HINVL*VLINWD*(1.- MOINWD) +
1 2./18.*MOINWD ! deduce H in vol and h2o
CINWD=(1.- HINWD)*(1.- MOINWD) ! C in wood
CHINWD=(1.- VLINWD)*(1.- MOINWD) ! H in wood
OINWD=16./18.*MOINWD ! assume all O in H2O
C.... Heat of combustion:
HVOL=HHH2O*HINVL + HCCO2*RHOAG
HCHAR=HCCO2
HWOOD=HVOL*VLINWD + HCHAR*(1.- VLINWD)
C.... Molecular weight of volatiles:
CALL GETSDR('WOODBURN','MLWTVL',AMWTVL)
IF(.NOT.NULLPR) THEN
CALL WRIT3R('VLINWD ',VLINWD,'HINVL ',HINVL,
1 'MLWTVL ',AMWTVL)
CALL WRIT3R('CINVL ',RHOAG,'HINWD ',HINWD,
1 'CINWD ',CINWD )
CALL WRIT3R('HWOOD ',HWOOD, 'HVOL ',HVOL,
1 'HCHAR ',HCHAR)
CALL WRIT2R('CHENWD ',CHINWD,'MOINWD ',MOINWD)
ENDIF
C.... Store reciprocal of the volatile gas molecular weight:
RMWTVL=1./(AMWTVL+1.E-10)
ENDIF ! ------------------------------end of data input for 'wood'
IF(.NOT.NULLPR) THEN
CALL WRITST
CALL WRITBL
ENDIF
C.... Note that there is no need to employ a constant specific heat; for
C coding which allows it to depend upon temperature and composition
C could also be introduced into this subroutine; but then the heats
C of combustion would also need to be modified accordingly.
c
C... get index for H0
L0H012=L0F(ITWO(IH01,IH02,IPHASE.EQ.1))
ENDIF
C.... Consistency check for for compressibility options:
IF(GRNDNO(5,DRH1DP)) THEN
IF(EQZ(RHO1C)) RHO1C=1./1.4
ENDIF
IF(GRNDNO(5,DRH2DP)) THEN
IF(EQZ(RHO2C)) RHO2C=1./1.4
ENDIF
if(dbrho) call banner(2,'INIDEN',220306)
END
c-----------------------------------------------------------------------
c