PEGN 624 Project
S2009
Jeff Brown
2009-04-06

Definitions

Packages

Needs["PlotLegends`"]

Definitions

Initialization

Data

Grid Constants

Δx = 100 ;(*ft*)

Δy = 100; (*ft*)

Δz = 20;(*ft*)

Δt = 1; (* day *)

VR = Δx*Δy*Δz

VRdt = VR/Δt

Tolerances

ϵP = 1.* 10^(-3)
ϵSo = ϵSg = 1.* 10^(-4)
ϵTs = ϵT = 1.* 10^(-1)
ϵG = 1.*10^(-2)

Rock Properties

km = 3000;  (*md*)

ϕm=0.3;

convertk = 0.006328

convertq = 5.6146

convertT = 459.67

convertU=((144/778.26))

Molecular Weight

MWo=250;

MWg=MWw=18;

Thermal Properties

ρR = 170 (* lbm/ft^3 *)

CR = 0.206 (* btu/(lbm °F) *)

kT = 1.6 * 24 (* btu/(day ft °F) *)

TR=100;

Ts

Correlations: T in °F; P in psia

calcTs[P_]:=116(P)^0.224

calcTs[P_,T_]:=calcTs[P]

calcTs[100]

calcTs[600]

Viscosities

Correlations: T in °F

calcμo[T_]:=3.62431*10^(-4)Exp[8485.4/(T+convertT)]

calcμo[P_,T_]:=calcμo[T]

{calcμo[100],calcμo[400]}

calcμw[T_]:=(2.185/((0.04012 T + 5.1547*(10^((-6)))*(T^2)-1)))

calcμw[P_,T_]:=calcμw[T]

{calcμw[100],calcμw[400]}

calcμg[T_]:=8.822*10^(-6) (T+convertT)^1.116

calcμg[P_,T_]:=calcμg[T]

{calcμg[100],calcμg[400]}

Molar Density

Correlations: T in °F; P in psia

calcξo[P_,T_]:=1.683*(0.1987-7.5885*10^(-5)*(T+convertT))*Exp[10^(-5)*(P-14.7)]

{calcξo[100,100],calcξo[100,400],calcξo[600,100],calcξo[600,400]}

calcξw[P_,T_]:=22.14*(0.202-7.725*10^(-5)*(T+convertT))*Exp[4*10^(-6)*(P-118)]

{calcξw[100,100],calcξw[100,400],calcξw[600,100],calcξw[600,400]}

May need to use other correlation where P > 1000.

calcξg[P_]:=((P^0.9585)/(363.9*18))

calcξg[P_,T_]:=calcξg[P]

{calcξg[100],calcξg[600]}

Molar Density Derivatives

calcDξoDP[P_,T_]:=10^(-5)*calcξo[P,T]

{calcDξoDP[100,100],calcDξoDP[100,400],calcDξoDP[600,100],calcDξoDP[600,400]}

D[calcξo[P,T],T]

calcDξoDT[P_,T_]:=-0.00012771445500000002` (E^((-14.7`+P)/100000))

{calcDξoDT[100,100],calcDξoDT[100,400],calcDξoDT[600,100],calcDξoDT[600,400]}

calcDξwDP[P_,T_]:=4*10^(-6)*calcξw[P,T]

{calcDξwDP[100,100],calcDξwDP[100,400],calcDξwDP[600,100],calcDξwDP[600,400]}

D[calcξw[P,T],T]

calcDξwDT[P_,T_]:=-0.0017103150000000002` (E^((-118+P)/250000))

{calcDξwDT[100,100],calcDξwDT[100,400],calcDξwDT[600,100],calcDξwDT[600,400]}

calcDξgDP[P_]:=0.9585*((P^(0.9585-1))/(363.9*18))

calcDξgDP[P_,T_]:=calcDξgDP[P]

{calcDξgDP[100],calcDξgDP[600]}

calcDξgDT[P_,T_]:=0

Specific Gravities

convertdensity = 0.4335275 (*psi/ft->g/cm^3*)

baseρw=62.42796

calcγo[P_,T_]:=((MWo calcξo[P,T])/(baseρw))*convertdensity

{calcγo[100,100],calcγo[100,400],calcγo[600,100],calcγo[600,400]}

calcγw[P_,T_]:=((MWw calcξw[P,T])/(baseρw))*convertdensity

{calcγw[100,100],calcγw[100,400],calcγw[600,100],calcγw[600,400]}

calcγg[P_,T_]:=((MWw calcξg[P,T])/(baseρw))*convertdensity

{calcγg[100,100],calcγg[100,400],calcγg[600,100],calcγg[600,400]}

Heat Capacities

Correlations: T in °F; P in psia

calcCpo[T_]:=57.88+6.03*10^(-2)(T+convertT)

calcCpo[P_,T_]:=calcCpo[T]

{calcCpo[100],calcCpo[400]}

calcCpw[T_]:=41.809-7.8289*10^(-2)(T+convertT)+6.2788*10^(-5)(T+convertT)^2

calcCpw[P_,T_]:=calcCpw[T]

{calcCpw[100],calcCpw[400]}

calcCpg[T_]:=7.136+1.4667*10^(-3)(T+convertT)+2.55*10^(-8)(T+convertT)^2

calcCpg[P_,T_]:=calcCpg[T]

{calcCpg[100],calcCpg[400]}

Enthalpies

Correlations: T in °F; P in psia

Integrate[calcCpo[To],{To,32,77}]

Integrate[calcCpo[To],{To,77,T}]

calcHo[T_]:=3999.8002950000005`+-6769.813127000001`+85.59810100000001` T+0.030150000000000003` (T^2)

calcHo[P_,T_]:=calcHo[T]

Integrate[calcCpw[To],{To,32,T}]

calcHwless[T_]:=-600.9972648597292`+19.088780570813203` T-0.010282740039999995` (T^2)+0.000020929333333333336` (T^3)

calcHwless[P_,T_]:=calcHwless[T]

Integrate[calcCpw[To],{To,32,calcTs[P]}]

calcHwmore[P_]:=-600.9972648597292`+2214.2985462143315` (P^(0.224`))-138.36454997823992` (P^(0.448`))+32.66851268266667` (P^(0.672`))

calcHwmore[P_,T_]:=calcHwmore[P]

calcHgmore[P_,T_]:=calcHwmore[P]+1318(((T/116))^4.464)^-0.08774

Internal Energies

calcUo[P_,T_]:=calcHo[T] - convertU (P/(MWo calcξo[P,T]))

{calcUo[100,100],calcUo[100,400],calcUo[600,100],calcUo[600,400]}

calcUwless[P_,T_]:=calcHwless[T] - convertU (P/(MWw calcξw[P,T]))

{calcUwless[100,100],calcUwless[100,400],calcUwless[600,100],calcUwless[600,400]}

calcUwmore[P_,T_]:=calcHwmore[P] - convertU (P/(MWw calcξw[P,T]))

{calcUwmore[100,500],calcUwmore[100,1000],calcUwmore[600,500],calcUwmore[600,1000]}

calcUgmore[P_,T_]:=calcHgmore[P,T]-convertU(P/(MWg calcξg[P]))

{calcUgmore[100,500],calcUgmore[100,1000],calcUgmore[600,500],calcUgmore[600,1000]}

Internal Energies Derivatives

D[calcUo[P,T],T]

calcDUoDT[P_,T_]:=85.59810100000001`+0.060300000000000006` T-((3.337102883149041`*^-8 (E^((14.7`-P)/100000)) P)/(((0.1987`-0.00007588500000000001` (459.67`+T)))^2))

{calcDUoDT[100,100],calcDUoDT[100,400],calcDUoDT[100,1000],calcDUoDT[600,100],calcDUoDT[600,400],calcDUoDT[600,1000]}

D[calcUo[P,T],P]

calcDUoDP[P_,T_]:=-((0.0004397579077747962` (E^((14.7`-P)/100000)))/(0.1987`-0.00007588500000000001` (459.67`+T)))+((4.397579077747962`*^-9 (E^((14.7`-P)/100000)) P)/(0.1987`-0.00007588500000000001` (459.67`+T)))

{calcDUoDP[100,100],calcDUoDP[100,400],calcDUoDP[100,1000],calcDUoDP[600,100],calcDUoDP[600,400],calcDUoDP[600,1000]}

D[calcUwless[P,T],T]

calcDUwlessDT[P_,T_]:=85.59810100000001`+0.060300000000000006` T-((3.337102883149041`*^-8 (E^((14.7`-P)/100000)) P)/(((0.1987`-0.00007588500000000001` (459.67`+T)))^2))

{calcDUwlessDT[100,100],calcDUwlessDT[100,400],calcDUwlessDT[600,100],calcDUwlessDT[600,400]}

D[calcUwless[P,T],P]

calcDUwlessDP[P_,T_]:=-((0.00046428821563847613` (E^((118-P)/250000)))/(0.202`-0.00007725000000000001` (459.67`+T)))+((1.8571528625539045`*^-9 (E^((118-P)/250000)) P)/(0.202`-0.00007725000000000001` (459.67`+T)))

{calcDUwlessDP[100,100],calcDUwlessDP[100,400],calcDUwlessDP[600,100],calcDUwlessDP[600,400]}

D[calcUwmore[P,T],T]

calcDUwmoreDT[P_,T_]:=-((3.586626465807229`*^-8 (E^((118-P)/250000)) P)/(((0.202`-0.00007725000000000001` (459.67`+T)))^2))

{calcDUwmoreDT[100,500],calcDUwmoreDT[100,1000],calcDUwmoreDT[600,500],calcDUwmoreDT[600,1000]}

D[calcUwmore[P,T],P]

calcDUwmoreDP[P_,T_]:=((496.0028743520103`)/(P^(0.776`)))-((61.98731839025149`)/(P^(0.552`)))+((21.953240522752004`)/(P^(0.32799999999999996`)))-((0.00046428821563847613` (E^((118-P)/250000)))/(0.202`-0.00007725000000000001` (459.67`+T)))+((1.8571528625539045`*^-9 (E^((118-P)/250000)) P)/(0.202`-0.00007725000000000001` (459.67`+T)))

{calcDUwmoreDP[100,500],calcDUwmoreDP[100,1000],calcDUwmoreDP[600,500],calcDUwmoreDP[600,1000]}

D[calcUgmore[P,T],T]

calcDUgmoreDT[P_,T_]:=-((3322.203982900616` (T^(3.4640000000000004`)))/((((T^(4.464`))))^(1.08774`)))

{calcDUgmoreDT[100,500],calcDUgmoreDT[100,1000],calcDUgmoreDT[600,500],calcDUgmoreDT[600,1000]}

D[calcUgmore[P,T],P]

calcDUgmoreDP[P_,T_]:=-((2.7942672114717433`)/(P^(0.9585`)))+((496.0028743520103`)/(P^(0.776`)))-((61.98731839025149`)/(P^(0.552`)))+((21.953240522752004`)/(P^(0.32799999999999996`)))

{calcDUgmoreDP[100,500],calcDUgmoreDP[100,1000],calcDUgmoreDP[600,500],calcDUgmoreDP[600,1000]}

Compressibilities

cϕ=10*10^(-6); (* psi^-1*)

Source Terms

flowsource = 35000./18

heatsource = 1000.*flowsource

Matrix Initialization

Mzero=Table[Table[0,{4}],{4}]

IMAX=3

A11=Mzero

A21=Mzero;
A12=

Capillary Pressure, Rel Perm

Rel Perm

Saturations

Swrm=0.20;
Sowrm=0.30;
Sogrm = 0.20;
Sgcm = 0.05;

Rel Perm

krwmstar=0.2;
nwm=3.0;

krowmstar=0.6;
nowm=4.0;

krgmstar = 0.5;
ngm=3.0;

computekrwm[Sw_]:=
Which[Sw≤Swrm,0,
Sw≥(1-Sowrm),krwmstar,
True,
krwmstar*((((Sw-Swrm))/((1-Swrm-Sowrm))))^nwm
]

computekrowm[So_]:=
Which[So≤Sowrm,0,
So≥(1-Swrm),krowmstar,
True,
krowmstar*((((So-Sowrm))/((1-Swrm-Sowrm))))^nowm
]

computekrgm[Sg_]:=
Which[Sg≤Sgcm,0,
Sg≥(1-Swrm-Sogrm),krgmstar,
True,
krgmstar*((((Sg-Sgcm))/((1-Swrm-Sogrm-Sgcm))))^ngm
]

Plot[{computekrwm[Sw],computekrowm[1-Sw]},{Sw,0,1},Frame->True,FrameLabel->{"(Subscript[S, w])","(Subscript[k, r])"},PlotLabel->"Rel Perm, Oil-Water System, Matrix",
PlotStyle->{Blue,Green},ImageSize->600,PlotLegend->{"(Subscript[k, rwm])", "(Subscript[k, rowm])"},LegendPosition->{.9,0},LegendSize->{.5,.3}]

Plot[{computekrgm[Sg],computekrowm[1-Sg-Swrm]},{Sg,0,1},Frame->True,FrameLabel->{"(Subscript[S, g])","(Subscript[k, r])"},PlotLabel->"Rel Perm, Gas-Oil System at Swr, Matrix",
PlotStyle->{Blue,Green},ImageSize->600,PlotLegend->{"(Subscript[k, rgm])","(Subscript[k, rom])"},LegendPosition->{.9,0},LegendSize->{.5,.3}]

Plot[{computekrgm[Sg],computekrwm[1-Sg-Sowrm]},{Sg,0,1},Frame->True,FrameLabel->{"(Subscript[S, g])","(Subscript[k, r])"},PlotLabel->"Rel Perm, Gas-Oil System at Sowrm, Matrix",
PlotStyle->{Blue,Green},ImageSize->600,PlotLegend->{"(Subscript[k, rgm])","(Subscript[k, rwm])"},LegendPosition->{.9,0},LegendSize->{.5,.3}]

Capillary Pressure

Capillary Pressure

Pcowstar = 5;
Pcogstar = 10;
Pcogentry=3;
ncw = 1.76;
ncg = 1.20;

Pcogm1[Sg_]:=Pcogentry+Pcogstar*((((Sg-Sgcm))/((1-Swrm-Sogrm-Sgcm))))^ncg

Pcowm1[Sw_]:=Pcowstar*((((1-Sowrm-Sw))/((1-Sowrm-Swrm))))^ncw

computePcogm[Sg_]:=
Module[{},
Which[Sg≤Sgcm,Pcogentry,
Sg≥(1-Swrm-Sogrm),Pcogentry+Pcogstar,
True,Pcogm1[Sg]
]]

computePcowm[Sw_]:=
Module[{},
Which[Sw≤Swrm,Pcowstar,
Sw≥(1-Sowrm),0,
True,Pcowm1[Sw]
]]

Plot[computePcogm[Sg],{Sg,0,1},Frame->True,FrameLabel->{"(Subscript[S, g])","(Subscript[P, cgom])"},PlotLabel->"Capillary Pressure, Gas-Oil System with Residual Water, Matrix",
ImageSize->500,PlotRange->All,AxesOrigin->{0,0}]

Plot[computePcowm[Sw],{Sw,0,1},Frame->True,FrameLabel->{"(Subscript[S, w])","(Subscript[P, cwom])"},PlotLabel->"Capillary Pressure, Water-Oil System with Residual Water, Matrix",
ImageSize->500,PlotRange->All,AxesOrigin->{0,0}]

Transmissibilities & Fractional Flows

Transmissibilities

use upstream only here since k, Δx are constant

calcTw[P_,T_,So_,Sg_,Sw_]:=
convertk*km*((computekrwm[Sw])/(calcμw[P,T]))(((Δy Δz))/(Δx))

calcTo[P_,T_,So_,Sg_,Sw_]:=
convertk*km*((computekrowm[So])/(calcμo[P,T]))(((Δy Δz))/(Δx))

calcTg[P_,T_,So_,Sg_,Sw_]:=
convertk*km*((computekrgm[Sg])/(calcμg[P,T]))(((Δy Δz))/(Δx))

a1={500,100,0.8,0,0.2}

{calcTw@@a1,calcTo@@a1,calcTg@@a1}

a2={500,500,0.7,0.1,0.2}

{calcTw@@a2,calcTo@@a2,calcTg@@a2}

Thermal Transmissibility

TE = kT(((Δy Δz))/(Δx))

Thermal Boundaries

calcqb[T_]:=-kT*Δx*Δy*(((T-TR))/(Sqrt[Pi ((kT/(ρR*CR))) Δt]))

TE = kT(((Δy Δz))/(Δx))

Fractional Flow

use upstream only here since k, Δx are constant

calcλw[P_,T_,So_,Sg_,Sw_]:=
((computekrwm[Sw])/(calcμw[P,T]))

calcλo[P_,T_,So_,Sg_,Sw_]:=
((computekrowm[So])/(calcμo[P,T]))

calcλg[P_,T_,So_,Sg_,Sw_]:=
((computekrgm[Sg])/(calcμg[P,T]))

calcλtot[P_,T_,So_,Sg_,Sw_]:=
calcλw[P,T,So,Sg,Sw]+calcλo[P,T,So,Sg,Sw]+calcλg[P,T,So,Sg,Sw]

calcffw[P_,T_,So_,Sg_,Sw_]:=
calcλw[P,T,So,Sg,Sw]/calcλtot[P,T,So,Sg,Sw]

calcffo[P_,T_,So_,Sg_,Sw_]:=
calcλo[P,T,So,Sg,Sw]/calcλtot[P,T,So,Sg,Sw]

calcffg[P_,T_,So_,Sg_,Sw_]:=
calcλg[P,T,So,Sg,Sw]/calcλtot[P,T,So,Sg,Sw]

a1={500,100,0.8,0,0.2}

{calcffw@@a1,calcffo@@a1,calcffg@@a1}

a2={500,500,0.7,0.1,0.2}

{calcffw@@a2,calcffo@@a2,calcffg@@a2}

Thermodynamic Constraints

calcG[P_,T_,So_,Sg_,Sw_]:=
((((T/116)))^4.464)*((Sw/(Sw+ϵG)))-P

calcDGDP[P_,T_,So_,Sg_,Sw_]:=
-1

calcDGDT[P_,T_,So_,Sg_,Sw_]:=
(4.464/116)((((T/116)))^3.464)*((Sw/(Sw+ϵG)))

calcDGDSo[P_,T_,So_,Sg_,Sw_]:=
((((T/116)))^4.464)*(((-ϵG)/(((Sw+ϵG))^2)))

At T=100F, P=500; n, ℓ

Conditions

a1n={P1n=500 (*P*),T1n=100 (*T*),So1n=0.8 (*So*),Sg1n=0 (*Sg*),Sw1n=0.2 (*Sw*)}

a11n=Take[a1n,2]

a1ℓ={P1ℓ=500 (*P*),T1ℓ=100 (*T*),So1ℓ=0.8 (*So*),Sg1ℓ=0 (*Sg*),Sw1ℓ=0.2 (*Sw*)}

a11ℓ=Take[a1ℓ,2]

T1ℓ

Ts1ℓ=calcTs[P1ℓ,T1ℓ]

Thus T1 < Ts

va1n={a1n,a1n,a1n}

Evaluate for water at node 1

ξtot = calcξw@@a11n;

qtot=flowsource/ξtot

Transmissibility

Tw1=Apply[calcTw,va1n,{1}]

Tg1=Apply[calcTg,va1n,{1}]

To1=Apply[calcTo,va1n,{1}]

Water Equation

Term 1,9,11  (δP)

W1=(calcTw@@a1n)*(calcξw@@a11n)

W9=-W1

W11=W1

Term 2,10,12  (δP)

W2=(calcTg@@a1n)*(calcξg@@a11n)

W10=-W2

W12=-W2

Term 3 (δP)

W3=-VRdt*ϕm*cϕ*Sw1ℓ*calcξw[P1ℓ,T1ℓ]+VRdt*ϕm*Sw1ℓ*calcDξwDP[P1ℓ,T1ℓ]

Term 4 (δT)

W4=-VRdt*ϕm*Sw1ℓ*calcDξwDT[P1ℓ,T1ℓ]

Term 5 (δSo)

W5=VRdt*ϕm*calcξw[P1ℓ,T1ℓ]

W6=0

Term 7 (ℓ)

W7=VRdt*(ϕm*Sw1ℓ*calcξw[P1ℓ,T1ℓ])

Term 8 (n)

W8=-VRdt*(ϕm*Sw1n*calcξw[P1n,T1n])

Term 13, 14 :  0 since ΔD = 0

W13=W14=0

Term 15, cell 1 injection

W15a=-flowsource

Term 15, cell N production

W15c=-(-qtot*(calcffw@@a1n)*(calcξw@@a11n)-qtot*(calcffg@@a1n)*(calcξg@@a11n))

Combine Water Equation

A11[[1,1]]=-2W1+-2W2

A21

WLLleft={W1+W2,0,0,0}

WLLright={W1+W2,0,0,0}

WLLcenter={-2W1+-2W2,0,0,0}
WRLcenter={W3,W4,W5,W6}

WRRcenter = W7+W8

WLRcenter = W9

Oil Equation

Term 1,7  (δP)

O1=(calcTo@@a1n)*(calcξo@@a11n)

O7=-O1

Term 2 (δP)

O2=-VRdt*ϕm*So1ℓ*(cϕ*calcξo[P1ℓ,T1ℓ]+calcDξoDP[P1ℓ,T1ℓ])

Term 3 (δT)

O3=-VRdt*ϕm*So1ℓ*calcDξoDT[P1ℓ,T1ℓ]

Term 4 (δSo)

O4=-VRdt*ϕm*calcξo[P1ℓ,T1ℓ]

Term 5 (ℓ)

O5=VRdt*(ϕm*So1ℓ*calcξo[P1ℓ,T1ℓ])

Term 6 (n)

O6=-VRdt*(ϕm*So1n*calcξo[P1n,T1n])

Term 8 :  0 since ΔD = 0

O8=0

Term 9, cell N production

O9c=-(-qtot*(calcffo@@a1n)*(calcξo@@a11n))

Energy Equation, LHS

Term 1  (δP)

E1=(calcTw@@a1n)*(calcξw@@a11n)*(calcHwless@@a11n)

E19=-E1

E23=E1

Term 2  (δP)

E2=(calcTg@@a1n)*(calcξg@@a11n)*(calcHgmore@@a11n)

E20=-E2

E24=-E2

Term 3  (δP)

E3=(calcTo@@a1n)*(calcξo@@a11n)*(calcHo@@a11n)

E21=-E3

Term 4  (δT)

E4=TE

E22=-E4

Term 5 (δP)

E5=-VRdt*ϕm*Sg1ℓ*(cϕ*calcξg[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcDξgDP[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcξg[P1ℓ,T1ℓ]*(calcDUgmoreDP@@a11ℓ))

Term 6 (δP)

E6=-VRdt*ϕm*Sw1ℓ*(cϕ*calcξw[P1ℓ,T1ℓ]*(calcUwless@@a11ℓ)+calcDξwDP[P1ℓ,T1ℓ]*(calcUwless@@a11ℓ)+calcξw[P1ℓ,T1ℓ]*(calcDUwlessDP@@a11ℓ))

Term 7 (δP)

E7=-VRdt*ϕm*So1ℓ*(cϕ*calcξo[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcDξoDP[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcξo[P1ℓ,T1ℓ]*(calcDUoDP@@a11ℓ))

Term 8 (δP)

E8=VRdt*ϕm*cϕ*ρR*CR*T1ℓ

Term 9 (δT)

E9=-VRdt*ϕm*Sg1ℓ*(calcDξgDT[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcξg[P1ℓ,T1ℓ]*(calcDUgmoreDT@@a11ℓ))

Term 10 (δT)

E10=-VRdt*ϕm*Sw1ℓ*(calcDξwDT[P1ℓ,T1ℓ]*(calcUwless@@a11ℓ)+calcξw[P1ℓ,T1ℓ]*(calcDUwlessDT@@a11ℓ))

Term 11 (δT)

E11=-VRdt*ϕm*So1ℓ*(calcDξoDT[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcξo[P1ℓ,T1ℓ]*(calcDUoDT@@a11ℓ))

Term 12 (δT)

E12=-VRdt*(1-ϕm)*ρR*CR

Term 13 (δSg)

E13=VRdt*ϕm*(calcξw[P1ℓ,T1ℓ]calcUwless[P1ℓ,T1ℓ]-calcξg[P1ℓ,T1ℓ]calcUgmore[P1ℓ,T1ℓ])

Term 14 (δSo)

E14=VRdt*ϕm*(calcξw[P1ℓ,T1ℓ]calcUwless[P1ℓ,T1ℓ]-calcξo[P1ℓ,T1ℓ]calcUo[P1ℓ,T1ℓ])

Energy Equation, RHS

Term 15 (ℓ)

E15 = VRdt*ϕm*(Sw1ℓ*calcξw[P1ℓ,T1ℓ]*calcUwless[P1ℓ,T1ℓ]+
So1ℓ*calcξo[P1ℓ,T1ℓ]*calcUo[P1ℓ,T1ℓ]+
Sg1ℓ*calcξg[P1ℓ,T1ℓ]*calcUgmore[P1ℓ,T1ℓ])

Term 16 (n)

E16 = -VRdt*ϕm*(Sw1n*calcξw[P1n,T1n]*calcUwless[P1n,T1n]+
So1n*calcξo[P1n,T1n]*calcUo[P1n,T1n]+
Sg1n*calcξg[P1n,T1n]*calcUgmore[P1n,T1n])

Term 17 (ℓ)

E17=VRdt*(1-ϕm)*ρR*CR*T1ℓ

Term 18 (n)

E18=-VRdt*(1-ϕm)*ρR*CR*T1n

Term 25,26,27 =0 since ΔD=0

E25=E26=E27=0

Term 28-29, cell 1 injection

E28a=-heatsource

Term 28-29, cell N production

E28c=-(-qtot)*((calcffw@@a1n)*calcξw[P1n,T1n]*calcHwless[P1n,T1n]+
(calcffo@@a1n)*calcξo[P1n,T1n]*calcHo[P1n,T1n]+
(calcffg@@a1n)*calcξg[P1n,T1n]*calcHgmore[P1n,T1n])

Term 30, to boundaries

E30=-(calcqb[T1n]+calcqb[T1n])

Thermodynamic Constraints

Term 1; δP

G1=-1;

Term 2; δSo

G2=calcDGDSo[P1ℓ,Ts1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]

Term 3; δTs

G3=calcDGDT[P1ℓ,Ts1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]

Term 4; (ℓ)

G4=-(calcG[P1ℓ,Ts1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]+computePcowm[Sw1n])

Assemble Equations

At T=500F, P=500; n, ℓ

Conditions

a1n={P1n=500 (*P*),T1n=500 (*T*),So1n=0.7 (*So*),Sg1n=0.1 (*Sg*),Sw1n=0.2 (*Sw*)}

a11n=Take[a1n,2]

a1ℓ={P1ℓ=500 (*P*),T1ℓ=500 (*T*),So1ℓ=0.7 (*So*),Sg1ℓ=0.1 (*Sg*),Sw1ℓ=0.2 (*Sw*)}

a11ℓ=Take[a1ℓ,2]

T1ℓ

Ts1ℓ=calcTs[P1ℓ,T1ℓ]

Thus T1 >=Ts

Water Equation

Term 1,9,11  (δP)

W1=(calcTw@@a1n)*(calcξw@@a11n)

W9=-W1

W11=W1

Term 2,10,12  (δP)

W2=(calcTg@@a1n)*(calcξg@@a11n)

W10=-W2

W12=-W2

Term 3 (δP)

W3=-(VRdt*ϕm*cϕ*Sw1ℓ*calcξw[P1ℓ,T1ℓ]+VRdt*ϕm*Sw1ℓ*calcDξwDP[P1ℓ,T1ℓ]+
VRdt*ϕm*cϕ*Sg1ℓ*calcξg[P1ℓ,T1ℓ]+VRdt*ϕm*Sg1ℓ*calcDξgDP[P1ℓ,T1ℓ])

Term 4 (δT)

W4=-(VRdt*ϕm*Sw1ℓ*calcDξwDT[P1ℓ,T1ℓ]+VRdt*ϕm*Sg1ℓ*calcDξgDT[P1ℓ,T1ℓ])

Term 5 (δSo)

W5=VRdt*ϕm*calcξw[P1ℓ,T1ℓ]

W6=VRdt*ϕm*(calcξw[P1ℓ,T1ℓ]-calcξg[P1ℓ,T1ℓ])

Term 7 (ℓ)

W7=VRdt*ϕm(Sw1ℓ*calcξw[P1ℓ,T1ℓ]+Sg1ℓ*calcξg[P1ℓ,T1ℓ])

Term 8 (n)

W8=-VRdt*ϕm(Sw1n*calcξw[P1n,T1n]+Sg1n*calcξg[P1n,T1n])

Term 13, 14 :  0 since ΔD = 0

W13=W14=0

Term 15, cell 1 injection

W15a=-flowsource

Term 15, cell N production

W15c=-(-flowsource*(calcffw@@a1n)*(calcξw@@a11n)-flowsource*(calcffg@@a1n)*(calcξg@@a11n))

Oil Equation

Term 1,7  (δP)

O1=(calcTo@@a1n)*(calcξo@@a11n)

O7=-O1

Term 2 (δP)

O2=-VRdt*ϕm*So1ℓ*(cϕ*calcξo[P1ℓ,T1ℓ]+calcDξoDP[P1ℓ,T1ℓ])

Term 3 (δT)

O3=-VRdt*ϕm*So1ℓ*calcDξoDT[P1ℓ,T1ℓ]

Term 4 (δSo)

O4=-VRdt*ϕm*calcξo[P1ℓ,T1ℓ]

Term 5 (ℓ)

O5=VRdt*(ϕm*So1ℓ*calcξo[P1ℓ,T1ℓ])

Term 6 (n)

O6=-VRdt*(ϕm*So1n*calcξo[P1n,T1n])

Term 8 :  0 since ΔD = 0

O8=0

Term 9, cell N production

O9c=-(-flowsource*(calcffo@@a1n)*(calcξo@@a11n))

Energy Equation, LHS

Term 1  (δP)

E1=(calcTw@@a1n)*(calcξw@@a11n)*(calcHwless@@a11n)

E19=-E1

E23=E1

Term 2  (δP)

E2=(calcTg@@a1n)*(calcξg@@a11n)*(calcHgmore@@a11n)

E20=-E2

E24=-E2

Term 3  (δP)

E3=(calcTo@@a1n)*(calcξo@@a11n)*(calcHo@@a11n)

E21=-E3

Term 4  (δT)

E4=TE

E22=-E4

Term 5 (δP)

E5=-VRdt*ϕm*Sg1ℓ*(cϕ*calcξg[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcDξgDP[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcξg[P1ℓ,T1ℓ]*(calcDUgmoreDP@@a11ℓ))

Term 6 (δP)

E6=-VRdt*ϕm*Sw1ℓ*(cϕ*calcξw[P1ℓ,T1ℓ]*(calcUwmore@@a11ℓ)+calcDξwDP[P1ℓ,T1ℓ]*(calcUwmore@@a11ℓ)+calcξw[P1ℓ,T1ℓ]*(calcDUwmoreDP@@a11ℓ))

Term 7 (δP)

E7=-VRdt*ϕm*So1ℓ*(cϕ*calcξo[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcDξoDP[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcξo[P1ℓ,T1ℓ]*(calcDUoDP@@a11ℓ))

Term 8 (δP)

E8=VRdt*ϕm*cϕ*ρR*CR*T1ℓ

Term 9 (δT)

E9=-VRdt*ϕm*Sg1ℓ*(calcDξgDT[P1ℓ,T1ℓ]*(calcUgmore@@a11ℓ)+calcξg[P1ℓ,T1ℓ]*(calcDUgmoreDT@@a11ℓ))

Term 10 (δT)

E10=-VRdt*ϕm*Sw1ℓ*(calcDξwDT[P1ℓ,T1ℓ]*(calcUwmore@@a11ℓ)+calcξw[P1ℓ,T1ℓ]*(calcDUwmoreDT@@a11ℓ))

Term 11 (δT)

E11=-VRdt*ϕm*So1ℓ*(calcDξoDT[P1ℓ,T1ℓ]*(calcUo@@a11ℓ)+calcξo[P1ℓ,T1ℓ]*(calcDUoDT@@a11ℓ))

Term 12 (δT)

E12=-VRdt*(1-ϕm)*ρR*CR

Term 13 (δSg)

E13=VRdt*ϕm*(calcξw[P1ℓ,T1ℓ]calcUwmore[P1ℓ,T1ℓ]-calcξg[P1ℓ,T1ℓ]calcUgmore[P1ℓ,T1ℓ])

Term 14 (δSo)

E14=VRdt*ϕm*(calcξw[P1ℓ,T1ℓ]calcUwmore[P1ℓ,T1ℓ]-calcξo[P1ℓ,T1ℓ]calcUo[P1ℓ,T1ℓ])

Energy Equation, RHS

Term 15 (ℓ)

E15 = VRdt*ϕm*(Sw1ℓ*calcξw[P1ℓ,T1ℓ]*calcUwmore[P1ℓ,T1ℓ]+
So1ℓ*calcξo[P1ℓ,T1ℓ]*calcUo[P1ℓ,T1ℓ]+
Sg1ℓ*calcξg[P1ℓ,T1ℓ]*calcUgmore[P1ℓ,T1ℓ])

Term 16 (n)

E16 = -VRdt*ϕm*(Sw1n*calcξw[P1n,T1n]*calcUwmore[P1n,T1n]+
So1n*calcξo[P1n,T1n]*calcUo[P1n,T1n]+
Sg1n*calcξg[P1n,T1n]*calcUgmore[P1n,T1n])

Term 17 (ℓ)

E17=VRdt*(1-ϕm)*ρR*CR*T1ℓ

Term 18 (n)

E18=-VRdt*(1-ϕm)*ρR*CR*T1n

Term 25,26,27 =0 since ΔD=0

E25=E26=E27=0

Term 28-29, cell 1 injection

E28a=-heatsource

Term 28-29, cell N production

E28c=-(-heatsource)*((calcffw@@a1n)*calcξw[P1n,T1n]*calcHwless[P1n,T1n]+
(calcffo@@a1n)*calcξo[P1n,T1n]*calcHo[P1n,T1n]+
(calcffg@@a1n)*calcξg[P1n,T1n]*calcHgmore[P1n,T1n])

Term 30, to boundaries

E30=-(calcqb[T1n]+calcqb[T1n])

Thermodynamic Constraints

Term 1; δP

G1=-1;

Term 2; δT

G2=calcDGDT[P1ℓ,T1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]

Term 3; δSo

G3=calcDGDSo[P1ℓ,T1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]

Term 4; δSg

G4=G3

Term 4; (ℓ)

G5=-(calcG[P1ℓ,T1ℓ,So1ℓ,Sg1ℓ,Sw1ℓ]+computePcowm[Sw1n])