C FINITE DIFFERENCE SOLN TO 1-D, 2-PHASE FLOW; IMPES FORMULATION C INP 5 = INPUT FROM KEYBOARD; 10 = INPUT FROM 'INPUT.DAT' C RAT END POINT MOBILITY RATIO, INVADING/DISPLACED C KRD RELATIVE PERMEABILITY OF DISPLACED PHASE C KRI RELATIVE PERMEABILITY OF INVADING PHASE C EX RELATIVE PERMEABILITY EXPONENT C PCNO NPc C C1 C1 or entry J function in Tomeer model, negative for imbibition C PCA C2 IN THOMEER MODEL C SIC INITIAL SATURATION C ICC O: FLOW OF BOTH PHASES FROM LEFT TO RIGHT C 1: ALLOW COUNTER-CURRENT FLOW OF INVADING PHASE C 2: NO INVASION AT X=0; ONLY SPONTANUOUS, COUNTER-CURRENT FLOW C DTDX DT/DX C DSMAX IF NONZERO, MAXIMUM SAT. CHANGE PER TIME STEP C NX NUMBER OF GRID BLOCKS C NOUT NUMBER OF TIME STEPS BETWEEN PROFILES OUTPUT C TSTOP TIME TO STOP C NSTOP NUMBER OF TIME STEPS TO END C NF NUMBER OF TIME STEPS BETWEEN HISTORY OUTPUT C INFLOW B.C. UNIT FLUX OF INVADING PHASE OR ZERO FLUX IF ICC=2 C OUTFLOW B.C. POT(NX+1/2)=0.0 ,PC(NX+1/2)=0.0 C DIMENSION S(101),DS(101),PC(101),POT(101),TOTP(101),TMOB(101) " ,KRD(101),KRI(101) " ,COEFA(101),COEFB(101),COEFC(101),RIGHT(101) COMMON RAT,EXD,EXI,PCNO,PCA,C1 REAL KRD,KRI,KRDF,KRIF OPEN(10,FILE='INPUT.DAT') OPEN(11,FILE='PROD1.DAT') OPEN(12,FILE='PROD2.DAT') OPEN(13,FILE='PROD3.DAT') OPEN(14,FILE='PROD4.DAT') OPEN(15,FILE='PROD5.DAT') OPEN(16,FILE='PROD6.DAT') OPEN(17,FILE='TOTP.DAT') OPEN(18,FILE='POT.DAT') OPEN(21,FILE='PROF1.DAT') OPEN(22,FILE='PROF2.DAT') OPEN(23,FILE='PROF3.DAT') OPEN(24,FILE='PROF4.DAT') OPEN(25,FILE='PROF5.DAT') OPEN(26,FILE='PROF6.DAT') INP=10 ! 5 for keyboard; 10 for INPUT.DAT file IFILE=10 1 IFILE=IFILE+1 IFILE2=IFILE+10 WRITE(6,2) 2 FORMAT(/) WRITE(6,3) 3 FORMAT ('?RAT,EX,PCNO,C1,PCA,SIC,ICC = ') READ(INP,*) RAT,EX,PCNO,C1,PCA,SIC,ICC EXD=EX EXI=EX DSMAX=0. WRITE(6,4) 4 FORMAT ('?DTDX,NX = ') READ(INP,*) DTDX,NX DO 5 I=1,NX 5 S(I)=SIC POT(NX)=0.0 IF(C1.LT.0.0) POT(NX)=-1.E-6 DX=1.0/FLOAT(NX) DT=DTDX*DX NN=0 NNP=0 CUM=0. TIM=0. 26 CONTINUE WRITE(6,6) 6 FORMAT ('?TSTOP,NOUT,NF = ') READ(INP,*) TSTOP,NOUT,NF NSTOP=TSTOP*NX/DTDX +0.01 C **************** DO 20 NT=1,NSTOP NN=NN+1 DTDX2=DTDX/DX C SINGLE POINT MOBILITY C UPSTREAM FOR I+1/2 IS ALWAYS AT I FOR DISPLACED PHASE C UPSTREAM FOR INVADING PHASE MUST BE DETERMINED C CALCULATE PRESSURE COEFA(1)=0.0 COEFC(NX)=0.0 POT(NX+1)=0.0 PC(NX+1)=0.0 BCFAC=1.0 DO 30 I=1,NX PC(I)=PCF(S(I)) KRD(I)=KRDF(S(I)) C FIND UPSTREAM DIRECTION FOR INVADING PHASE IF ICC=1 OR 2 IF(ICC.EQ.0) THEN ! Flow only from left to right KRI(I)=KRIF(S(I)) ELSE ! Look for upstream direction IF(POT(I).GE.POT(I+1)) THEN KRI(I)=KRIF(S(I)) ELSE IF(I.LT.NX) THEN KRI(I)=KRIF(S(I+1)) ELSE KRI(I)=1.0 ENDIF ENDIF ENDIF TMOB(I)=KRD(I)/RAT+KRI(I) IF(I.GT.1) COEFA(I)=TMOB(I-1) IF(I.EQ.NX) BCFAC=2.0 COEFC(I)=BCFAC*TMOB(I) COEFB(I)=-(COEFA(I)+COEFC(I)) 30 RIGHT(I)=0.0 C PREVENT BACKFLOW OF INVADING PHASE IF ICC=0 IF(ICC.EQ.0.AND.POT(NX).LT.0.0) THEN KRI(NX)=0.0 TMOB(NX)=KRD(NX)/RAT+KRI(NX) COEFC(NX)=2.0*TMOB(NX) COEFB(NX)=-(COEFA(NX)+COEFC(NX)) ENDIF IF(ICC.EQ.2)THEN FLUXIN=0. ! ZERO FLUX AT X=0 ELSE FLUXIN=1. ! UNIT FLUX AT X=0 ENDIF RIGHT(1)=-DX*FLUXIN+KRD(1)*(PC(2)-PC(1))/RAT BCFAC=1.0 DO 31 I=2,NX IF(I.EQ.NX) BCFAC=2.0 C INCLUDE CAPILLARY PRESSURE 31 RIGHT(I)=RIGHT(I) "+(KRD(I)*BCFAC*(PC(I+1)-PC(I))-KRD(I-1)*(PC(I)-PC(I-1)))/RAT CALL TRIDAG(COEFA,COEFB,COEFC,RIGHT,POT,NX) C CALCULATE SATURATION DS(1)=DTDX2*(KRI(1)*(POT(2)-POT(1)) +DX*FLUXIN) BCFAC=1.0 DO 40 I=2,NX IF(I.EQ.NX) BCFAC=2.0 40 DS(I)=DTDX2*(KRI(I)*BCFAC*(POT(I+1)-POT(I)) " -KRI(I-1)*(POT(I)-POT(I-1))) SUM=0. DO 13 I=1,NX S(I)=S(I)+DS(I) TOTP(I)=POT(I)-PC(I) 13 SUM=SUM+S(I) EFF=SUM/FLOAT(NX)-SIC CUT=KRI(NX)*POT(NX)/(KRI(NX)*POT(NX)+KRD(NX)*TOTP(NX)/RAT) CUM=CUM+(1.0-CUT)*DT RATE=KRD(NX)*TOTP(NX)/(0.5*DX*RAT) TIM=TIM+DT NNP=NNP+1 C Print and file production IF(NNP.LT.NF) GO TO 22 WRITE(6,23) TIM,CUT,EFF,CUM,RATE WRITE(IFILE,23) TIM,CUT,EFF,CUM,RATE NNP=0 22 CONTINUE 23 FORMAT(1X,F10.6,4F10.4) C Print and file profile IF(NN.LT.NOUT) GO TO 20 NN=0 WRITE(6,14) TIM 14 FORMAT(/' TIM = ',F10.3) WRITE(6,15) (S(I),I=1,NX) 15 FORMAT(1X,10F6.3) DX=1.0/FLOAT(NX) DO 16 I=1,NX X=DX*(FLOAT(I)-0.5) WRITE(17,17) X,TOTP(I) WRITE(18,17) X,POT(I) 16 WRITE(IFILE2,17) X,S(I) 17 FORMAT(2F10.5) C C TIME STEP SIZE C IF(DSMAX.EQ.0) GO TO 19 C MAXIMUM SATURATION CHANGE PER TIME STEP DTSEL=0.1 IF(TIM.GT.0.1) DTSEL=0.50*TIM DSAT=0.0 DO 18 I=1,NX ABSAT=ABS(DS(I)) IF(ABSAT.GT.DSAT) DSAT=ABSAT 18 CONTINUE DTSEL=MIN(DTSEL,DT*DSMAX/DSAT) IF(NT.GT.1) THEN DTOLD=DT DT=0.5*(DT+DTSEL) DTDX=DT/DX ENDIF 19 CONTINUE C IF(TIM.GT.TSTOP) GO TO 24 20 CONTINUE C END OF LOOP ON TIME 24 WRITE(6,25) 25 FORMAT(1X'MORE ? '\) READ(5,*) MORE IF(MORE.EQ.1) GO TO 1 IF(MORE.EQ.2) GO TO 26 STOP END C SUBROUTINE FOR IMPES FUNCTION KRDF(S) COMMON RAT,EXD,EXI REAL KRDF IF(S.GE.1.0) THEN KRDF=0.0 RETURN ENDIF KRDF=(1.0-S)**EXD RETURN END FUNCTION KRIF(S) COMMON RAT,EXD,EXI REAL KRIF IF(S.LE.0.0) THEN KRIF=0.0 RETURN ENDIF KRIF=S**EXI RETURN END FUNCTION PCF(S) COMMON RAT,EXD,EXI,PCNO,PCA,C1 IF(C1.GE.0.0) THEN C DRAINAGE, for C1>0 IF(S.GT.1.E-4.AND.S.LT.0.999) THEN PCF=PCNO*C1*EXP(-PCA/ALOG(S)) ELSE IF(S.LT.1.E-4) THEN PCF=PCNO*C1 ELSE PCF=PCNO*C1*EXP(-PCA/ALOG(0.999)) ENDIF ENDIF ELSE C IMBIBITION, for C1<0 IF(S.GT.1.E-4.AND.S.LT.0.999) THEN PCF=PCNO*C1*EXP(-PCA/ALOG(1.0-S)) ELSE IF(S.GE.0.999) THEN PCF=PCNO*C1 ELSE PCF=PCNO*C1*EXP(-PCA/ALOG(0.999)) ENDIF ENDIF ENDIF RETURN END SUBROUTINE TRIDAG(A,B,C,R,U,N) PARAMETER (NMAX=100) DIMENSION GAM(NMAX),A(N),B(N),C(N),R(N),U(N) IF(B(1).EQ.0.)PAUSE BET=B(1) U(1)=R(1)/BET DO 11 J=2,N GAM(J)=C(J-1)/BET BET=B(J)-A(J)*GAM(J) IF(BET.EQ.0.)PAUSE U(J)=(R(J)-A(J)*U(J-1))/BET 11 CONTINUE DO 12 J=N-1,1,-1 U(J)=U(J)-GAM(J+1)*U(J+1) 12 CONTINUE RETURN END