Appendix G; Listing of GENIUS
C.... File name ................... GENTRA.FOR .................. 120816
SUBROUTINE GENTRA
C------------------------------------------------------------------------------C
C*
C* This subroutine is part of the GENTRA particle-tracker option
C* of PHOENICS
C
C------------------------------------------------------------------------------C
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdbfc'
INCLUDE 'grdear'
INCLUDE 'bfcear'
INCLUDE 'tracmn'
INCLUDE 'moncom'
COMMON/GENI/IDUM1(10),NWHOLE,IDUM2(31),NFTOT,IDUM3(2),LOOPZ,
1 IDUM4(14)
COMMON/GENCHQ/CHQ(8)
LOGICAL GENCAL,PSTLSW,QNE,QEQ
CHARACTER*80 BUFF(3)
INTEGER HEATZ
C.... ISWCNT is a safeguard for the end of time step
SAVE GENCAL, PSTLSW, LGSWEP
C===============================================================================
C
C GROUP 1. Run title
C =========
IF(IGR.EQ.1) THEN
C * ----------------- SECTION 2 ------ GENTRA Preparation
IF(ISC.EQ.1) THEN
CALL GENPRE
GENCAL=.FALSE.
PI=3.1415927
IF(LG(10)) THEN
DO 101 IS=1,8
101 CHQ(IS)=0.0
ENDIF
C * ----------------- SECTION 2 ------
ELSEIF(ISC.EQ.2) THEN
C
C.... NFTOT is the last used storage in F array
IPRT0=MAX0(NFTOT,NBFTOT,NF1TOT)
C..... 0 location for first particle
NOTGXM= .NOT.(TSTSWP.EQ.12345.OR.TSTSWP.EQ.10001.OR.
+ TSTSWP.LT.0)
ENDIF
C
C GROUP 9. Properties of the medium (or media)
C ===================================
ELSEIF(IGR.EQ.9) THEN
C * ------------------- SECTION 1: Density of phase 1 (DEN1)
IF(ISC.EQ.1.AND.QEQ(RHO1,GRND)) THEN
cccc IF(LG(30).AND.ISTEP.GT.1) CALL FN2(DEN1,C1,RG(1),RHO2-RG(1))
! lg(30) relates only to library case g722
! it sets the density = rg(1) i.e. density of air
! + c1*(rho2 - rg(1)) i.e. vol fraction of liquid *
! density difference
! In-Form can now handle this more directly
c
IF(STORE(VAPO).AND.STORE(DEN1).AND.
1 (STORE(JTEM1).OR.STORE(H1))) THEN
C.... Gas constant
FN2A=8314.0/GMWCON
FN2B=8314.0/GMWVAP-FN2A
L0FVAP=L0F(VAPO)
IF(STORE(JTEM1))L0FTEM=L0F(JTEM1)
IF(STORE(H1))L0FH1=L0F(H1)
L0FDEN=L0F(DEN1)
L0FPRE=L0F(P1)
IF(STORE(VPOR))L0FVPO=L0F(VPOR)
IF(STORE(PROPS))L0FPRP=L0F(PROPS)
DO 9010 J=1,NY
DO 9010 I=1,NX
JNYIM1=J+NY*(I-1)
F(L0FDEN+JNYIM1)= 1.0
C.... Is this is a blocked cell?
CALL SUB2R(SOLREG,0.0,CHKBLK,1.1)
IF(STORE(VPOR)) CHKBLK=F(L0FVPO+JNYIM1)
IF(STORE(PROPS))SOLREG=F(L0FPRP+JNYIM1)
IF(SOLREG.LT.REAL(PRINDX).AND.CHKBLK.GT.GPOROS) THEN
C.... the vapour mass fraction,
VAPOUR=F(L0FVAP+JNYIM1)
C.... Find the local temperature,
IF(STORE(JTEM1)) THEN
CTEMPR=F(L0FTEM+JNYIM1)
ELSEIF(STORE(H1)) THEN
C.... Currently this is only for constant Cp. The user can
C change it, so that Cp will be a function of local
C and vapour mass fraction. Iteration may be needed for
C that purpose.
CTEMPR=TMP1A+TMP1B*F(L0FH1+JNYIM1)
ENDIF
C.... the gas constant for the mixture
GASCOS=FN2A+FN2B*VAPOUR
C.... the pressure
PRESSR=F(L0FPRE+JNYIM1)
C.... and finally the density
F(L0FDEN+JNYIM1)=(PRESSR+PRESS0)/(GASCOS*CTEMPR)
ENDIF
9010 CONTINUE
ENDIF
C * ------------------- SECTION 10: Temperature
ELSEIF(ISC.EQ.10.AND.INT(TMP1).EQ.INT(GRND)) THEN
IF(STORE(JTEM1).AND.SOLVE(H1)) THEN
L0FTEM=L0F(JTEM1)
L0FH1=L0F(H1)
IF(STORE(VAPO))L0FVAP=L0F(VAPO)
DO 9110 J=1,NY
DO 9110 I=1,NX
JNYIM1=J+NY*(I-1)
CENTHP=F(L0FH1+JNYIM1)
CTEMPR=F(L0FTEM+JNYIM1)
VAPOUR=0.0
IF(STORE(VAPO)) VAPOUR=F(L0FVAP+JNYIM1)
GCPGAS=(1.0-VAPOUR)*GPROPS(16, CTEMPR, GCPCON)+
+ VAPOUR*GPROPS(15, CTEMPR, GCPVAP)
F(L0FTEM+JNYIM1)=GPROPS(6,CENTHP,GRND1)
9110 CONTINUE
ENDIF
ENDIF
C
ELSEIF(IGR.EQ.11) THEN ! the following appears to be
c IF(LG(30)) THEN ! intended for g722 only.
c IF(INDVAR.NE.NPHI) RETURN ! what is nphi? Probably density;
c L0FXG=L0F(XG2D) ! but only if that is, nas in 227,
c L0FYG=L0F(YG2D) ! the first stored
c L0VAL=L0F(VAL) ! It appears that density equals
c DO 1100 IX=1,NX ! RG(1)
c DO 1100 IY=1,NY ! unless xg+1.63177*yg <
c I=IY + (IX-1)*NY ! xulast rg(6) when it equals RHO2
c F(L0VAL+I)=RG(1) ! In case g722,
c XG=F(L0FXG+I) ! rg(1)=1.189, rg(6=)=0.058
c YG=F(L0FYG+I) ! rho2=998.2
c ! altogether a curious mixture. Presumably the 1.63177 and 0.058 are related
c ! in some way; but what?
c IF((XG + 1.63177*YG).LT.XULAST*RG(6)) F(L0VAL+I)= RHO2
c 1100 CONTINUE
c ENDIF
C GROUP 13. Boundary conditions and special sources
C =======================================
ELSEIF(IGR.EQ.13) THEN
C * ------------------ SECTION 12: GENTRA sources
IF(ISC.EQ.12) THEN
IF(NPATCH(1:6).EQ.'GENMAS'.OR.NPATCH(1:6).EQ.'GENPAT') THEN
IF(ISWEEP.GE.GSWEP1) THEN
IF(NPATCH(1:6).EQ.'GENMAS') THEN
C.... Mass source
IF(INDVAR.EQ.P1.AND.MASS>0) THEN
CALL FN0(VAL,MASS)
CALL GETCOV('GENMAS',1,GCOV,GVAL)
IF(QNE(GCOV,FIXFLU)) THEN
GMULT=1./GCOV
CALL FN25(VAL,GMULT)
ENDIF
ENDIF
ELSEIF(NPATCH(1:6).EQ.'GENPAT') THEN
C.... Other interfacial sources
IF((INDVAR.EQ.U1.OR.INDVAR.EQ.LBNAME('UC1')).AND.
1 MOMX>0) THEN
CALL FN0(VAL,MOMX)
ELSEIF((INDVAR.EQ.V1.OR.INDVAR.EQ.LBNAME('VC1')).AND.
1 MOMY>0) THEN
CALL FN0(VAL,MOMY)
ELSEIF((INDVAR.EQ.W1.OR.INDVAR.EQ.LBNAME('WC1')).AND.
1 MOMZ>0) THEN
CALL FN0(VAL,MOMZ)
ELSEIF((INDVAR.EQ.H1.OR.INDVAR.EQ.JTEM1).AND.
1 HEAT>0) THEN
CALL FN0(VAL,HEAT)
ELSEIF((INDVAR.EQ.VAPO).AND.MASS>0) THEN
CALL FN0(VAL,MASS)
ENDIF
ENDIF
IF(.NOT.STEADY) CALL FN25(VAL,1.0/DT)
ELSE
CALL FN1(VAL,0.0)
ENDIF
ENDIF
ENDIF
C
C
C GROUP 19. Special calls to GROUND from EARTH
C ==================================
C
ELSEIF(IGR.EQ.19) THEN
C * ----------------- SECTION 1 ---- START OF TIME STEP.
IF(ISC.EQ.1) THEN
LGSWEP=0
PSTLSW=.FALSE.
IF(ISTEP.EQ.FSTEP) THEN
C.... Domain size
CALL SUB3R (XLASTM,XULAST,YLASTM,YVLAST,ZLASTM,ZWLAST)
CELMIN=AMIN1(XLASTM/FLOAT(NX),YLASTM/FLOAT(NY),
$ ZLASTM/FLOAT(NZ))
C.... Find the axis
CALL FNDAXI
ENDIF
CALL GENIUS (1,2)
C * ------------------- SECTION 2 ---- START OF SWEEP.
ELSEIF(ISC.EQ.2) THEN
C.. A. Index of first actual sweep...........................
IF(GFASWP.EQ.-9999) THEN
GFASWP=ISWEEP
GSWEP1=MAX0(GFASWP,GSWEP1)
ENDIF
C.. B. Last sweep?...........................................
C At least two sweeps are needed for one time step
LSTSWP=(ISWEEP.GE.LSWEEP).OR.(RESMET.AND.ISWEEP.GT.1)
GENCAL=(ISWEEP.GE.GSWEP1).AND.
+ ((MOD(ISWEEP,GSWEPF).EQ.0).OR.(GSWEP1.EQ.ISWEEP).OR.LSTSWP)
C.. D. Reset sources ........................................
IF(GENCAL) THEN
DO 1921 IZZ=1,NZ
C.... Setting interphase sources to 0 in the first GENTRA sweep
IF(ISWEEP.EQ.GSWEP1) THEN
IF(.NOT.LG(11).AND.NOTGXM.AND.IZZ.EQ.1) THEN
WRITE(BUFF(1),'(A)')
+ 'GENTRA resetting interphase sources to 0'
CALL PRINT_CHECK(BUFF,1,LUPRO)
ENDIF
IF(STORE(MOMX)) THEN
MOMXZ=ANYZ(MOMX,IZZ)
CALL FN1(MOMXZ,0.0)
ENDIF
IF(STORE(MOMY)) THEN
MOMYZ=ANYZ(MOMY,IZZ)
CALL FN1(MOMYZ,0.0)
ENDIF
IF(STORE(MOMZ)) THEN
MOMZZ=ANYZ(MOMZ,IZZ)
CALL FN1(MOMZZ,0.0)
ENDIF
IF(STORE(HEAT)) THEN
HEATZ=ANYZ(HEAT,IZZ)
CALL FN1(HEATZ,0.0)
ENDIF
IF(STORE(MASS)) THEN
MASSZ=ANYZ(MASS,IZZ)
CALL FN1(MASSZ,1.0E-20)
ENDIF
ENDIF
IF(((GRESTI.NE.0).AND.LSTSWP).AND.STORE(REST)) THEN
IF(LOOPZ.EQ.1) THEN
JREST=ANYZ(REST,IZZ)
CALL FN1(JREST,0.0)
ENDIF
ENDIF
IF(STEADY.AND.(LSTSWP.OR.ENUFSW).AND.
1 (STOMCO.OR.STOVFR)) THEN
IF(LOOPZ.EQ.1) THEN
IF(STOMCO) THEN
JPMCO=ANYZ(PMCO,IZZ)
CALL FN1(JPMCO,0.0)
ENDIF
IF(STOVFR) THEN
JPVFR=ANYZ(PVFR,IZZ)
CALL FN1(JPVFR,0.0)
ENDIF
ENDIF
ENDIF
1921 CONTINUE
C
ENDIF
C
IF(.NOT.STEADY.AND.LSTSWP.AND..NOT.PSTLSW) THEN
DO 1922 IZZ=1,NZ
IF(LOOPZ.EQ.1) THEN
IF(STOMCO) THEN
JPMCO=ANYZ(PMCO,IZZ)
CALL FN1(JPMCO,0.0)
ENDIF
IF(STOVFR) THEN
JPVFR=ANYZ(PVFR,IZZ)
CALL FN1(JPVFR,0.0)
ENDIF
ENDIF
1922 CONTINUE
ENDIF
C
C * ------------------- SECTION 3 ---- START OF IZ SLAB.
ELSEIF(ISC.EQ.3) THEN
C ...............................................................
C PSICEL is the main particle-tracking module. It will
C calculate the particle trajectories and the interphase
C sources of momentum, heat and mass at each cell.
C Its CALL is conditioned to GENCAL, set above
C................................................................
C PSICEL protected against double calling in the same sweep
IF(.NOT.PSTLSW.AND.IZ.EQ.1.AND.GENCAL.AND.
+ ISWEEP.NE.LGSWEP) THEN
IF(.NOT.LG(11).AND.NOTGXM) THEN
WRITE(BUFF(1),'(A,I4)')
+ 'GENTRA now tracking particles at sweep',ISWEEP
CALL PRINT_CHECK(BUFF,1,LUPRO)
ENDIF
if(dbggen) then
call writ8('GENTRA..')
call writ1r('TIME= ',TIM)
call writ1i('ISWEEP',ISWEEP)
endif
CALL PSICEL
IF(.NOT.LG(11).AND.NOTGXM) THEN
IF(LUPRO.NE.6) THEN
WRITE(LUPRO,'(A)')'GENTRA returns control to Earth'
ELSE
CALL PUT_LINE('GENTRA returns control to Earth',.true.)
ENDIF
ENDIF
IF(LSTSWP) PSTLSW=.TRUE.
LGSWEP=ISWEEP
ENDIF
C * ------------------- SECTION 6 ---- FINISH OF IZ SLAB.
cc ELSEIF(ISC.EQ.6) THEN
C * ------------------- SECTION 7 ---- FINISH OF SWEEP.
ELSEIF(ISC.EQ.7) THEN
C... Test whether the RESREF criterion was met before LSWEEP
CALL LASTSW
C... Compute mixture density
IF((STOMXD.OR.STOMFR).AND.(LSTSWP.OR.ENUFSW)) THEN
DO 1972 IZZ=1,NZ
IF(LOOPZ.EQ.1) THEN
IF(STOVFR.AND.STOMCO) THEN
JPMCO=ANYZ(PMCO,IZZ)
JPVFR=ANYZ(PVFR,IZZ)
JRHMX=ANYZ(RHMX,IZZ)
C... Particle mixture density - RHMX = (1.-PVF)*RHC + PMCO
IF(DEN1.GT.0) THEN
JRHC =ANYZ(DEN1,IZZ)
CALL FN21(JRHMX,JPVFR,JRHC,0.0,-1.0)
CALL FN34(JRHMX,JRHC,1.0)
CALL FN34(JRHMX,JPMCO,1.0)
ELSE
CALL FN10(JRHMX,JPMCO,JPVFR,RHO1,1.0,-RHO1)
ENDIF
ELSE
call writ40('Mixture density requires STORE(PMCO) ')
call writ40('and STORE(PVFR) ')
call writ40('run terminated')
CALL SET_ERR(509, 'Error. See result file.',1)
ENDIF
C.. Particle mass fraction
IF(STOMCO.AND.STOMXD) THEN
JPMFR=ANYZ(PMFR,IZZ)
CALL FN15(JPMFR,JPMCO,JRHMX,0.0,1.0)
ELSE
call writ40('Particle mass fraction requires ')
call writ40('STORE(RHMX) & STORE(PMCO) ')
call writ40('run terminated')
CALL SET_ERR(510, 'Error. See result file.',1)
ENDIF
ENDIF
1972 CONTINUE
ENDIF
C * ------------------- SECTION 8 ---- FINISH OF TIME STEP.
ELSEIF(ISC.EQ.8) THEN
IF(LG(30)) THEN
LXU2D=L0F(XU2D)
LYV2D=L0F(YV2D)
L0C1=L0F(C1)
IFACE= IG(30)
CALL TRKDEN(F(L0C1+1),LXU2D,LYV2D,IFACE)
ENDIF
IF(GENCAL) THEN
C.... write GENTRA restart file when PHI dumped
IF(.NOT.STEADY.AND.IDISPA.GT.0)
1 CALL WRSTRTS(ISTEP,IDISPA,IDISPB,IDISPC,PNAM(1:1))
IF(ISTEP.EQ.LSTEP) THEN
C.... write GENTRA restart file
IF(.NOT.STEADY.AND.NTRACK.GT.0) CALL WRSTRT
C.... GENTRA check
IF(LG(10)) THEN
RNTRAK=REAL(NTRACK)
RESIDL=0.0
ICK=6
IF(TSOLVE)ICK=7
IF(VAPSOL)ICK=8
DO 115 IS=1,ICK
IRG=50+IS
RESIDL=RESIDL+ABS((CHQ(IS)-RG(IRG))/(RG(IRG)+1.E-7))
115 CONTINUE
IF(RESIDL.GT.0.2) THEN
WRITE(BUFF(1),'(A)')'CHECK !!'
WRITE(BUFF(2),1150)(50+IS,CHQ(IS),IS=1,6)
WRITE(BUFF(3),1160)(50+IS,CHQ(IS),IS=7,8)
CALL PRINT_CHECK(BUFF,3,LUPRO)
ENDIF
ENDIF
IF(LG(11)) THEN
IF(NOTGXM) CALL GRCLZZ
CALL CLOSZZ(44)
IF(LG(12)) CALL CLOSZZ(43)
ENDIF
ENDIF
ENDIF
ENDIF
C===============================================================================
ENDIF
1150 FORMAT(2('RG(',I2,')=',1PE9.2,';'),'RG(',I2,')=',1PE9.2)
1160 FORMAT('RG(',I2,')=',1PE9.2,';RG(',I2,')=',1PE9.2)
END
C
C ........1.........2.........3.........4.........5.........6.........7.........
C
C-------------------------------------------------------------------------------
c
C Subroutine GENIUS:
C GENtra Interface for User Sequences
C === = = =
C
C Examples included:
C GENEX1: Time-step statistics
C GENEX2: Automatic generation of USE file for PHOTON
C-------------------------------------------------------------------------------
SUBROUTINE GENIUS (IGENGR, IGENSC)
C.....Earth and GENTRA data imported via COMMONs in INCLUDE
C files
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
INCLUDE 'bfcear'
INCLUDE 'tracmn'
C===============================================================================
C
C
C GROUP 1: Preliminaries
C =============
IF(IGENGR.EQ.1) THEN
C * -------------- Section 1: Beginning of current PHOENICS session
C The property index used to identify solid region (INTEGER)
C in conjugate-heat transfer cases. Consult CHAM if uncertain.
PRINDX=100
IF(IGENSC.EQ.1) THEN
C * -------------- Section 2: Beginning of current Eulerian time step
ELSEIF(IGENSC.EQ.2) THEN
ENDIF
C
C GROUP 2: Start of new track
C ==================
ELSEIF(IGENGR.EQ.2) THEN
C
C GROUP 3: Start of new Largrangian time-step for the current track
C ========================================================
C
ELSEIF(IGENGR.EQ.3) THEN
IF(LG(30)) THEN
IF(IP.EQ.1) THEN
UCPARN= 0.0
VCPARN= 0.73*UCGASN
ELSEIF(IP.EQ.NTRACK) THEN
UCPARN= UCGASN
VCPARN= 0
ENDIF
ENDIF
C
C GROUP 4: Particle reaches cell boundary
C ==============================
ELSEIF(IGENGR.EQ.4) THEN
C ----------------------------------------------------------------
C IGENSC values as follows:
C 1-Exit 2-Wall 3-Axis or symmetry surface 4-New cell
C ----------------------------------------------------------------
C
C GROUP 5: End of the current Lagrangian time step
C ========================================
ELSEIF(IGENGR.EQ.5) THEN
C
C GROUP 6: End of current track
C ====================
ELSEIF(IGENGR.EQ.6) THEN
C
C GROUP 7: GENTRA returns control to Earth
C ===============================
ELSEIF(IGENGR.EQ.7) THEN
C
C GROUP 8: Special calls
C =============
ELSEIF(IGENGR.EQ.8) THEN
IF(IGENSC.EQ.1) THEN
C * ----------- Section 1: Momentum equation
C Give constants GVCSAA, GVCSCX, GVCSCY and GVCSCZ
C in the energy equation
C
C
C / Up / / Ug / / Up / / GVCSCX /
C d | | | | | | | |
C ----| Vp | = GVCSBB*| Vg | - GVCSAA*| Vp | + | GVCSCY |
C dt | | | | | | | |
C / Wp / / Wg / / Wp / / GVCSCZ /
C
C
C Up
C Vp : The particle velocity
C Wp
C
C Ug
C Vg : The gas velocity
C Wg
C
ELSEIF(IGENSC.EQ.2) THEN
C * ----------- Section 2: Energy equation
C Give constants GHCSAA, GHCSBB and GHCSCC
C in the energy equation
C
C dTp
C ------- = GHCSBB*Tgas - GHCSAA*Tp + GHCSCC
C dt
C
C M: Particle mass
C T: Temperature
C
C The default setting:
C A=B1=interface heat transfer coefficient
C B1=Latent heat du to evaporation
C Disregarding default setting by reseting A, B1 and B2 to zero.
C
ELSEIF(IGENSC.EQ.3) THEN
C * ------------- Section 3: Particle mass equation
C
C dD**2
C ------- = [GMCSCC]-[GMCSAA]*(D**2)
C dt
C
C D: Particle diameter
C
C
ELSEIF(IGENSC.EQ.4) THEN
C * ------------- Section 4: Formulation for solidification
C
C d(Mass frac. of solid)
C ----------------------- = [GHCSAA]
C d(Temper. of part.)
C
C
ELSEIF(IGENSC.EQ.5) THEN
C * ------------- Section 5: User's equations
ENDIF
C
C GROUP 9: Particle inlet condition as function of time (TIM)
C ==================================================
ELSEIF(IGENGR.EQ.9) THEN
C... PRVLIN=?
C
C... * ------------- Section 1: Particle X-coordinate
IF(IGENSC.EQ.1) THEN
PRVLIN=2.*TIM
C
C... * ------------- Section 2: Particle Y-coordinate
ELSEIF(IGENSC.EQ.2) THEN
PRVLIN=1.0
C
C... * ------------- Section 3: Particle Z-coordinate
ELSEIF(IGENSC.EQ.3) THEN
PRVLIN=0.2
C
C... * ------------- Section 4: Particle U-velocity (Cartesian components)
ELSEIF(IGENSC.EQ.4) THEN
PRVLIN=0.1
C
C... * ------------- Section 5: Particle V-velocity (Cartesian component)
ELSEIF(IGENSC.EQ.5) THEN
PRVLIN=0.1
C
C... * ------------- Section 6: Particle W-velocity (Cartesian component)
ELSEIF(IGENSC.EQ.6) THEN
PRVLIN=2.0
C
C... * ------------- Section 7: Particle diameter
ELSEIF(IGENSC.EQ.7) THEN
PRVLIN=0.001
C
C... * ------------- Section 8: Liquid density of particle
ELSEIF(IGENSC.EQ.8) THEN
PRVLIN=1000.0
C
C... * ------------- Section 9: Mass flow pass particle inlet (Nozzle etc.)
ELSEIF(IGENSC.EQ.9) THEN
PRVLIN=0.1
C
C... * ------------- Section 10: Number of identical particle parcels.
ELSEIF(IGENSC.EQ.10) THEN
C
C... * ------------- Section 11: Particle temperature
ELSEIF(IGENSC.EQ.12) THEN
PRVLIN=350.0
C
C... * ------------- Section 12: Density of the solid
ELSEIF(IGENSC.EQ.13) THEN
ENDIF
C
C================================ END ========================================
ENDIF
END
C
c
FUNCTION GPROPS(FUNAME, PARAMT, DEFVAL)
C----------------------------------------------------------------
C Functions to calculate liquid and gas properties for GENTRA
C Only one parameter is permited for each function. others
C are passed through 'TRACMN'
C
C FUNAME: character string name of the function
C PARAMT: real parameter of the function
C DEFVAL: real default value
C----------------------------------------------------------------
C
INCLUDE 'farray'
INCLUDE 'tab_mem'
INCLUDE 'satear'
INCLUDE 'tracmn'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INTEGER FUNAME
LOGICAL FINDGR,GRN
CHARACTER*80 BUFF
if(dbglev.and.dbgrnd) then
call writ1i('findx:',funame)
call writ2r('paramt',paramt,'defval',defval)
endif
C
C---------------------------------------------------------------
C Default value or Q1-set flag
C---------------------------------------------------------------
GPROPS=DEFVAL
C
C================= GROUND 1: WATER AND AIR ================
C
IF(GRNDNO(1,DEFVAL)) THEN
C
IF(FUNAME.EQ.1) THEN
C---------------------------------------------------------------
C 1 Drag coefficient
C---------------------------------------------------------------
C gprops = drag coefficient
C defval = q1-set flag
C paramt = particle Reynolds number
C..... Law for spherical particles
GPROPS=(24.0/PARAMT)*(1.0 + 0.15*PARAMT**0.687)+
+ (0.42/(1. + (4.25E + 4*PARAMT**(-1.16))))
ELSEIF(FUNAME.EQ.2) THEN
C---------------------------------------------------------------
C 2 Nusselt number
C---------------------------------------------------------------
C gprops = Nusselt number
C defval = q1-set flag
C paramt = particle Reynolds number
C
C.... Prandtl number
IF(STORE(H1)) THEN
PRNDT=PRNDTL(H1)
ELSEIF(STORE(JTEM1)) THEN
PRNDT=PRNDTL(JTEM1)
ENDIF
IF(GRN(-ABS(PRNDT)) .OR. PRNDT.LT.0.0) PRNDT= 0.7
GPROPS=2.0+0.6*SQRT(PARAMT)*PRNDT**0.3333
C.... Froessling Number
GPRFFF=1.0
IF(VAPSOL.AND.GPRYVS.LT.1.0) THEN
GBM=(GPRYVS-GCONGS)/(1.0-GPRYVS)
IF(GBM.GT.1.E-05)GPRFFF=ALOG(1.0+GBM)/GBM
ENDIF
GPROPS=GPRFFF*GPROPS
ELSEIF(FUNAME.EQ.3) THEN
C---------------------------------------------------------------
C 3 Thermal conductivity of cont. phase.(without vapour)
C---------------------------------------------------------------
C gprops = thermal conductivity of pure continuous phase
C paramt = temperature of the continuous phase
ELSEIF(FUNAME.EQ.4) THEN
C---------------------------------------------------------------
C 4 Thermal conductivity of vapour
C---------------------------------------------------------------
C gprops = thermal conductivity of vapour
C paramt = particle temperature
C
GPROPS=10.0**(9.1E-04*(PARAMT - 373.0) + 1.39)/1000.0
ELSEIF(FUNAME.EQ.5) THEN
C---------------------------------------------------------------
C 5 Latent heat of solidification
C---------------------------------------------------------------
C gprops = latent heat of solidification
C paramt = particle temperature
C
ELSEIF(FUNAME.EQ.6) THEN
C---------------------------------------------------------------
C 6 Temperature of cont. phase as a function of enthalpy
C---------------------------------------------------------------
C gprops = temperature of the continuious phase
C paramt = enthalpy of the continuous phase
GPROPS=PARAMT/GCPGAS
ELSEIF(FUNAME.EQ.7) THEN
C---------------------------------------------------------------
C 7 Enthalpy of cont. phase as a function of temperature
C---------------------------------------------------------------
C gprops = enthalpy of the continuous phase
C paramt = temperature of the continuious phase
GPROPS=(PARAMT-TMP1A)/TMP1B
ELSEIF(FUNAME.EQ.8) THEN
C---------------------------------------------------------------
C 8 Heat capacity (Cp) of the particle
C (Cp of liquid for melt/solidif. particle)
C---------------------------------------------------------------
C gprops = heat capacity of the particle
C (Cp of liquid for melt/solidif. particle)
C paramt = particle temperature (Kelvin)
XPARAM=PARAMT/273.15
GPROPS=-3892.438*XPARAM**3 + 14895.93*XPARAM**2
+ -18779.6*XPARAM + 11993.33
ELSEIF(FUNAME.EQ.9) THEN
C---------------------------------------------------------------
C 9 Latent heat of evaporation
C---------------------------------------------------------------
C gprops = latent heat of evaporation
C paramt = temperature of the particle
C. Derived from curve fit on steam tables (T in Kelvin)
GPROPS=6.2909E+06 - 2.7457E+4*PARAMT + 65.991*PARAMT**2
+ - 5.7653E - 2*PARAMT**3
ELSEIF(FUNAME.EQ.10) THEN
C---------------------------------------------------------------
C 10 Solid fraction
C---------------------------------------------------------------
C gprops = solid fraction
C paramt = temperature of the particle
GPROPS=0.0
GBASE= (GTLIQD-PARAMT)/(GTLIQD-GTSOLD)
IF(GBASE.GT.0.0)GPROPS=GBASE**GMINDX
IF(PARAMT.LE.GTSOLD) GPROPS=1.0
ELSEIF(FUNAME.EQ.11) THEN
C---------------------------------------------------------------
C 11 Index of solid-fraction formula
C---------------------------------------------------------------
C gprops = index of solidification formula
C paramt = local pressure of the continuous phase
C
ELSEIF(FUNAME.EQ.12) THEN
C---------------------------------------------------------------
C 12 Density of the particle
C---------------------------------------------------------------
C gprops = density of particle (for liquid only if GTYPE=50)
C paramt = temperature of the particle
GPROPS=1000.0
ELSEIF(FUNAME.EQ.13) THEN
C---------------------------------------------------------------
C 13 Density of the solid phase of the particle
C---------------------------------------------------------------
C gprops = density of solid
C paramt = temperature of the solid
C
ELSEIF(FUNAME.EQ.14) THEN
C---------------------------------------------------------------
C 14 Saturation press. of vapour.
C---------------------------------------------------------------
C gprops = saturation pressure
C paramt = temperature of the liquid
C
C (Vapour-pressure correlation of Bain [1964]
C Convert from degree Kelvin to Degrees Celcius)
TCLCIS=AMAX1(PARAMT-273.15,1.0)
C Convert Tsat from degC to degR
TDEGR=1.8*TCLCIS+492.
IF(TDEGR.LE.672.0) THEN
C Correlation for T <= 100 degC
TDEGR1=73.32642 - 8.2*ALOG(TDEGR)
+ + 0.003173*TDEGR - 13023.8/TDEGR
ELSE
C Correlation for 100 degC =< T <= Tcrit
ZZZ=TDEGR*TDEGR-951588.
C.... No super-critical state is included
HHH=AMAX1(1165.09 - TDEGR,0.0)
CONS2=2.624453E-12*ZZZ*ZZZ
IF(CONS2.GT.50.) CONS2=50.
CONS3=-0.00631141*HHH**1.25
IF(CONS3.LT.-50.) THEN
CONS3=-50.
ELSEIF(CONS3.GT.50.) THEN
CONS3=50.
ENDIF
TDEGR2=0.00017741*ZZZ*(EXP(CONS2)-1.)/TDEGR
TDEGR3=-0.01013139*EXP(CONS3)
TDEGR1=15.182911-8310.453/TDEGR+TDEGR2+TDEGR3
ENDIF
C Convert Psat from lbf/in2 to N/m2 by multiplying by 6894.76
IF(TDEGR1.GT.50.) TDEGR1=50.
IF(TDEGR1.LT.-50.) TDEGR1=-50.
GPROPS=EXP(TDEGR1)*6894.76
ELSEIF(FUNAME.EQ.15) THEN
C---------------------------------------------------------------
C 15 Heat capacity (Cp) of vapour
C---------------------------------------------------------------
C gprops = heat capacity of vapour
C paramt = temperature
C Cp is calculated under saturation condition
C The temperature should be higher than 100 degree C.
C TEMQ=AMAX1(PARAMT-273.15,100.0)
C GPROPS=1.2745E-07*TEMQ**5-1.2475E-04*TEMQ**4
C + +4.7457E-02*TEMQ**3-8.7110*TEMQ**2+773.98*TEMQ-2.4505E+04
C Cp,vap formula valid for 0.0 < TdegC < 300 polynomial fit to
C data in "Thermodyanamic & Transport Properties of Fluids",
C G. Rogers & Y.Mayhew, pg 10, Basil Blackwell, 4th Edtn. (1987).
TEMQ=PARAMT-273.15
IF(TEMQ.LE.220.) THEN
AP0 = 1856.0
AP1 = 1.157
AP2 = -2.662E-2
AP3 = 4.996E-4
AP4 = -1.802E-6
AP5 = -1.158E-8
AP6 = 1.289E-10
AP7 = -2.943E-13
GPROPS = AP0 + AP1*TEMQ + AP2*TEMQ**2 + AP3*TEMQ**3 +
+ AP4*TEMQ**4 + AP5*TEMQ**5 +
+ AP6*TEMQ**6 + AP7*TEMQ**7
ELSE
GPROPS=13320.-105.4*TEMQ+0.2708*TEMQ**2
ENDIF
ELSEIF(FUNAME.EQ.16) THEN
C---------------------------------------------------------------
C 16 Cp of the continuous phase (without vapour).
C---------------------------------------------------------------
C gprops = heat capacity of pure continuous phase
C paramt = temperature of the gas
C
ELSEIF(FUNAME.EQ.17) THEN
C---------------------------------------------------------------
C 17 Solidus temperature
C---------------------------------------------------------------
C gprops = saturation temperature of solidification
C paramt = pressure of cont. phase
C
ELSEIF(FUNAME.EQ.18) THEN
C---------------------------------------------------------------
C 18 Liquidus temperature
C---------------------------------------------------------------
C gprops = saturation temperature of cont. phase
C paramt = pressure of cont. phase
C
ELSEIF(FUNAME.EQ.19) THEN
C---------------------------------------------------------------
C 19 Particle enthalpy
C (liquid enthalpy for melt/solidif. particle)
C---------------------------------------------------------------
C gprops = particle enthalpy
C (liquid enthalpy for melt/solid. part.)
C paramt = particle temperature
C
XPARAM=PARAMT/273.15
GPROPS=-1.80172E+06-265805.0*XPARAM**4+1.35627E+06*XPARAM**3
+ -2.56482E+06*XPARAM**2+3.27598E+06*XPARAM
C
ELSEIF(FUNAME.EQ.20) THEN
C---------------------------------------------------------------
C 20 Vapour saturation temperature as function of pressure
C---------------------------------------------------------------
C gprops = saturation temperature
C paramt = pressure of the continuous phase
C
XPARAM=ALOG(PARAMT)-5.0
GPROPS=0.31911*XPARAM**3+3.1032*XPARAM**2+27.287*XPARAM+370.8
C
ELSEIF(FUNAME.EQ.21) THEN
C---------------------------------------------------------------
C 21 Cp of solid for solidifying/melting particle
C---------------------------------------------------------------
C gprops = heat capacity of solid
C paramt = particle temperature
C
C
ENDIF
C======================== END OF GROUND 1 ====================
C
C
C======================== GROUND 2: (Examples) ====================
C
C Example functions used for testing and validation.
C
ELSEIF(GRNDNO(2,DEFVAL)) THEN
IF(FUNAME.EQ.1) THEN
C---------------------------------------------------------------
C 1 Drag coefficient
C---------------------------------------------------------------
C gprops = drag coefficient
C defval = q1-set flag
C paramt = particle Reynolds number
C..... Drag Law supplied as table of Cd vs Re values
GPROPS=INTERPFC(ITABNUM,1,2,PARAMT)
ELSEIF(FUNAME.EQ.5) THEN
C---------------------------------------------------------------
C 5 Latent heat of solidification
C---------------------------------------------------------------
C gprops = latent heat of solidification
C paramt = particle temperature
C
GPROPS=-1.3E+06 + 1000.0*PARAMT
ELSEIF(FUNAME.EQ.8) THEN
C---------------------------------------------------------------
C 8 Heat capacity (Cp) of the particle
C---------------------------------------------------------------
C gprops = heat capacity of the liquid
C paramt = liquid temperature (Kelvin)
GPROPS=1000.0 + 0.1*PARAMT
ENDIF
C
C======================== GROUND X: ( ) ====================
C ELSEIF(GRNDNO(X,DEFVAL)) THEN
C-----------------------------------------------------------------
C Properties for other materials can be added in exactly the
C same way as for water. X is from 1 to 10 and GRNDX should be
C given to the request property as default in Q1.
C-----------------------------------------------------------------
C ENDIF
ENDIF
C=========================== END =============================
C
IF(FINDGR(GPROPS)) THEN
WRITE(BUFF,'(A,I2,A)')'GPROPS ',FUNAME,
+ ' is required but no constant or function has been supplied !'
CALL PRINT_CHECK(BUFF,1,LUPRO)
CALL SET_ERR(293, 'Error. See result file.',1)
C CALL EAROUT(2)
ENDIF
if(dbglev.and.dbgrnd) call writ1r('gprops',gprops)
END
Top Previous Next