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