cgxifric.htm cC... File name ... GXIFRIC.HTM .... 091101 FUNCTION FRICCO(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 /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP C C.... Note that: c (1) The dimensions of FRICCO are force in cell exerted by one phase on c the other divided by their differnce of velocity c (2) Interphase friction coefficient (R! means the larger of R, RLOLIM) c IF(IGRND.EQ.-1) THEN ! PRPRTY is the value of CFIPS IF(.NOT.STORE(10)) THEN ! specified in the Q1 file, with dimensions FRICCO=PRPRTY*CELVOL(I)! force per unit volume / velocity difference ELSE ! If R2 is stored (the usual case) fricco = cfips*m1*r1*r2 ! if cfips > 0.0, else s12= -cfips*m2*r1*r2 FRICCO=ABS(PRPRTY)*CELMAS(I,IPHASE)*AMAX1(RLOLIM,F(L0R+I)) ENDIF ELSEIF(IGRND.EQ.1) THEN ! Fip= Const * Rho1 * R1 * R2! * Vol: FRICCO= CFIPC*CELMAS(i,1)*AMAX1(RLOLIM,F(L0R+I)) ELSEIF(IGRND.EQ.2) THEN ! Fip= Const * Rho1 * R1 * R2!* Vol * RelVel: VSLP = AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1)) CALL STORVAL(LBVREL,L0VSLP,VSLP,I) FRICCO= CFIPC*CELMAS(I,1)*AMAX1(RLOLIM,F(L0R+I))*VSLP IF(L0LEN.NE.0) THEN ! divide by length scale if solved for FRICCO=FRICCO / F(L0LEN+I) ENDIF ELSEIF(IGRND.EQ.3) THEN ! Fip= Const * Rho1 * R1 * R2! * Vol * ! RelVel ** Pwr CALL STORVAL(LBVREL,L0VSLP,VSLP,I) FRICCO= CFIPC*CELMAS(i,1)*AMAX1(RLOLIM,F(L0R+I))* 1 VSLP**CFIPB ELSEIF(IGRND.EQ.4) THEN ! Fip= Const * Rho1 * R1 * R2! * Vol * ! RelVel ** Pwr1 * Len ** Pwr2 VSLP = AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1)) CALL STORVAL(LBVREL,L0VSLP,VSLP,I) FRICCO= CFIPC*CELMAS(I,1)*AMAX1(RLOLIM,F(L0R+I))* 1 VSLP**CFIPB*F(L0LEN+I)**CFIPD ELSEIF(IGRND.EQ.5) THEN ! Fip= Const * Rho2 * R2 * R1! * Vol * ! RelVel**Pwr1*Len**Pwr2 CALL STORVAL(LBVREL,L0VSLP,VSLP,I) FRICCO= CFIPC*CELMAS(I,2)*AMAX1(RLOLIM,F(L0R+I))* 1 VSLP**CFIPB*F(L0LEN+I)**CFIPD ELSEIF(IGRND.EQ.6) THEN C.... Fip= Const*Volume: FRICCO= CFIPC*CELVOL(i) ELSEIF(IGRND.EQ.7 .OR. IGRND.EQ.8) THEN C.... Fip= 0.75*Rhc*Rc*Rd*Cd*ABS(RelVel)*Vol/Dp. CFIPS=GRND7 assumes C Phase 1 is the 'carrier', and Phase 2 is dispersed, GRND8 assumes C reverse. (NOTE, 1.5 is taken into account in APROJ): VSLP= AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1)) CALL STORVAL(LBVREL,L0VSLP,VSLP,I) C.... Compute particle size Dp= Dp0*(R2/RS)**0.3333, where Dp0=CFIPB C is the initial particle diameter for the SHADOW method, and C particle Reynolds number ReD= Dp*RelVel/Enul: DIAM= ABS(CFIPB) IF(SOLVRS) 1 DIAM= DIAM*AMAX1(0., AMIN1(1.,F(L0RD+I)/F(L0RS+I)))**0.3333 REYD= DIAM*VSLP/(F(L0VISL+I)+TINY) + TINY IF(IDRAGCO.EQ.7) THEN ! Particle-Fluidization Drag Model RCON= F(L0R +I) RDIS= F(L0RD+I) REYD= RCON*REYD IF(RCON.LT.0.8) THEN C.... Ergun correlation for dense fluidized region: CDP= 0.0 FRICCO= RDIS*(150.*F(L0VISL+I)*RDIS/(RCON*DIAM)+1.75*VSLP) ELSE C.... Subcritical correlation for dilute fluidized region (also suitable C for use in pneumatic-conveying applications): CDP= AMAX1(0.44, 24.*(1.+0.15*REYD**0.687)/REYD) FRICCO= 0.75*CDP*AMAX1(RLOLIM,RDIS)*VSLP/RCON**1.65 ENDIF FRICCO= FRICCO*F(L0DEN+I)*CELVOL(I)/DIAM ELSE C.... Compute projected area Apr= 1.5*Vol*/Dp: APRJ= 1.5*CELVOL(I)*AMAX1(F(L0RD+I),RLOLIM)/DIAM C.... Compute drag coefficient: IF(IDRAGCO.EQ.0) THEN C.... Dispersed-Solid Drag Model/Standard drag curve ( see 'Bubbles, C Drops & Particles', R.Clift, J.R.Grace, M.E.Weber, pg111-112, C Table 5.1 eqn (10) & Table 5.2 eqns (H),(I) & (J), Academic Press, C (1978). IF(REYD.LE.3.38E5) THEN CDP= 24.*(1.+ 0.15*REYD**0.687)/REYD + 1 0.42/(1.+ 4.25E4/(REYD**1.16)) ELSEIF(REYD.LE.4.E5) THEN CDP= 29.78 - 5.3*ALOG10(REYD) ELSEIF(REYD.LE.1.E6) THEN CDP= 0.1*ALOG10(REYD) - 0.49 ELSE CDP= 0.19-8.0E4/REYD ENDIF ELSEIF(IDRAGCO.EQ.1) THEN ! Dispersed-Solid Drag Model/Stokes ! Regime: [Re<1] CDP= 24./(REYD+TINY) ELSEIF(IDRAGCO.EQ.2) THEN ! Dispersed-Solid Drag Model/Turbulent ! Regime: [1.E3 100] (see Kuo & Wallis [1988]) CDP= 6.3/REYD**0.385 ELSEIF(IDRAGCO.EQ.6) THEN ! Ellipsoidal-Bubble "Clean-water" ! Drag correlation ! [ 0.1 < Eo < 40 ] ! (see Clift et al [1988]) DELRHO= ABS(F(L0DEN+I) - F(L0DEND+I)) EOTVOS= 9.81/CFIPC*DELRHO*DIAM**2 ! Eotvos number: CDP = 0.622/(1./EOTVOS + 0.235*F(L0DEN+I)/DELRHO) CALL STORVAL(LBEOT,L0EOT,EOTVOS,I) ENDIF FRICCO= 0.5*CDP*APRJ*F(L0DEN+I)*VSLP C.... Multiply by volume fraction of continuous phase if CFIPB > 0. IF(MULBYR) FRICCO= FRICCO*AMAX1(RLOLIM, F(L0R+I)) ENDIF CALL STORVAL(LBSIZE,L0DIAM,DIAM,I) CALL STORVAL(LBREYD,L0REYD,REYD,I) CALL STORVAL(LBCD ,L0CD, CDP ,I) CALL STORVAL(LBAPRJ,L0APRJ,APRJ,I) ENDIF C.... Account for droplet-size variation. c call writ2l('sizvar ',sizvar,'prtsiz ',prtsiz) c call writ3i('igrnd ',igrnd,'l0r2 ',l0r2,'l0rs ',l0rs) IF(PRTSIZ) THEN IF(L0R2*L0RS.NE.0) FRICCO= FRICCO* 1 AMAX1(TINY, AMIN1(1.,F(L0R2+I)/F(L0RS+I)))**(-0.6667) ENDIF END c----------------------------------------------------------------------- SUBROUTINE STORVAL(LBVALUE,L0VAL,VALUE,I) c.... This subroutine places in 3D or slab-wise stores values of various c physical quantities which are needed by more than one of the c property or inter-phase-transport modules. c INCLUDE '/phoenics/d_includ/farray' IF(LBVALUE.NE.0) THEN ! quantity is stored 3D L0VALUE = L0F(LBVALUE) F(L0VALUE+I) = VALUE ELSE F(L0VAL+I) = VALUE ! quantity is stored slabwise ENDIF END c----------------------------------------------------------------------- SUBROUTINE SLBFRC(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/CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP COMMON/GENI/IGF1(2),NXNYST,IGF2(52),IPRPS,IGF3(4) COMMON/NAMFN/NAMFUN,NAMSUB LOGICAL dbgloc,SLD CHARACTER*6 NAMFUN,NAMSUB C NAMSUB= 'slbfrc' if(flag.or.dbgloc) then call banner(1,namsub,200503) call writ1i('igrnd ',igrnd) endif IGR= 10 ISC= 1 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 ENDIF C.... For a constant property go to the slab loop CALL PRPADR C.... Set constants and other auxiliary variables: C----------------------------------------------------------------------- L0VOL= L0F(LVOL) IF(IPILOPT.EQ.-1) THEN IF(CFIPS.GE.0.0) THEN IPHASE=1 L0MAS= L0F(LM1) L0R = L0F(R2) ELSE IPHASE=2 L0MAS= L0F(LM2) L0R = L0F(R1) ENDIF c!!!! note that the go to 700, below, might as well be placed here ELSEIF(IPILOPT.LT.5) THEN IPHASE=1 L0MAS= L0F(LM1) L0R = L0F(R2) ELSEIF(IPILOPT.EQ.5) THEN IPHASE=2 L0MAS= L0F(LM2) L0R = L0F(R1) ENDIF IF(IPILOPT.EQ.-1) GO TO 700 IF(IPILOPT.NE.1 .AND. IPILOPT.NE.6) THEN C.... Set addresses to calculate slip velocity (NOTE!, CONST(1)=CFIPA C stores minimum slip velocity; if 'VREL' is stored, it will be C filled with calculated slip velocities): IF(IPILOPT.GE.3 .AND. IPILOPT.LE.5) THEN L0LEN= L0F(LEN1) ELSEIF(IPILOPT.EQ.7 .OR. IPILOPT.EQ.8) THEN C.... Determine indices for continuous & dispersed phases. GRND7 assumes C Phase 1 is the 'carrier', and Phase 2 is dispersed; GRND8 assumes C the reverse. IF(IPILOPT.EQ.7) THEN L0DEN = L0F(DEN1) L0R = L0R1 L0DEND= L0F(DEN2) L0RD = L0R2 ELSEIF(IPILOPT.EQ.8) THEN L0DEN = L0F(DEN2) L0R = L0R2 L0DEND= L0F(DEN1) L0RD = L0R1 ENDIF IF(LBSIZE.NE.0) L0DIAM= L0F(LBSIZE) IF(LBREYD.NE.0) L0REYD= L0F(LBREYD) IF(LBCD.NE.0) L0CD = L0F(LBCD) IF(LBAPRJ.NE.0) L0APRJ= L0F(LBAPRJ) L0VISL= L0F(VISL) L0VOL = L0F(LVOL) c.... Const(2) may be set in a variety of ways, as indicated by the equivalence c in prpcmn. Which setting is relevant here is not clear cccc equivalence (const(2),rhobg,enulbg,prlhbg, cccc 1 phnhbg,tmpbg,elbg,cinhbg,cvmbg,dvobg,drhbg,cpbg) c.... PROPS(6) would be name than CONST(6). c.... Note that CONST(..) can be filled with values from the PROPS file c in gxprutil c MULBYR= CONST(2).GT.0.0 IF(SOLVRS) L0RS= L0F(RS) C.... Select drag coeff. model: IDRAGCO= NINT(CFIPD) IF(IDRAGCO.EQ.4) THEN IF(LBWEB.NE.0) L0WEB= L0F(LBWEB) ELSEIF(IDRAGCO.EQ.6) THEN IF(LBEOT.NE.0) L0EOT= L0F(LBEOT) ENDIF ENDIF ENDIF C!!!! The implementation of the correction due to droplet- C!!!! size variation needs clarification. It always assumes that R2 C!!!! is the dispersed phase IF(SIZVAR) L0RS= L0F(RS) 700 IGRND=IPILOPT C.... Loop over slab to set cell properties: IF(IPRPS.EQ.0) THEN ! One material only DO 60 I= 1,NXNYST 60 F(KPROP+I)= FRICCO(I) ELSE C.... exclude solids DO 70 I= 1,NXNYST IF(SLD(I)) THEN F(KPROP+I)= TINY ELSE F(KPROP+I)= FRICCO(I) ENDIF 70 CONTINUE ENDIF C---------------------------------------------------------------------- C.... Corrections, debug print-out, and other property adjustments: IF(IPILOPT.EQ.7 .OR. IPILOPT.EQ.8) THEN IF(IDRAGCO.NE.7) THEN IF(LBAPRJ.NE.0) THEN C.... Store projected area per unit volume for print-out etc DO 130 I= 1,NXNYST 130 IF(.NOT.SLD(I)) F(L0APRJ+I)= F(L0APRJ+I)/CELVOL(I) ENDIF ENDIF ENDIF C.... Call GREX to correct a property set above 800 IF(USEGRX) THEN CALL GREX3 ENDIF C.... Call ALTPRP for an alternative property setting IF(USEALT) THEN CALL ALTPRP ENDIF C.... Call GROUND for the user to correct a property set above IF(USEGRD) THEN IF(IPILOPT.GT.0) THEN CALL GROUND ENDIF ENDIF NAMSUB= 'slbfrc' if(flag.or.dbgloc) call banner(2,namsub,0) END c---------------------------------------------------------------- SUBROUTINE INIFRC INCLUDE '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/satear' INCLUDE '/phoenics/d_includ/satgrd' INCLUDE '/phoenics/d_includ/prpcmn' LOGICAL GRNDNO,SOLVE,FLUID1,QEQ,EQZ,GRN c DIMENSION NAMPRP(30) c EQUIVALENCE (NAMPRP(1),DENST1) if(flag) call banner(1,'INIFRC',220104) C C.... Initialize FNVSLP: CALL FNVSLP(0,0,0.0) C.... Density is a linear function of two or three scalars IF((GRNDNO(1,RHO2).OR.GRNDNO(2,RHO2)).AND.IBUOYB.NE.0) THEN IF(IBUOYA.NE.0) THEN CALL WRIT40('3-scalar density formula is available not ') CALL WRIT40('when ONEPHS = F ! ') IBUOYA= 0 ENDIF ENDIF SOLVRS= SOLVE(11) SIZVAR= PRTSIZ .AND. SOLVRS .AND. SOLVE(10) IGRCFIP= 0 IF(GRN(CFIPS)) IGRCFIP= NINT(ABS(CFIPS-GRND))/10 IGRCINT=0 DO 30 MPH= 1,NPHI IF(SOLVE(MPH)) THEN IF(GRN(CINT(MPH))) IGRCINT= NINT(ABS(CINT(MPH)-GRND))/10 IF(IGRCINT.EQ.7 .OR. IGRCINT.EQ.8) GO TO 31 ENDIF 30 CONTINUE C.... Initialization of the slip velocity calculations: 31 CONTINUE C.... Initialization of CFIPS calculations: IF(IGRCFIP.EQ.2.OR.IGRCFIP.EQ.4) THEN IF(LBVREL.EQ.0) CALL GXMAKE(L0VSLP,NX*NY,'VSLP') ELSEIF(IGRCFIP.GE.3 .AND. IGRCFIP.LE.5) THEN IF(LEN1.LE.0 .AND..NOT.GRN(EL1)) THEN CALL WRIT40('CFIPS=GRND3-GRND5 needs storage of LEN1!') CALL SET_ERR(303, * 'Error. CFIPS=GRND3-GRND5 needs storage of LEN1.',1) C CALL EAROUT(1) ENDIF ELSEIF(IGRCFIP.EQ.7 .OR. IGRCFIP.EQ.8) THEN IF(SOLVE(11) .AND. IGRCFIP.EQ.8 .AND..NOT.FLUID1(11)) THEN CALL WRIT40('CFIPS=GRND8 requires TERMS(RS,P,P,P,P,Y,P)') CALL SET_ERR(304, * 'Error. CFIPS=GRND8 requires TERMS(RS,P,P,P,P,Y,P).',1) C CALL EAROUT(1) ENDIF IF(QEQ(CFIPB,0.0)) THEN CALL WRIT40('CFIPS=GRND7/GRND8 require CFIPB = diameter') CALL WRIT40('of particle set in Q1 file! ') CALL SET_ERR(305, 'Error. See result file.',1) C CALL EAROUT(1) ENDIF IF(LBREYD.EQ.0) LBREYD= LBNAME('REYN') IF(LBVREL.EQ.0) CALL GXMAKE(L0VSLP,NX*NY,'VSLP') ! create slabwise IF(LBSIZE.EQ.0) CALL GXMAKE(L0DIAM,NX*NY,'DIAM') ! storage if 3D not IF(LBREYD.EQ.0) CALL GXMAKE(L0REYD,NX*NY,'REYD') ! provided IF(LBCD.EQ.0) CALL GXMAKE(L0CD,NX*NY,'CD ') IF(LBAPRJ.EQ.0) CALL GXMAKE(L0APRJ,NX*NY,'APRJ') IF(LBWEB.EQ.0) CALL GXMAKE(L0WEB,NX*NY,'WEBN') IF(LBEOT.EQ.0) CALL GXMAKE(L0EOT,NX*NY,'EOTN') ENDIF C.... Initialization of CIN calculations: IF(IGRCINT.EQ.7 .OR. IGRCINT.EQ.8) THEN IF(.NOT.(IGRCFIP.EQ.7.OR.IGRCFIP.EQ.8)) THEN ! create slabwise IF(LBSIZE.EQ.0) CALL GXMAKE(L0DIAM,NX*NY,'DIAM') ! storage if 3D not IF(LBREYD.EQ.0) CALL GXMAKE(L0REYD,NX*NY,'REYD') ! provided ENDIF ENDIF C.... Preliminaries for the bubble-induced ENUT contribution: IF(ENUT.LT.0.0 .AND. NINT(ENUTB).EQ.1) THEN IF(EQZ(ENUTC)) ENUTC= 0.6 IF(GRNDNO(3,ENUT)) KELIN= 3 ENDIF C.... Preliminaries for mass-transfer rate: IF(NINT(ABS(CMDOT-GRND))/10.EQ.3) THEN IF(LBMIXF.EQ.0) THEN CALL WRIT40('CMDOT=GRND3 requires storage of MIXF in Q1') CALL SET_ERR(306, * 'Error. CMDOT=GRND3 requires storage of MIXF in Q1.',1) C CALL EAROUT(1) ENDIF ENDIF if(flag) call banner(2,'INIFRC',0) END c