Problem Dec 1 and Dec 3

Class_problems/Problem_13new/IMPES.py

@@ -0,0 +1,142 @@
+"""
+reservoir simulation assignment 10
+1D reservoir simulation Q10: Main file (Implicit pressure explicit saturation)
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 12/4/2020
+"""
+
+#import inbuilt libraries
+import numpy as np
+from scipy.sparse import lil_matrix, csr_matrix, identity
+from scipy.sparse.linalg import inv
+from scipy.sparse.linalg import spsolve
+import matplotlib.pyplot as plt
+import time as timer
+from math import floor, ceil
+import warnings
+warnings.filterwarnings("ignore")
+
+#importing personal libraries
+from input_file_2D import inputfile
+from myarrays import myarrays
+from updatewells import updatewells
+from rel_perm import rel_perm
+
+#making simulation and IC classes
+class numerical:
+ def __init__(self):
+ self.Bw = []
+class reservoir:
+ def __init__(self):
+ self.dt = [] 
+class fluid:
+ def __init__(self):
+ self.dt = [] 
+class grid:
+ def __init__(self):
+ self.xmin = []
+class BC:
+ def __init__(self):
+ self.xmin = []
+class IC:
+ def __init__(self):
+ self.xmin = []
+class petro:
+ def __init__(self):
+ self.xmin = []
+class well:
+ def __init__(self):
+ self.xmin = []
+
+#loading inputfile
+inputfile(fluid,reservoir,petro,numerical,BC,IC,well)
+#Implicit pressure and explicit saturation for update
+#looping through time
+t = np.empty((100000)) #dimensional time
+t[0]= 0
+t_D = np.empty((100000))#non dimensional time
+t_D[0]= 0
+k = 0
+PV = 0
+P = np.copy(IC.P)
+Pw = np.copy(IC.Pw)
+Sw = np.array(np.copy(IC.Sw))
+Sw_hyst=np.empty((numerical.N,2))
+nmax = ceil(numerical.tfinal / numerical.dt)
+
+fw = np.empty((100000)) #fractional flow of wetting phase
+fw[0]= 0
+P_plot= np.zeros((numerical.N,10000)) #matrix to save pressure 
+P_plot[:,0] = IC.P[:,0] 
+Sw_plot= np.zeros((numerical.N, 10000)) #matrix to save pressure 
+Sw_plot[:,0]= IC.Sw[:,0] 
+
+while k <0: #(t[k] < numerical.tfinal): #non dimensional time marching 
+ P_old = np.copy(P) #Placeholdering the old array
+ Sw_old= np.copy(Sw) #Placeholdering the old array
+
+ #Calculating the arrays
+ Tw, To, T, d11, d12, d21, d22, D, G, Pc = myarrays(fluid,reservoir,petro,numerical,BC,P,Pw,Sw,Sw_hyst)
+
+ #updating the wells
+ well, Qw, Qo, Jw, Jo = updatewells(reservoir,fluid,numerical,petro,P,Sw,well)
+
+ J = (-d22 @ inv(d12)) @ Jw + Jo
+ Q = (-d22 @ inv(d12)) @ Qw + Qo + 800.0 * J @ np.ones((numerical.N,1)) #Pwf = 800 psi
+
+ if numerical.method == 'IMPES':
+ IM = T + J + D #implicit part coefficient in Eq. 3.44 
+ EX = D @ P_old + Q + G #explicit part or RHS of Eq. 3.44
+
+ P = np.transpose([spsolve(IM,EX)]) #solving IM*P = EX or Ax=B 
+ Sw = Sw + inv(d12) @ (-Tw @ (P - (fluid.rhow/144.0) * numerical.D - Pc) - d11 @ (P - P_old) + Qw + Jw @ (800.0 - P)) #explicit saturation
+
+ for i in range(0, numerical.N):
+ if Sw[i,0] > Sw_old[i,0] and Sw_hyst[i,1] == 0: # [i,1] is a flag
+ Sw_hyst[i,0] = Sw[i,0]
+ Sw_hyst[i,1] = 1.0
+
+ elif Sw[i,0] < Sw_old[i,0]:
+ Sw_hyst[i,0] = Sw[i,0] 
+
+ k = k+1 
+ P_plot[:,k] = P[:,0] 
+ Sw_plot[:,k]= np.array(Sw)[:,0]
+ t[k]= t[k-1] + numerical.dt
+ t_D[k]= well.constraint[0][0]*t[k-1]/(reservoir.L*reservoir.W*reservoir.h*reservoir.phi[0,0])
+
+ kblock = numerical.N-1
+ krw,kro = rel_perm(petro,Sw[kblock,0])
+ M = (kro*fluid.muw[kblock,0])/(krw*fluid.muo[kblock,0])
+ fw[k] = 1/(1+M)
+
+P_plot[np.argwhere(reservoir.permx < 0.01)] = np.nan
+
+#post process
+'''
+plt.figure()
+plt.plot(numerical.xD,C_plot[:,21],'r-',label=f'Time=%0.2f' % t_D[21])
+plt.plot(numerical.xD,C_plot[:,102],'b-',label=f'Time=%0.2f' % t_D[102])
+plt.plot(numerical.xD,C_plot[:,203],'k-',label=f'Time=%0.2f' % t_D[203])
+plt.legend(loc='best', shadow=False, fontsize='medium')
+plt.xlabel(r'$x_D$')
+plt.ylabel(r'$C_D$')
+'''
+
+if numerical.Ny==1:
+ plt.figure()
+ plt.plot(numerical.xc/reservoir.L,Sw_plot[:,47],label=r'$t_D=0.1$')
+ plt.plot(numerical.xc/reservoir.L,Sw_plot[:,47*2],label=r'$t_D=0.2$')
+ plt.plot(numerical.xc/reservoir.L,Sw_plot[:,47*3],label=r'$t_D=0.3$')
+ plt.xlabel(r'$x_D$')
+ plt.ylabel(r'$S_w$')
+ plt.legend(loc='best', shadow=False, fontsize='medium')
+ plt.savefig('SwvsT.png',bbox_inches='tight', dpi = 600)
+
+
+ plt.figure()
+ plt.plot(t_D[0:k],fw[0:k])
+ plt.ylabel(r'Water cut')
+ plt.xlabel(r'Pore volumes injected')
+ plt.savefig('watercutvsPVI.png',bbox_inches='tight', dpi = 600)

Class_problems/Problem_13new/Thalf.py

@@ -0,0 +1,47 @@
+"""
+reservoir simulation project 2 (2020): Problem 1
+2D reservoir simulation: Interblock transmissibility
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 10/31/2020
+"""
+
+from rel_perm import rel_perm
+
+#fluid, reservoir and simulation parameters 
+def Thalf(i,j,direction,fluid,reservoir,petro,numerical,P,Sw):
+
+ if direction == 'x': 
+ kAd = (2*reservoir.permx[i,0]*numerical.dy[i,0]*reservoir.h*reservoir.permx[j,0]*numerical.dy[j,0]*reservoir.h)/ \
+ (reservoir.permx[i,0]*numerical.dy[i,0]*reservoir.h*numerical.dx[j,0] + \
+ reservoir.permx[j,0]*numerical.dy[j,0]*reservoir.h*numerical.dx[i,0])
+ elif direction == 'y': 
+ kAd = (2*reservoir.permy[i,0]*numerical.dx[i,0]*reservoir.h*reservoir.permy[j,0]*numerical.dx[j,0]*reservoir.h)/ \
+ (reservoir.permy[i,0]*numerical.dx[i,0]*reservoir.h*numerical.dy[j,0] + \
+ reservoir.permy[j,0]*numerical.dx[j,0]*reservoir.h*numerical.dy[i,0])
+
+ #Calculating the potential
+ POT_i = P[i,0] - ( fluid.rhoo[i,0] / 144.0 ) * numerical.D[i,0]
+ POT_j = P[j,0] - ( fluid.rhoo[j,0] / 144.0 ) * numerical.D[j,0]
+
+ if POT_i >= POT_j:
+ krw,kro = rel_perm(petro,Sw[i,0])
+ Bw = fluid.Bw[i,0]
+ Bo = fluid.Bo[i,0] 
+ muw = fluid.muw[i,0] 
+ muo = fluid.muo[i,0] 
+
+ elif POT_i < POT_j:
+ krw,kro = rel_perm(petro,Sw[j,0])
+ Bw = fluid.Bw[j,0]
+ Bo = fluid.Bo[j,0] 
+ muw = fluid.muw[j,0] 
+ muo = fluid.muo[j,0] 
+
+ fluidhalfw = krw/(muw * Bw)
+ fluidhalfo = kro/(muo * Bo) 
+
+ Tw = fluidhalfw*kAd
+ To = fluidhalfo*kAd
+
+ return Tw, To;

Class_problems/Problem_13new/assignment1_reservior_init.py

@@ -0,0 +1,140 @@
+"""
+reservoir simulation assignment 1
+Initializing a reservoir
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 9/3/2020
+"""
+
+#Importing required libraries
+
+import numpy as np #import numpy library
+import matplotlib.pyplot as plt #library for plotting
+import build_gridfun2D #personal library for making grid
+plt.rcParams.update({'font.size': 22})
+from matplotlib import cm
+
+#Making grid class
+class grid:
+ def __init__(self):
+ self.xmin = []
+ self.xmax = []
+ self.Nx = []
+ self.N = []
+
+#Input parameters
+length = 6000 # length of the reservoir [feet] 
+breadth = 7500 # breadth of the reservoir [feet] 
+Nx = 80 # cells in x-direction
+Ny = 75 # cells in y-direction
+thic = 50 # thickness of the reservoir [feet] 
+D_woc = 7475.45 # depth of water oil contact line [feet] 
+P_w_woc = 4500 # pressure at water oil contact line [psi] 
+rho_w_sc = 62.4 # density of water at standard conditions [lbm/ft^3] 
+rho_o_sc = 53.0 # density of oil at standard conditions [lbm/ft^3]
+rho_g_sc = 0.0458171 # density of gas at standard conditions [lbm/ft^3]
+c_w = 2.87e-6 # compressibility of water [1/psi] 
+c_o = 3e-5 # compressibility of oil [1/psi] 
+B_w = 1.012298811 # formation volume factor of water 
+B_o = 1.04567 # formation volume factor of oil 
+Rs = 90.7388 # solution gas oil ratio [ft^3/bbl]
+Pe = 3.5 # capillary entry pressure at water oil contact line [psi] 
+s_wr = 0.2 # residual water saturation
+s_or = 0.4 # residual oil saturation
+lam = 2 # model parameter for capillary pressure
+
+#Calculate water and oleic phase densities
+rho_w = rho_w_sc / B_w
+rho_o = (rho_o_sc + rho_g_sc * Rs /5.61)/B_o #5.61 ft^3 in a barrel (bbl)
+
+#Building grid
+grid.xmin = 0.0
+grid.xmax = length
+grid.Nx = Nx
+
+grid.ymin = 0.0
+grid.ymax = breadth
+grid.Ny = Ny
+
+build_gridfun2D.build_grid(grid) #building grid
+[Xc,Yc] = np.meshgrid(grid.xc,grid.yc) #building the (x,y) matrix
+
+#Importing depth text file
+depth = np.loadtxt("Inc_Depth.dat")
+depth[depth==0.0] = np.nan #removing the zero depth for nice plotting
+
+# Variation of pressure and saturation with depth (ignoring compressibility)
+P_w = P_w_woc + rho_w / 144.0 * (depth - D_woc) #144 in^2 in 1 ft^2 
+P_o = (P_w_woc + Pe) + rho_o / 144.0 * (depth - D_woc) 
+s_w = s_wr + (1-s_wr)*((P_o-P_w)/Pe) **(-lam) #From Corey-Brooks draining curve
+s_w[depth>=D_woc] = 1.0
+s_o = 1-s_w
+
+#converting to a column vector
+depth_col= np.reshape(np.transpose(depth), (grid.N,-1)) #building the single depth vector
+P_o_col = np.reshape(np.transpose(P_o), (grid.N,-1)) #building the single P_o vector
+P_w_col = np.reshape(np.transpose(P_w), (grid.N,-1)) #building the single P_w vector
+s_w_col = np.reshape(np.transpose(s_w), (grid.N,-1)) #building the single s_w vector
+s_o_col = np.reshape(np.transpose(s_o), (grid.N,-1)) #building the single s_o vector
+
+#Plotting starts here
+
+#(a) Plotting the pressure variation with depth
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.plot(P_w_col,depth_col,'b.',label=r'$P_{water}$ [psi]')
+plot = plt.plot(P_o_col,depth_col,'r.',label=r'$P_{oil}$ [psi]')
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.xlabel(r'$Pressure$ [psi]')
+plt.ylabel(r'$Depth, D$ [m]')
+plt.gca().invert_yaxis()
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'Pressures.png',bbox_inches='tight', dpi = 600)
+
+#(b) Plotting the water saturation variation with depth
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.plot(s_w_col,depth_col,'b.')
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.xlabel(r'$s_{w}$')
+plt.ylabel(r'$Depth, D$ [m]')
+plt.gca().invert_yaxis()
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'Saturation.png',bbox_inches='tight', dpi = 600)
+
+#(c) Plotting the 2D contour of depth
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+ax1= plt.contourf(Xc, Yc,depth,cmap=cm.coolwarm, antialiased=True)
+plt.title('Depth of the reservoir')
+plt.axis('scaled')
+plt.xlabel(r'$x [feet]$')
+plt.ylabel(r'$y [feet] $')
+clb = plt.colorbar(ax1)
+clb.set_label(r'$Depth,D [feet]$', labelpad=-40, y=1.1, rotation=0)
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'reservoir_depth.png',bbox_inches='tight', dpi = 600)
+
+#(d) Plotting the 2D contour of water saturation
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+ax1= plt.contourf(Xc, Yc,s_w,cmap=cm.coolwarm, antialiased=True)
+plt.title('Water saturation in the reservoir')
+plt.axis('scaled')
+plt.xlabel(r'$ x [feet]$')
+plt.ylabel(r'$y [feet] $')
+clb = plt.colorbar(ax1)
+clb.set_label(r'$s_w$', labelpad=-40, y=1.1, rotation=0)
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'reservoir_water_saturation.png',bbox_inches='tight', dpi = 600)
+
+#(e) Plotting the 2D contour of oil pressure
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+ax1= plt.contourf(Xc, Yc,P_o,cmap=cm.coolwarm, antialiased=True)
+plt.title('Oil pressure in the reservoir')
+plt.axis('scaled')
+plt.xlabel(r'$x [feet]$')
+plt.ylabel(r'$y [feet] $')
+clb = plt.colorbar(ax1)
+clb.set_label(r'$P_{oil} [psi]$', labelpad=-40, y=1.1, rotation=0)
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'reservoir_oil_pressure.png',bbox_inches='tight', dpi = 600)

Class_problems/Problem_13new/cap_press.py

@@ -0,0 +1,38 @@
+"""
+reservoir simulation project 2 (2020): Problem 1
+2D Multiphase reservoir simulation: 
+Capillary pressure and relative permeability: Capillary pressure
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 12/04/2020
+"""
+import numpy as np
+
+def cap_press(petro,Sw,Sw_hyst):
+
+ S = (Sw - petro.Swr)/(1.0 - petro.Sor - petro.Swr) #Normalized saturation {Eq. 1.29}
+ Se = (Sw - petro.Swr)/(1.0 - petro.Swr) #{Eq. 1.28a}
+
+ #Corey-Brooks model
+ Pcd = petro.Pe * Se **(-1.0 / petro.lamda) #Drainage capillary pressure (used for initialization) {Eq. 1.28a} 
+ Pci = petro.Pe * (S **(-1.0 / petro.lamda) - 1.0) #Imbibition capillary pressure (used for initialization) {Eq. 1.28b}
+
+ #Capillary pressure scanning curve
+ epspc = 1e-5
+ Sw_max = 1.0 - petro.Sor
+ #f = ((Sw_max - Sw_hyst + epspc)/(Sw_max - Sw_hyst)) * ((Sw - Sw_hyst) / (Sw - Sw_hyst + epspc))
+ f = 1.0
+ Pc = f * Pci + (1.0 - f) * Pcd
+
+ #Calculate derivative
+ S2 = (Sw + 0.001 - petro.Swr) / (1.0 - petro.Swr - petro.Sor)
+ Se2 = (Sw + 0.001 - petro.Swr) / (1.0 - petro.Swr)
+
+ Pcd2 = petro.Pe * Se2 **(-1.0 / petro.lamda)
+ Pci2 = petro.Pe * (S2 **(-1.0 / petro.lamda) - 1.0 )
+ #f2 = ((Sw_max - Sw_hyst + epspc) / (Sw_max - Sw_hyst)) * ((Sw + 0.001 - Sw_hyst) / (Sw + 0.001 - Sw_hyst + epspc))
+ f2 = 1.0
+ Pc2 = f2 * Pci + (1.0 - f2) * Pcd
+ Pcprime = (Pc2 - Pc)/0.001
+
+ return Pc,Pcprime;

Class_problems/Problem_13new/depth3x3.txt

@@ -0,0 +1 @@
+1883.2 2275.5 2300.6 1883.2 2275.5 2300.6 1883.2 2275.5 2300.6