cgxidif.htm cc... File name ... GXIDIF.HTM ... 080910 REAL FUNCTION INTDIF(I) INCLUDE '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/grdloc' INCLUDE '/phoenics/d_includ/satgrd' INCLUDE '/phoenics/d_includ/satear' INCLUDE '/phoenics/d_includ/grdear' INCLUDE '/phoenics/d_includ/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 GRN c----------------------------------------------------------------------- INTDIF=-999.0 C.... Bulk-to-interface diffusive-transfer coefficient: IF(IGRND.EQ.-1) THEN INTDIF=PRPRTY ELSEIF(IGRND.EQ.1) THEN INTDIF=0.0 ELSEIF(IGRND.EQ.7 .OR. IGRND.EQ.8) THEN C.... Use stored slip velocity, or calculate it (CINH1C stores minimum C slip velocity): IF(STVSLP) THEN VSLP= F(L0VSLP+I) ELSE VSLP= AMAX1(CINH1C, VSLPCL(I,XCYCZ,SLUNX1)) IF(LBVREL.NE.0) F(L0VSLP+I)= VSLP ENDIF IF(L0DIAM.NE.0.AND.L0REYD.NE.0) THEN DIAM= F(L0DIAM+I) REYD= F(L0REYD+I) ELSE C.... Compute particle size Dp= Dp0*(Rd/RS)**0.3333, where Dp0=CINH2C C is the initial particle diameter for the SHADOW method, and C particle Reynolds number ReD= Dp*RelVel/Enul: DIAM= ABS(CINH2C) IF(SOLVRS) C 1 DIAM= DIAM*AMAX1(0., AMIN1(1.,F(L0RD+I)/F(L0RS+I)))*0.3333 1 DIAM= DIAM*AMAX1(0., AMIN1(1.,F(L0RD+I)/F(L0RS+I)))**0.3333 REYD= DIAM*VSLP/(F(L0VISL+I)+TINY) F(L0DIAM+I)= DIAM F(L0REYD+I)= REYD ENDIF C.... Get Prandtl number and diffusivity: PRLM= PRNDTL(INDVAR) IF(GRN(-ABS(PRLM))) THEN C.... GROUND-set Prandtl number (PRLM=GRND) or diffusivity (PRLM=-GRND): IF(PRLM.LT.0.0) THEN PRLM= F(L0PRL+I) DIFF= F(L0VISL+I)/PRLM ELSE DIFF= F(L0PRL+I) PRLM= F(L0VISL+I)/DIFF ENDIF ELSE C.... Constant Prandtl (PRLM>0.0) or diffusivity (PRLM<0.0): IF(PRLM.GT.0.0) THEN DIFF= F(L0VISL+I)/PRLM ELSE DIFF= - PRLM PRLM= F(L0VISL+I)/DIFF ENDIF ENDIF IF(IGRND.EQ.7) THEN C.... Calculate Nusselt number from Clift-Grace-Weber correlation: FAC= ((1.0 + 1.0/(REYD*PRLM))*PRLM)**0.3333 IF(REYD.LE.400) THEN DNUSS= 1.0 + FAC*REYD**0.41 ELSEIF(REYD.LE.2000.) THEN DNUSS= 1.0 + 0.752*FAC*REYD**0.472 ELSE DNUSS= 1.0 + FAC*(0.44*SQRT(REYD) + 0.034*REYD**0.71) ENDIF ELSE C.... Calculate Nusselt number from Ranz-Marshall correlation: DNUSS= 2.0 + 0.6*SQRT(REYD)*PRLM**0.3333 ENDIF IF(LBNUSS.NE.0) F(L0NUSS+I)= DNUSS C.... Compute the surface area of particles (spheres) in each cell: ASD= 6.*F(L0RD+I)*CELVOL(I)/DIAM C.... Set interface transfer coefficient Ci= RHOCG*(Enul*Nu/Prl)/Dp*As: INTDIF= CINHFG*F(L0DEN+I)*DIFF*DNUSS/DIAM*ASD ELSEIF(IGRND.EQ.9) THEN C.... Constant diffusive-transfer coefficient: INTDIF= CINHAG ENDIF END C----------------------------------------------------------------------- SUBROUTINE SLBIDF(IPILOPT,dbgloc) INCLUDE '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/grdloc' INCLUDE '/phoenics/d_includ/satgrd' INCLUDE '/phoenics/d_includ/satear' INCLUDE '/phoenics/d_includ/grdear' INCLUDE '/phoenics/d_includ/prpcmn' COMMON/VMSCMN/FL1CON 1 /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1, 1 L0UD2,L0UD3,L0UD4 1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE C 1 /TSKEM/ GKTDKP,LBKP,LBKT,LBET,LBVOSQ,LBOMEG COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG 1 /TSKEMR/ GKTDKP 1 /VELCMN/L0UVW(6),L0UVW2(6) 1 /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 REAL INTDIF LOGICAL DBGLOC,SLD,GRN,CINT78,FL1CON CHARACTER*6 NAMFUN,NAMSUB C NAMSUB= 'SLBIDF' if(flag.or.dbgloc) call banner(1,namsub,220306) IGR= 10 ISC= IPROP-20 C.... Call GROUND for the user set property: IF(IGRND.EQ.0) THEN IF(USEGRD) THEN CALL GROUND ENDIF GO TO 800 ENDIF C.... For a constant property go to the slab loop IF(IGRND.EQ.-1) GO TO 700 CALL PRPADR C.... Set constants and other auxiliary variables: C----------------------------------------------------------------------- C.... Bulk-to-interface diffusive-transfer coeff.: IF(IPHASE.EQ.1) THEN CINHAG = CINH1A c CONST(2)= CINH1B c CONST(3)= CINH1C c ELSE CINHAG = CINH2A c CONST(2)= CINH2B c CONST(3)= CINH2C ENDIF CINT78= .FALSE. IF(IPILOPT.EQ.7 .OR. IPILOPT.EQ.8) THEN L0VISL= L0F(VISL) IF(GRN(-ABS(PRNDTL(INDVAR)))) L0PRL= L0F(IPRL) IF(LBSIZE.NE.0) L0DIAM= L0F(LBSIZE) IF(LBREYD.NE.0) L0REYD= L0F(LBREYD) IF(LBNUSS.NE.0) L0NUSS= L0F(LBNUSS) C.... First check the option used to calculate Cfip. If it is GRND7 or C GRND8, use values of slip velocity, particle sizes and Reynolds C number calculated for Cfip. IF(CFIPS.GE.0) THEN IGRNDCF=-1 ELSE IGRNDCF= NINT(ABS(CFIPS-GRND))/10 ENDIF CFIP78= IGRNDCF.EQ.7 .OR. IGRNDCF.EQ.8 STVSLP= CFIP78 .AND. LBVREL.NE.0 IF(STVSLP) L0VSLP= L0F(LBVREL) C.... Determine indices for continuous and dispersed phases: IF(CFIP78) THEN IF(IGRNDCF.EQ.7) THEN L0DEN= L0F(DEN1) L0RD = L0R2 ELSE L0DEN= L0F(DEN2) L0RD = L0R1 ENDIF ELSEIF(IPHASE.EQ.1) THEN L0DEN= L0F(DEN1) L0RD = L0R2 ELSE L0DEN= L0F(DEN2) L0RD = L0R1 ENDIF INEIGH = ITWO(INDVAR+1,INDVAR-1,IPHASE.EQ.1) PHNHFG= 1.0 IF(CINT(INEIGH).GT.1.E10) THEN PHNHFG = 0.5 CINT78 = IPHASE.EQ.1 ENDIF ENDIF C----------------------------------------------------------------------- C.... Loop over slab to get and set cell properties: 700 IGRND=IPILOPT IF(IPRPS.EQ.0) THEN C.... One material only DO 60 I= 1,NXNYST 60 F(KPROP+I) = INTDIF(I) ELSE C.... exclude solids DO 70 I= 1,NXNYST IF(SLD(I)) THEN F(KPROP+I)= TINY ELSE F(KPROP+I)= INTDIF(I) ENDIF 70 CONTINUE ENDIF C---------------------------------------------------------------------- C.... Corrections, debug print-out, and other property adjustments: C.... Call GREX to correct a property set above: 800 IF(USEGRX) CALL GREX3 C.... Call ALTPRP for an alternative property setting IF(USEALT) CALL ALTPRP C.... Call GROUND for the user to correct a property set above: IF(USEGRD) THEN IF(IPILOPT.GT.0) CALL GROUND ENDIF NAMSUB= 'slbidf' if(flag.or.dbgloc) call banner(2,namsub,0) END c