SUBROUTINE RELP(SL,SG,REPL,REPG,REPO,NMAT) C C-----THIS ROUTINE COMPUTES RELATIVE PERMEABILITIES FOR LIQUID C AND GASEOUS PHASES. C C-----ADDED FUNCTION FOR 3PHASE RELATIVE PERMEABILITIES, IRP=9 C------USE PARKER ET AL. (1987) METHOD FALTA 5/89 CCCCC COMMON/P3/DELX(1) COMMON/RPCAP/IRP(27),RP(7,27),ICP(27),CP(7,27),IRPD,RPD(7), 1 ICPD,CPD(7) C SAVE ICALL DATA ICALL/0/ ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'RELP 1.1 ! 12 NOVEMBER 1993',6X, X'RELATIVE PERMEABILITIES AS FUNCTIONS', X' OF SATURATION (INCL 3-PHASE)'/ X47X,'IRP = 6 OPTION ("STONE I") MODIFIED TO FORCE KRN --> 0', X' AS SN --> SNR') C SO=1.0E0-SL-SG REPO=0.0E0 GO TO (10,11,12,12,13,14,15,16,17,14,18),IRP(NMAT) C C RELATIVE OIL AND WATER PERMEABILITIES FOR BUCKLEY-LEVERETT PROBLEM C USED BY FAUST, 1985 C 18 CONTINUE REPG=0.0E0 REPL=((SL-0.16E0)**2)/0.64E0 REPO=((0.8E0-SL)**2)/0.64E0 IF(REPO.LT.0.0E0) REPO=0.0E0 IF(REPL.LT.0.0E0) REPL=0.0E0 RETURN C 10 CONTINUE C-----LINEAR FUNCTIONS. C SLO=SL IF(RP(5,NMAT).GT.0.) SLO=SO REPL=(SLO-RP(1,NMAT))/(RP(3,NMAT)-RP(1,NMAT)) IF (SLO.GE.RP(3,NMAT)) REPL=1.0E0 IF (SLO.LE.RP(1,NMAT)) REPL=0.0E0 REPG=(SG-RP(2,NMAT))/(RP(4,NMAT)-RP(2,NMAT)) IF (SG.GE.RP(4,NMAT)) REPG=1.0E0 IF (SG.LE.RP(2,NMAT)) REPG=0.0E0 IF(RP(5,NMAT).GT.0.) THEN REPO=REPL REPL=0.0E0 END IF RETURN C C-------3PHASE RELATIVE PERMEABILITY OF PARKER ET AL., 1987 C 17 CONTINUE SM=RP(1,NMAT) AN=RP(2,NMAT) AM=1.0E0-1.0E0/AN AMINV=1.0E0/AM C SGBAR=(SG-0.0001E0)/(1.0E0-SM) IF(SGBAR.LT.0.0E0) SGBAR=0.0E0 IF(SGBAR.GT.1.0E0) SGBAR=1.0E0 SWBAR=(SL+0.0001E0-SM)/(1.0E0-SM) IF(SWBAR.LT.0.0E0) SWBAR=0.0E0 IF(SWBAR.GT.1.0E0) SWBAR=1.0E0 STBAR=(SL+0.0001E0+SO-SM)/(1.0E0-SM) IF(STBAR.LT.0.0E0) STBAR=0.0E0 IF(STBAR.GT.1.0E0) STBAR=1.0E0 C RTC=(1.0E0-STBAR**AMINV)**AM RLC=(1.0E0-SWBAR**AMINV)**AM CC PRINT 100, STBAR-SWBAR,RTC,RLC,SL,SG CC100 FORMAT(1X,5(E14.8)) REPG=(SGBAR**0.5E0)*RTC*RTC REPL=(SWBAR**0.5E0)*(1.0E0-RLC)*(1.-RLC) IF(SO.LE.0.) THEN REPO=0. ELSE REPO=((STBAR-SWBAR)**0.5E0)*(RLC-RTC)*(RLC-RTC) END IF C IF(REPG.LT.0.0E0) REPG=0.0E0 IF(REPL.LT.0.0E0) REPL=0.0E0 IF(REPO.LT.0.0E0) REPO=0.0E0 IF(REPG.GT.1.0E0) REPG=1.0E0 IF(REPL.GT.1.0E0) REPL=1.0E0 IF(REPO.GT.1.0E0) REPO=1.0E0 RETURN C C 11 CONTINUE C-----RELATIVE PERMEABILITY OF PICKENS ET AL. C REPG=1.0E0 REPL=(1.0E0-SG)**RP(1,NMAT) C RETURN 12 CONTINUE C-----COREY'S OR GRANT'S CURVES. C SSTAR=(SL-RP(1,NMAT))/(1.-RP(1,NMAT)-RP(2,NMAT)) REPL=SSTAR**4 REPG=(1.0E0-SSTAR**2)*(1.0E0-SSTAR)**2 IF (SG.GE.RP(2,NMAT)) GO TO 50 REPG=0.0E0 REPL=1.0E0 GO TO 102 50 IF (SG.LT.(1.0E0-RP(1,NMAT))) GO TO 102 REPL=0.0E0 REPG=1.0E0 102 CONTINUE IF (IRP(NMAT).EQ.4) REPG=1.0E0-REPL RETURN C 13 CONTINUE C-----ALL PHASES ARE PERFECTLY MOBILE. C REPL=1.0E0 REPG=1.0E0 REPO=1.0E0 C RETURN C 14 CONTINUE CC---IF IRP=6, KRG IS APPROX SG**N IF IRP=10 KRG IS APPROX 1-KRW C CCC---THREE PHASE RELATIVE K USING STONE'S 1970 METHOD MODIFIED BY CCC---AZIZ (1979) C RP1=RESID WATER SAT; C RP2=RESID OIL SAT; C RP3=RESID GAS SAT; C RP4=EXPONENT FOR KR. CCC---R. FALTA 6/29/89---------- NEXP=RP(4,NMAT) SWR=RP(1,NMAT) SGR=RP(3,NMAT) SNR=RP(2,NMAT) C C-----AQUEOUS PHASE RELATIVE PERMEABILITY IF (SL.GT.RP(1,NMAT)) THEN SS=(SL-RP(1,NMAT))/(1.0E0-RP(1,NMAT)) REPL=SS**NEXP ELSE REPL=0. ENDIF C C-----GAS PHASE RELATIVE PERMEABILITY REPG=0. IF(SG.LE.SGR) GOTO 301 IF (IRP(NMAT).EQ.10) THEN XSG=(1.0E0-SG-SWR)/(1.0E0-SWR) REPG=1.0E0-XSG**NEXP ELSE SSGR=(SG-SGR)/(1.0E0-SWR) REPG=(SSGR)**NEXP ENDIF IF(REPG.GT.1.) REPG=1. C C-----NAPL RELATIVE PERMEABILITY 301 CONTINUE REPO=0.0E0 IF(SO.LE.SNR) RETURN C SSUB=1.0E0-SWR-SNR SSO=(SO-SNR)/SSUB SSW=0.0E0 IF (SL.GT.SWR) SSW=(SL-SWR)/SSUB SSG=SG/SSUB C RKOCW=(1.0E0-SWR )**NEXP RKOW=(1.0E0-SL)**NEXP CCC RKOG=(1.0E0-SG-SWR )**NEXP C.....6-4-93: MODIFIED RKOG, SO THAT RKOG --> 0 FOR SG --> 1-SWR-SNR. RKOG=MAX(0.,1.-SG-SWR-SNR) RKOG=RKOG**NEXP C.....END OF MODIFICATION. C BW=RKOW/(RKOCW*(1.0E0-SSW)) BG=RKOG/(RKOCW*(1.0E0-SSG)) REPO=RKOCW*SSO*BW*BG IF (SO.GE.SNR .AND. SO.LT.(SNR+.005)) THEN FACT=((SO-SNR)/.005) REPO=REPO*FACT END IF RETURN C....................................................................... SUBROUTINE PCAP(T,SL,SG,DW,PC,PCO,NMAT) C C-----THIS ROUTINE COMPUTES CAPILLARY PRESSURE AS FUNCTION OF LIQUID C SATURATION SL AND TEMPERATURE T. C C------THREE PHASE CAPILLARY PRESSURES FROM PARKER ET AL., 1987 ---- C---------FALTA 5/89-----------ICP(NMAT)=8 -- C COMMON/RPCAP/IRP(27),RP(7,27),ICP(27),CP(7,27),IRPD,RPD(7), XICPD,CPD(7) C SAVE ICALL DATA ICALL/0/ ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'PCAP 1.1 ! 1 APRIL 1993',6X, X'CAPILLARY PRESSURE AS FUNCTION OF SATURATION (INCL 3-PHASE)') IF(ICP(NMAT).NE.8) PCO=0.0E0 SO=1.0E0-SL-SG GOTO(10,11,12,13,14,15,16,17,18),ICP(NMAT) C C C NO CAPILLARY PRESSURE C 18 CONTINUE PC=0.0E0 pco=0. RETURN C 10 CONTINUE C-----LINEAR FUNCTION. SLO=SL IF(CP(4,NMAT).NE.0.0E0) SLO=SO PC=-CP(1,NMAT)*(CP(3,NMAT)-SLO)/(CP(3,NMAT)-CP(2,NMAT)) IF(SLO.GE.CP(3,NMAT)) PC=0.0E0 IF(SLO.LE.CP(2,NMAT)) PC=-CP(1,NMAT) IF(CP(4,NMAT).NE.0.0E0) PCO=PC IF(CP(4,NMAT).NE.0.0E0) PC=0.0E0 RETURN C CC----THREE PHASE CAPILLARY PRESSURE MODIFIED FROM PARKER ET AL METHOD-- 17 CONTINUE SM=CP(1,NMAT) AN=CP(2,NMAT) ALAO=CP(3,NMAT) ALOW=CP(4,NMAT) SWBAR=(SL+0.0001E0-SM)/(1.0E0-SM) IF(SWBAR.GT.1.0E0) SWBAR=1.0E0 IF(SWBAR.LT.0.0E0) SWBAR=0.0E0 STBAR=(SL+0.0001E0+SO-SM)/(1.0E0-SM) IF(STBAR.GT.1.0E0) STBAR=1.0E0 IF(STBAR.LT.0.0E0) STBAR=0.0E0 AMINV=1.0E0/(1.0E0-1.0E0/AN) ANINV=1.0E0/AN DGAO=-DW*9.806E0/ALAO DGOW=-DW*9.806E0/ALOW IF(SWBAR.LT.0.1E0) GO TO 700 PCOW=0.0E0 IF(SWBAR.LT.1.0E0) PCOW=DGOW*(SWBAR**(-AMINV)-1.0E0)**ANINV GO TO 710 700 PCOW1=DGOW*(0.1E0**(-AMINV)-1.0E0)**ANINV PCOW2=DGOW*(0.101E0**(-AMINV)-1.0E0)**ANINV PCOW3=DGOW*(0.099E0**(-AMINV)-1.0E0)**ANINV SLOPE=(PCOW2-PCOW3)/0.002E0 DELT=0.1E0-SWBAR PCOW=PCOW1-SLOPE*DELT 710 CONTINUE IF(STBAR.LT.0.1E0) GO TO 720 PCO=0.0E0 IF(STBAR.LT.1.0E0) PCO=DGAO*(STBAR**(-AMINV)-1.0E0)**ANINV GO TO 730 720 PCO1=DGAO*(0.1E0**(-AMINV)-1.0E0)**ANINV PCO2=DGAO*(0.101E0**(-AMINV)-1.0E0)**ANINV PCO3=DGAO*(0.099E0**(-AMINV)-1.0E0)**ANINV SLOPE=(PCO2-PCO3)/0.002E0 DELT=0.1E0-STBAR PCO=PCO1-SLOPE*DELT 730 PC=PCO+PCOW RETURN C 11 CONTINUE C-----CAPILLARY PRESSURE FUNCTION OF PICKENS ET AL, AS GIVEN IN C J. HYDROLOGY 40, 243-264, 1979. SLX=MAX(SL,1.001E0*CP(2,NMAT)) IF(SLX.GT.0.999E0*CP(3,NMAT)) SLX=0.999E0*CP(3,NMAT) A=(1.0E0+SLX/CP(3,NMAT))*(CP(3,NMAT)-CP(2,NMAT))/ X(CP(3,NMAT)+CP(2,NMAT)) B=(1.0E0-SLX/CP(3,NMAT)) PC=-CP(1,NMAT)*LOG(A*(1.0E0+SQRT(1.0E0-B*B/(A*A)))/B)** X(1.0E0/CP(4,NMAT)) IF(SL.GT.0.999E0*CP(3,NMAT)) PC=PC*(1.0E0-SL)/0.001E0 RETURN C C 12 CONTINUE C-----CAPILLARY PRESSURE FUNCTION AS USED IN THE TRUST-PROGRAM, WHICH C WAS DEVELOPED BY T.N. NARASIMHAN AT LAWRENCE BERKELEY LABORATORY. C IF(SL.NE.1.0E0) GOTO 120 PC=0.0E0 RETURN C 120 SLX=SL IF(CP(5,NMAT).EQ.0.0E0)SLX=MAX(SL,1.001E0*CP(2,NMAT)) PC=-ABS(CP(5,NMAT)) IF(SLX.GT.CP(2,NMAT)) XPC=-CP(4,NMAT)-CP(1,NMAT)*((1.0E0-SLX)/(SLX-CP(2,NMAT))) X**(1.0E0/CP(3,NMAT)) IF(CP(5,NMAT).NE.0.0E0)PC=MAX(PC,-ABS(CP(5,NMAT))) IF(SL.GT.0.999E0) PC=PC*(1.0E0-SL)/0.001E0 RETURN C 13 CONTINUE C-----CAPILLARY PRESSURE OF YOLO CLAY AFTER CHRIS MILLY, C WATER RES. RES., VOL. 18 NO.3 (JUNE 1982), PP. 489-498. C IF(SL-CP(1,NMAT).GE.0.371E0) GOTO 130 SLX=MAX(SL,1.001E0*CP(1,NMAT)) EX=(0.371E0/(SLX-CP(1,NMAT))-1.0E0)**0.25E0 EX=2.26E0*EX-2.0E0 PC=-9.7783E3*10.0E0**EX RETURN C 130 PC=-97.783E0 RETURN 14 CONTINUE 15 CONTINUE C-----LEVERETT'S J-FUNCTION. SS=0.0E0 SLO=SL IF(CP(3,NMAT).NE.0.0E0) SLO=SO IF(SLO.GT.CP(2,NMAT)) SS=(SLO-CP(2,NMAT))/(1.0E0-CP(2,NMAT)) OSS=1.0E0-SS F=1.417E0*OSS-2.120E0*OSS**2+1.263E0*OSS**3 CALL SIGMA(T,ST) PC=-CP(1,NMAT)*ST*F IF(CP(3,NMAT).NE.0.0E0) PCO=PC IF(CP(3,NMAT).NE.0.0E0) PC=0.0E0 RETURN 16 CONTINUE C-----CAPILLARY FUNCTION OF VAN GENUCHTEN, SOIL SCI. SOC. AM. J. 44, C PP.892-898, 1980. C IF(SL.NE.1.0E0)GO TO 160 PC=0.0E0 RETURN C 160 SLX=SL IF(SLX.GE.CP(5,NMAT)) GOTO 161 IF(CP(4,NMAT).EQ.0.0E0)SLX=MAX(SL,1.001E0*CP(2,NMAT)) PC=-ABS(CP(4,NMAT)) IF(SLX.GT.CP(2,NMAT)) XPC=-1.0E0/ABS(CP(3,NMAT))*(((SL-CP(2,NMAT))/ #(CP(5,NMAT)-CP(2,NMAT))) X**(-1.0E0/CP(1,NMAT))-1.0E0)**(1.0E0-CP(1,NMAT)) IF(CP(4,NMAT).NE.0.0E0) PC=MAX(PC,-ABS(CP(4,NMAT))) IF(SL.GT.0.999E0) PC=PC*(1.0E0-SL)/0.001E0 RETURN 161 PC=0.0E0 RETURN C END