cgxifric.htm c c
C... File name ... GXIFRIC.HTM .... 261121
      FUNCTION FRICCO(I)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE '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==-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==1) THEN  ! Fip= Const * Rho1 * R1 * R2! * Vol:
        FRICCO= CFIPC*CELMAS(i,1)*AMAX1(RLOLIM,F(L0R+I))
      ELSEIF(IGRND==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/=0) THEN  ! divide by length scale if solved for
          FRICCO=FRICCO / F(L0LEN+I)
        ENDIF
      ELSEIF(IGRND==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==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==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==6) THEN
C.... Fip= Const*Volume:
        FRICCO= CFIPC*CELVOL(i)
      ELSEIF(IGRND==7 .OR. IGRND==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==7) THEN ! Particle-Fluidization Drag Model
          RCON= F(L0R +I)
          RDIS= F(L0RD+I)
          REYD= RCON*REYD
          IF(RCON<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
C
        ELSEIF(IDRAGCO==8) THEN ! Particle-Cluster 4-Zone Fluidization Drag Model
C Li et al, "Drag models for simulating gas-solid flow in turbulent
C             fluidisation of FCC particles", Particuology,7, 269-277, (2009)
          RCON= F(L0R +I) ; RDIS= F(L0RD+I)
          RDISL= AMAX1(RLOLIM,RDIS)
          REYD= RCON*REYD
          DIAMC= ABS(CFIPC) ! Cluster diameter
          REYCL= RCON*DIAMC*VSLP/(F(L0VISL+I)+TINY) + TINY
          R1LIM= 0.8 ; R2LIM= 0.933 ; R3LIM= 0.99
C Dense regime - Ergun correlation - RCON <= 0.8
          CDP= 0.0
          BETA1= RDIS*(150.*F(L0VISL+I)*RDIS/(RCON*DIAMC)+1.75*VSLP)*
     1           F(L0DEN+I)/DIAMC
C Sub-dense regime RCON > 0.8 & RCON <= 0.933)
          CDP=   AMAX1(0.44, 24.*(1.+0.15*REYCL**0.687)/REYCL)
          BETA2= (5./72.)*CDP*RDISL*RCON*VSLP/(RDISL**0.293)*
     1            F(L0DEN+I)/DIAMC
C Sub-dilute regime Wen & Yu -  RCON => 0.933 & RCON <= 0.99)
          CDP= AMAX1(0.44, 24.*(1.+0.15*REYD**0.687)/REYD)
          BETA3= 0.75*CDP*RDISL*VSLP/(RCON**1.65)*
     1           F(L0DEN+I)/DIAM
C Dilute regime Schiller & Neumann   RCON > 0.99
          CDP= AMAX1(0.44, 24.*(1.+0.15*REYD**0.687)/REYD)
          BETA4= 0.75*CDP*RDISL*RCON*VSLP*F(L0DEN+I)/DIAM
C Stitching function for smooth transition between regimes
          PI  = 3.141592654
          PHI1=ATAN(150.*1.75*(RCON-R1LIM))/PI+0.5
          PHI2=ATAN(150.*1.75*(RCON-R2LIM))/PI+0.5
          PHI3=ATAN(150.*1.75*(RCON-R3LIM))/PI+0.5
          BETA=(1-PHI1)*BETA1 + PHI1*((1-PHI2)*BETA2
     1        + PHI2*((1-PHI3)*BETA3 + PHI3*BETA4))
          FRICCO= BETA*CELVOL(I)
          CALL STORVAL(LBREYC,L0REYC,REYCL,I)
        ELSE
C
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==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<=3.38E5) THEN
              CDP= 24.*(1.+ 0.15*REYD**0.687)/REYD +
     1                   0.42/(1.+ 4.25E4/(REYD**1.16))
            ELSEIF(REYD<=4.E5) THEN
              CDP= 29.78 - 5.3*ALOG10(REYD)
            ELSEIF(REYD<=1.E6) THEN
              CDP= 0.1*ALOG10(REYD) - 0.49
            ELSE
              CDP= 0.19-8.0E4/REYD
            ENDIF
          ELSEIF(IDRAGCO==1) THEN ! Dispersed-Solid Drag Model/Stokes
                                     ! Regime: [Re<1]
            CDP= 24./(REYD+TINY)
          ELSEIF(IDRAGCO==2) THEN ! Dispersed-Solid Drag Model/Turbulent
                                     ! Regime: [1.E38.) THEN
                CDP= 2.666
              ELSE
                WEBX= 2065.1/WEBER**2.6
                IF(REYD>WEBX) CDP= 0.3333*WEBER
              ENDIF
              CALL STORVAL(LBWEB,L0WEB,WEBER,I)
            ENDIF
          ELSEIF(IDRAGCO==5) THEN
C.... Spherical-Bubble "Dirty-water" Drag correlation
C                         [ Re> 100]  (see Kuo & Wallis [1988])
            CDP= 6.3/REYD**0.385
          ELSEIF(IDRAGCO==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/=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 'farray'
      IF(LBVALUE/=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 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE 'prpcmn'
      COMMON/CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      COMMON/GENI/IGF1(2),NXNYST,IGF2(52),IPRPS,IGF3(4)
      COMMON /F0/KFIL1(109),KZXCY,KFIL2(194)
      COMMON/VELCMN/L0UVW(6),L0UVW2(6)
      COMMON/NAMFN/NAMFUN,NAMSUB
      LOGICAL dbgloc,SLD,NEZ
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB= 'slbfrc'
      if(flag.or.dbgloc) then
        call banner(1,namsub,050918)
        call writ1i('igrnd   ',igrnd)
      endif
      IGR= 10; ISC= 1
C.... Call GROUND for the user-set property:
      IF(IPILOPT==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
      IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
      IF(SOLVE(3))  L0UVW(3)=L0F(3); IF(SOLVE(4))  L0UVW2(3)=L0F(4)
      IF(SOLVE(5))  L0UVW(1)=L0F(5); IF(SOLVE(6))  L0UVW2(1)=L0F(6)
      IF(SOLVE(7))  L0UVW(5)=L0F(7); IF(SOLVE(8))  L0UVW2(5)=L0F(8)
      IF(STORE(9))  L0R1=L0F(9);     IF(SOLVE(10)) L0R2=L0F(10)
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
      L0VOL= L0F(LVOL)
      IF(IPILOPT==-1) THEN
        IF(CFIPS>=0.0) THEN
          IPHASE=1; L0MAS= L0F(LM1); L0R  = L0R2
        ELSE
          IPHASE=2; L0MAS= L0F(LM2); L0R  = L0R1
        ENDIF
c!!!! note that the go to 700, below, might as well be placed here
      ELSEIF(IPILOPT<5) THEN
        IPHASE=1; L0MAS= L0F(LM1); L0R  = L0F(R2)
      ELSEIF(IPILOPT==5) THEN
        IPHASE=2; L0MAS= L0F(LM2); L0R  = L0F(R1)
      ENDIF
      IF(IPILOPT==-1) GO TO 700
      IF(IPILOPT/=1 .AND. IPILOPT/=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>=3 .AND. IPILOPT<=5) THEN
          L0LEN= L0F(LEN1)
        ELSEIF(IPILOPT==7 .OR. IPILOPT==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==7) THEN
            L0DEN = L0F(DEN1); L0DEND= L0F(DEN2)
            L0R   = L0R1; L0RD  = L0R2
          ELSEIF(IPILOPT==8) THEN
            L0DEN = L0F(DEN2); L0DEND= L0F(DEN1)
            L0R   = L0R2; L0RD  = L0R1
          ENDIF
          IF(LBSIZE/=0) L0DIAM= L0F(LBSIZE)
          IF(LBREYD/=0) L0REYD= L0F(LBREYD)
          IF(LBCD/=0)   L0CD  = L0F(LBCD)
          IF(LBAPRJ/=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)>0.0
          IF(SOLVRS) L0RS= L0F(RS)
C.... Select drag coeff. model:
          IDRAGCO= NINT(CFIPD)
          IF(IDRAGCO==4) THEN
            IF(LBWEB/=0) L0WEB= L0F(LBWEB)
          ELSEIF(IDRAGCO==6) THEN
            IF(LBEOT/=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==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==7 .OR. IPILOPT==8) THEN
        IF(IDRAGCO/=7) THEN
          IF(LBAPRJ/=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>0) THEN
          CALL GROUND
        ENDIF
      ENDIF
      NAMSUB= 'slbfrc'
      if(flag.or.dbgloc) call banner(2,namsub,0)
      END
c----------------------------------------------------------------
      SUBROUTINE INIFRC
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'satgrd'
      INCLUDE '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/=0) THEN
        IF(IBUOYA/=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==7 .OR. IGRCINT==8) GO TO 31
        ENDIF
  30  CONTINUE
C.... Initialization of the slip velocity calculations:
  31  CONTINUE
C.... Initialization of CFIPS calculations:
      IF(IGRCFIP==2.OR.IGRCFIP==4) THEN
        IF(LBVREL==0) CALL GXMAKE(L0VSLP,NX*NY,'VSLP')
      ELSEIF(IGRCFIP>=3 .AND. IGRCFIP<=5) THEN
        IF(LEN1<=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==7 .OR. IGRCFIP==8) THEN
        IF(SOLVE(11) .AND. IGRCFIP==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==0) LBREYD= LBNAME('REYN')
        IF(LBREYC==0) LBREYC= LBNAME('REYC')
        IF(LBVREL==0) CALL GXMAKE(L0VSLP,NX*NY,'VSLP')  ! create slabwise
        IF(LBSIZE==0) CALL GXMAKE(L0DIAM,NX*NY,'DIAM')  ! storage if 3D not
        IF(LBREYD==0) CALL GXMAKE(L0REYD,NX*NY,'REYD')  ! provided
        IF(LBREYC==0.AND.CFIPD==8.) CALL GXMAKE(L0REYC,NX*NY,'REYC')
        IF(LBCD==0)   CALL GXMAKE(L0CD,NX*NY,'CD  ')
        IF(LBAPRJ==0) CALL GXMAKE(L0APRJ,NX*NY,'APRJ')
        IF(LBWEB==0) CALL GXMAKE(L0WEB,NX*NY,'WEBN')
        IF(LBEOT==0) CALL GXMAKE(L0EOT,NX*NY,'EOTN')
      ENDIF
C.... Initialization of CIN calculations:
      IF(IGRCINT==7 .OR. IGRCINT==8) THEN
        IF(.NOT.(IGRCFIP==7.OR.IGRCFIP==8)) THEN       ! create slabwise
          IF(LBSIZE==0) CALL GXMAKE(L0DIAM,NX*NY,'DIAM') ! storage if 3D not
          IF(LBREYD==0) CALL GXMAKE(L0REYD,NX*NY,'REYD') ! provided
        ENDIF
      ENDIF
C.... Preliminaries for the bubble-induced ENUT contribution:
      IF(ENUT<0.0 .AND. NINT(ENUTB)==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==3) THEN
        IF(LBMIXF==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