Appendix G; Listing of GENIUS

    c

    SUBROUTINE GENTRA

    C------------------------------------------------------------------------------C

    C*

    C* This subroutine is part of the GENTRA particle-tracker option

    C* of PHOENICS

    C

    C------------------------------------------------------------------------------C

    INCLUDE '/phoenics/d_includ/satear'

    INCLUDE '/phoenics/d_includ/grdloc'

    INCLUDE '/phoenics/d_includ/satgrd'

    INCLUDE '/phoenics/d_includ/grdbfc'

    INCLUDE '/phoenics/d_includ/grdear'

    INCLUDE '/phoenics/d_includ/bfcear'

    INCLUDE '/phoenics/d_includ/tracmn'

    INCLUDE '/phoenics/d_includ/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

    LASTF=IPRT0

    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

    IF(LG(30).AND.ISTEP.GT.1) CALL FN2(DEN1,C1,RG(1),RHO2-RG(1))

    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)

    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

    IF(LG(30)) THEN

    IF(INDVAR.NE.NPHI) RETURN

    L0FXG=L0F(XG2D)

    L0FYG=L0F(YG2D)

    L0VAL=L0F(VAL)

    DO 1100 IX=1,NX

    DO 1100 IY=1,NY

    I=IY + (IX-1)*NY

    F(L0VAL+I)=RG(1)

    XG=F(L0FXG+I)

    YG=F(L0FYG+I)

    IF((XG + 1.63177*YG).LT.XULAST*RG(6)) F(L0VAL+I)= RHO2

    1100 CONTINUE

    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.STORE(MASS)) 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).AND.STORE(MOMX)) THEN

    CALL FN0(VAL,MOMX)

    ELSEIF((INDVAR.EQ.V1).AND.STORE(MOMY)) THEN

    CALL FN0(VAL,MOMY)

    ELSEIF((INDVAR.EQ.W1).AND.STORE(MOMZ)) THEN

    CALL FN0(VAL,MOMZ)

    ELSEIF((INDVAR.EQ.H1.OR.INDVAR.EQ.JTEM1).AND.

    1 STORE(HEAT)) THEN

    CALL FN0(VAL,HEAT)

    ELSEIF((INDVAR.EQ.VAPO).AND.STORE(MASS)) 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.1) 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

    1921 CONTINUE

    ENDIF

    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 * ------------------- SECTION 8 ---- FINISH OF TIME STEP.

    ELSEIF(ISC.EQ.8) THEN

    IF(GHFILE(1:4).NE.'NONE') THEN

    IF(.NOT.STEADY.AND.(NX.EQ.1.OR.NY.EQ.1.OR.NZ.EQ.1)) THEN

    DO I= 1, NTRACK

    ILOC0= IPFSTR(I)+IPRSTN

    IWRT= ISTEP

    IF(I.EQ.1) IWRT= -ISTEP

    WRITE(LUHIS,'(I4,A,11(1PE13.3))')IWRT,' ',

    + TIM, (F(ILOC0+J),J=1, 10)

    ENDDO

    ENDIF

    ENDIF

    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.AND.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

    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..........8

    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 '/phoenics/d_includ/satear'

    INCLUDE '/phoenics/d_includ/grdloc'

    INCLUDE '/phoenics/d_includ/satgrd'

    INCLUDE '/phoenics/d_includ/grdear'

    INCLUDE '/phoenics/d_includ/grdbfc'

    INCLUDE '/phoenics/d_includ/bfcear'

    INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/satear'

    INCLUDE '/phoenics/d_includ/tracmn'

    INCLUDE '/phoenics/d_includ/grdloc'

    INCLUDE '/phoenics/d_includ/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 =<= 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.

    TEMQ=AMAX1(PARAMT-273.15,100.0)

    GPROPS=1.2745E-07*TEMQ**5-1.2475E-04*TEMQ**4

    + +4.7457E-02*TEMQ**3-8.7110*TEMQ**2+773.98*TEMQ-2.4505E+04

    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.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 GWYOUT(2)

    ENDIF

    if(dbglev.and.dbgrnd) call writ1r('gprops',gprops)

    END

    c


Top Previous Next