c
c -- file name gxsurprp.htm 110723
C.... SURHOL is called from SLBproperty to set properties when either the
C VOF, scalar equation method, or height of liquid methods are used for
C flows with inter-fluid boundaries.
c
SUBROUTINE SURHOL(KPROP,IPROP,IFILEP,FL1PRP)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'lbnamer'
COMMON /PRPCMR/CONST(6),WV0,WIDTH,GRDAV,ENT0,TABSPC,HLPC,HHPC
1 /CELPAR/IPHASE,IPROPC,IOPT,IFILEC,KPRP0
1 /GENI/NXNY,IFG1(54),IPRPS,IFG2(4)
LOGICAL VAC,POR,FL1PRP, THREEP
C
L0PRPS= L0F(IPRPS); IPROPC=IPROP; IFILEC=IFILEP; KPRP0=KPROP
ISURN = LBSURN
ISRN2 = LBSRN2
THREEP= ISRN2>0
IF(THREEP) THEN
I1M1=L0F(ISURN)
I2M1=L0F(ISRN2)
ENDIF
DO 20 I= 1,NXNY
IF(VAC(I) .OR. POR(I)) THEN
F(KPROP+I)= TINY
ELSE
IMAT= NINT(F(L0PRPS+I))
IF(IMAT==-1) THEN ! i.e. use domain-material properties
F(KPROP+I)=CELPRP(I)
ELSEIF(IMAT<0) THEN ! cell has been marked as having both materials
IF(THREEP)THEN
F(KPROP+I)= PRPM3(I,I1M1,I2M1) ! by setting PRPS negative
ELSE
F(KPROP+I)= PRPM(I,ABS(F(L0PRPS+I))) ! by setting PRPS negative
ENDIF
IF(IFILEP==4) F(KPROP+I)=-F(KPROP+I)
ELSE
INDTB= INDPRTB(IMAT,0)
PRPVAL= F(INDTB+IFILEP)
IF(PRPVAL>-TINY) THEN
F(KPROP+I)= PRPVAL
IF(IFILEP==4) F(KPROP+I)=-F(KPROP+I)
ELSE
IOPT = NINT(ABS(PRPVAL-GRND))/10
INDTB= INDPRTB(IMAT,IFILEP)
DO IC= 1,NINT(F(INDTB+1))
CONST(IC)= F(INDTB+IC+1)
ENDDO
CALL SETOPT
F(KPROP+I)= CELPRP(I)
ENDIF
ENDIF
ENDIF
20 CONTINUE
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SURHOLCP(KPROP,IPROP,IFILEP,FL1PRP)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'lbnamer'
COMMON /PRPCMR/CONST(6),WV0,WIDTH,GRDAV,ENT0,TABSPC,HLPC,HHPC
1 /CELPAR/IPHASE,IPROPC,IOPT,IFILEC,KPRP0
1 /GENI/NXNY,IFG1(54),IPRPS,IFG2(4)
LOGICAL VAC,POR,FL1PRP
C
L0PRPS= L0F(IPRPS); IPROPC=IPROP; IFILEC=IFILEP; KPRP0=KPROP
DO 20 I= 1,NXNY
IF(VAC(I) .OR. POR(I)) THEN
F(KPROP+I)= TINY
ELSE
IMAT= NINT(F(L0PRPS+I))
IF(IMAT==-1) THEN ! i.e. use domain-material properties
F(KPROP+I)=CELPRP(I)
ELSE
IF(IMAT<0.OR.IMAT==IPRPSA.OR.IMAT==IPRPSB.OR.
1 (LBSRN2>0.AND.IMAT==IPRPSC)) THEN
INDTB= INDPRTB(IPRPSB,0)
ELSE
INDTB= INDPRTB(IMAT,0)
ENDIF
PRPVAL= F(INDTB+IFILEP)
IF(PRPVAL>-TINY) THEN
F(KPROP+I)= PRPVAL
ELSE
IOPT = NINT(ABS(PRPVAL-GRND))/10
INDTB= INDPRTB(IMAT,IFILEP)
DO IC= 1,NINT(F(INDTB+1))
CONST(IC)= F(INDTB+IC+1)
ENDDO
CALL SETOPT
F(KPROP+I)= CELPRP(I)
ENDIF
ENDIF
ENDIF
20 CONTINUE
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE FNPRPS(K1,K2,A,B,I1,I2)
INCLUDE 'farray'
COMMON /IGE/IXF,IXL,IYF,IYL,IGFLL(21)
COMMON/LSG/LSGFL1(91),HOL,LSGFL2(4),SURF,LSGFL3(3)
LOGICAL LSGFL1,LSGFL2,LSGFL3,HOL,SURF
C
IF(IXL<0) RETURN
CALL L0F2(K1,K2,I,I2M1,IADD,'FNPRPS')
IF(HOL) THEN
RLOLM=A+1.E-5
RUPLM=B-1.E-5
ELSE
RLOLM=A+1.E-6
RUPLM=B-1.E-6
ENDIF
DO IX=IXF,IXL
I=I+IADD
DO IY=IYF,IYL
I=I+1
IF(F(I)<100) THEN
IF(F(I+I2M1)<=RLOLM) THEN
F(I)=I1
ELSEIF(F(I+I2M1)>=RUPLM) THEN
F(I)=I2
ELSE
F(I)=-I1*100-I2 -F(I+I2M1)
ENDIF
ENDIF
ENDDO
ENDDO
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... Integer 100 below is used to decode the mixture-property indices.
C It corresponds to the integer used in subroutine FNPRPS
c
FUNCTION PRPM(I,PRPFLG)
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /PRPCMR/CONST(6),WV0,WIDTH,GRDAV,ENT0,TABSPC,HLPC,HHPC
1 /CELPAR/IPHASE,IPROP,IOPT,IFILEP,KPROP
LOGICAL FL1PRP
C.... Decode indices and pick up properties
C.... First material:
IMAT= INT(PRPFLG/100); NGO=1; IMAT1=IMAT
10 INDTB= INDPRTB(IMAT,0)
PRPVAL= F(INDTB+IFILEP)
IF(PRPVAL<-TINY) THEN
IOPT = NINT(ABS(PRPVAL-GRND))/10
INDTB= INDPRTB(IMAT,IFILEP)
DO IC= 1,NINT(F(INDTB+1))
CONST(IC)= F(INDTB+IC+1)
ENDDO
CALL SETOPT
PRPVAL= CELPRP(I)
ENDIF
C.... Second material & CVOL
IF(NGO==1) THEN
PRP1= PRPVAL
CVOL= PRPFLG - 100*IMAT
IMAT= INT(CVOL); NGO=2; IMAT2=IMAT
CVOL= CVOL - IMAT
GO TO 10
ENDIF
PRP2= PRPVAL
C.... Mass-average CP for VOF no mass average of CP
IF(IFILEP==3) THEN
C.... the CP of phase (heavy) is equal to CP of phase (light)
C.... should come out with a constant cp of the mixture = cp of light fluid
! PRP1=PRP2
PRP2=PRP1
CVOL=0.0
ELSEIF(IFILEP==2)THEN
C.... Mixture of the dynamic viscosities
RHOL= F(INDPRTB(IPRPSA,0)+1)
RHOG= F(INDPRTB(IPRPSB,0)+1)
CRHO=CVOL*RHOL + (1.0-CVOL)*RHOG
IF(CRHO>0.0)THEN
PRP1=PRP1*RHOG/CRHO
PRP2=PRP2*RHOL/CRHO
ENDIF
ENDIF
C.... Compute mixture properties
PRPM= PRP1*(1.0-CVOL)+PRP2*CVOL
PMIN= MIN(PRP1,PRP2)
PMAX= MAX(PRP1,PRP2)
PRPM= MIN(PRPM,PMAX)
PRPM= MAX(PRPM,PMIN)
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE FNPRPS3(K1,K2,K3,A,B,I1,I2,I3)
C... A C B A C
INCLUDE 'farray'
COMMON /IGE/IXF,IXL,IYF,IYL,IGFLL(21)
COMMON/SLDPRP/SOLPRP,PORPRP,VACPRP,SOLCRT,SOLSAV
C
IF(IXL<0) RETURN
CALL L0F3(K1,K2,K3,I,I2M1,I3M1,IADD,'FNPRPS3')
RLOLM=A+1.E-6
RUPLM=B-1.E-6
DO 1 IX=IXF,IXL
I=I+IADD
DO 1 IY=IYF,IYL
I=I+1
IF(F(I)=RUPLM) THEN ! Full C
F(I)=I3
ELSE ! Mixture B & C
F(I)=-I1*10000-100*I2-I3-F(I+I2M1)-F(I+I3M1)
ENDIF
ELSEIF(F(I+I3M1)<=RLOLM) THEN ! Empty C
IF(F(I+I2M1)<=RLOLM) THEN ! Empty A
F(I)=I1 ! Must be full B
ELSEIF(F(I+I2M1)>=RUPLM) THEN ! Full A
F(I)=I2
ELSE ! Mixture A & B
F(I)=-I1*10000-100*I2-I3-F(I+I2M1)-F(I+I3M1)
ENDIF
ELSEIF(F(I+I2M1)>=RUPLM) THEN ! Full A
F(I)=I2
ELSEIF(F(I+I3M1)>=RUPLM) THEN ! Full C
F(I)=I3
ELSE ! Mixture A, B & C
F(I)=-I1*10000-100*I2-I3-F(I+I2M1)-F(I+I3M1) ! not needed as never decoded!
ENDIF
ENDIF
1 CONTINUE
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... Integer 100 below is used to decode the mixture-property indices.
C It corresponds to the integer used in subroutine FNPRPS
c
FUNCTION PRPM3(I,K1,K2)
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /PRPCMR/CONST(6),WV0,WIDTH,GRDAV,ENT0,TABSPC,HLPC,HHPC
1 /CELPAR/IPHASE,IPROP,IOPT,IFILEP,KPROP
LOGICAL FL1PRP
C.... Decode indices and pick up properties
C.... First material:
IMAT1= IPRPSB !Light Fluid
IMAT2= IPRPSA !Heavy Fluid 1
IMAT3= IPRPSC !Heavy Fluid 2
PRP1 = F(INDPRTB(IMAT1,0)+IFILEP)
PRP2 = F(INDPRTB(IMAT2,0)+IFILEP)
PRP3 = F(INDPRTB(IMAT3,0)+IFILEP)
CVOL2= F(I+K1)
CVOL3= F(I+K2)
CVOL2=AMAX1(0.0,CVOL2)
CVOL2=AMIN1(1.0,CVOL2)
CVOL3=AMAX1(0.0,CVOL3)
CVOL3=AMIN1(1.0,CVOL3)
DO ITER=1,5
CVOL1=AMAX1(0.0,1.0-CVOL2-CVOL3)
CVOL1=AMIN1(1.0,CVOL1)
CVOL2=AMAX1(0.0,1.0-CVOL1-CVOL3)
CVOL2=AMIN1(1.0,CVOL2)
CVOL3=AMAX1(0.0,1.0-CVOL2-CVOL1)
CVOL3=AMIN1(1.0,CVOL3)
ENDDO
IF(IFILEP==3) THEN
C.... the CP of phase (heavy) is equal to CP of phase (light)
C.... should come out with a constant cp of the mixture = cp of light fluid
PRP2=PRP1
PRP3=PRP1
CVOL1=1.0
CVOL2=0.0
CVOL3=0.0
ELSEIF(IFILEP==2)THEN
C.... Mixture of the dynamic viscosities to obtain a mixed Kinematic viscosity
RHOL2= F(INDPRTB(IMAT2,0)+1)
RHOL3= F(INDPRTB(IMAT3,0)+1)
RHOG= F(INDPRTB(IMAT1,0)+1)
CRHO=CVOL1*RHOG + CVOL2*RHOL2 +CVOL3*RHOL3
IF(CRHO>0)THEN
PRP1=PRP1*RHOG/CRHO
PRP2=PRP2*RHOL2/CRHO
PRP3=PRP3*RHOL3/CRHO
ENDIF
ENDIF
C.... Compute mixture properties
PRPM3= PRP1*CVOL1+PRP2*CVOL2+PRP3*CVOL3
PMIN= MIN(PRP1,PRP2,PRP3)
PMAX= MAX(PRP1,PRP2,PRP3)
PRPM3= MIN(PRPM3,PMAX)
PRPM3= MAX(PRPM3,PMIN)
END
c