import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline
# Productivity Index (darcy law)
def j_darcy(ko, h, bo, uo, re, rw, s, flow_regime = 'seudocontinuo'):
if flow_regime == 'seudocontinuo':
J_darcy = ko * h / (141.2 * bo * uo * (np.log(re / rw) - 0.75 + s))
elif flow_regime == 'continuo':
J_darcy = ko * h / (141.2 * bo * uo * (np.log(re / rw) + s))
return J_darcy
# Índice de Productividad
def j(q_test, pwf_test, pr, pb, ef=1, ef2=None):
if ef == 1:
if pwf_test >= pb:
J = q_test / (pr - pwf_test)
else:
J = q_test / ((pr - pb) + (pb / 1.8) * \
(1 - 0.2 * (pwf_test / pb) - 0.8 * (pwf_test / pb)**2))
elif ef != 1 and ef2 is None:
if pwf_test >= pb:
J = q_test / (pr - pwf_test)
else:
J = q_test / ((pr - pb) + (pb / 1.8) * \
(1.8 * (1 - pwf_test / pb) - 0.8 * ef * (1 - pwf_test / pb)**2))
elif ef !=1 and ef2 is not None:
if pwf_test >= pb:
J = (q_test / (pr - pwf_test) / ef) * ef2
else:
J = (q_test / ((pr - pb) + (pb / 1.8) * \
(1.8 * (1 - pwf_test / pb) - 0.8 * ef * \
(1 - pwf_test / pb)**2)) / ef) * ef2
return J
# Q(bpd) @ Pb
def Qb(q_test, pwf_test, pr, pb, ef=1, ef2=None):
qb = j(q_test, pwf_test, pr, pb, ef, ef2) * (pr - pb)
return qb
# AOF(bpd)
def aof(q_test, pwf_test, pr, pb, ef=1, ef2=None):
if (ef == 1 and ef2 is None):
if pr > pb: # Yac. subsaturado
if pwf_test >= pb:
AOF = j(q_test, pwf_test, pr, pb) * pr
elif pwf_test < pb:
AOF = Qb(q_test, pwf_test, pr, pb, ef=1) + ((j(q_test, pwf_test, pr, pb) * pb) / (1.8))
elif pr <= pb: # Yac. Saturado
AOF = q_test / (1 - 0.2 * (pwf_test / pr) - 0.8 * (pwf_test / pr)**2)
elif (ef < 1 and ef2 is None):
if pr > pb:
if pwf_test >= pb:
AOF = j(q_test, pwf_test, pr, pb, ef) * pr
elif pwf_test < pb:
AOF = Qb(q_test, pwf_test, pr, pb, ef) + ((j(q_test, pwf_test, pr, pb, ef) * pb) / (1.8)) * (1.8 - 0.8 * ef)
elif pr <= pb:
AOF = q_test / (1.8 * ef * (1 - pwf_test/pr) - 0.8 * ef**2 * (1 - pwf_test/pr)**2) * (1.8 * ef - 0.8 * ef**2)
elif (ef > 1 and ef2 is None):
if pr > pb:
if pwf_test >= pb:
AOF = j(q_test, pwf_test, pr, pb, ef) * pr
elif pwf_test < pb:
AOF = Qb(q_test, pwf_test, pr, pb, ef) + ((j(q_test, pwf_test, pr, pb, ef) * pb) / (1.8)) * (0.624 + 0.376 * ef)
elif pr <= pb:
AOF = q_test / (1.8 * ef * (1 - pwf_test/pr) - 0.8 * ef**2 * (1 - pwf_test/pr)**2) * (0.624 + 0.376 * ef)
elif (ef < 1 and ef2 >= 1):
if pr > pb:
if pwf_test >= pb:
AOF = j(q_test, pwf_test, pr, pb, ef, ef2) * pr
elif pwf_test < pb:
AOF = Qb(q_test, pwf_test, pr, pb, ef, ef2) + ((j(q_test, pwf_test, pr, pb, ef, ef2) * pb) / (1.8)) * (0.624 + 0.376 * ef2)
elif pr <= pb:
AOF = q_test / (1.8 * ef * (1 - pwf_test/pr) - 0.8 * ef**2 * (1 - pwf_test/pr)**2) * (0.624 + 0.376 * ef2)
elif (ef > 1 and ef2 <= 1):
if pr > pb:
if pwf_test >= pb:
AOF = j(q_test, pwf_test, pr, pb, ef, ef2) * pr
elif pwf_test < pb:
AOF = Qb(q_test, pwf_test, pr, pb, ef, ef2) + (j(q_test, pwf_test, pr, pb, ef, ef2) * pb / 1.8) * (1.8 - 0.8 * ef2)
elif pr <= pb:
AOF = q_test / (1.8 * ef * (1 - pwf_test/pr) - 0.8 * ef**2 * (1 - pwf_test/pr)**2) * (1.8 * ef - 0.8 * ef2**2)
return AOF
# Qo (bpd) @ Darcy Conditions
def qo_darcy(q_test, pwf_test, pr, pwf, pb):
qo = j(q_test, pwf_test, pr, pb) * (pr - pwf)
return qo
#Qo(bpd) @ vogel conditions
def qo_vogel(q_test, pwf_test, pr, pwf, pb):
if pr <= pb: # Yac. Saturado
qo = aof(q_test, pwf_test, pr, pb) * \
(1 - 0.2 * (pwf / pr) - 0.8 * ( pwf / pr)**2)
return qo
#Qo(bpd) @ vogel conditions
def qo_ipr_compuesto(q_test, pwf_test, pr, pwf, pb):
if pr > pb: # Yac. subsaturado
if pwf >= pb:
qo = qo_darcy(q_test, pwf_test, pr, pwf, pb)
elif pwf < pb:
qo = Qb(q_test, pwf_test, pr, pb) + \
((j(q_test, pwf_test, pr, pb) * pb) / (1.8)) * \
(1 - 0.2 * (pwf / pb) - 0.8 * ( pwf / pb)**2)
elif pr <= pb: # Yac. Saturado
qo = aof(q_test, pwf_test, pr, pb) * \
(1 - 0.2 * (pwf / pr) - 0.8 * ( pwf / pr)**2)
return qo
# Datos de un yacimiento saturado
pr = 2400 #psia
pb = 2500 #psia
pwf = 1000 #psia
q_test = 100 #stb/d
pwf_test = 1800 #psia
# Crear Dataframe
df = pd.DataFrame()
df['Pwf(psia)'] = np.array([2400, 2000, 1500, 1000, 500, 0])
df['Qo(bpd)'] = df['Pwf(psia)'].apply(lambda x: qo_vogel(q_test, pwf_test, pr, x, pb))
# Plot
fig, ax = plt.subplots(figsize=(20, 10))
x = df['Qo(bpd)']
y = df['Pwf(psia)']
X_Y_Spline = make_interp_spline(x, y) # This step is used to smooth the curve
X_ = np.linspace(x.min(), x.max(), 500)
Y_ = X_Y_Spline(X_)
ax.plot(X_, Y_, c='g')
ax.set_xlabel('Qo(bpd)')
ax.set_ylabel('Pwf(psia)')
ax.set_title('IPR')
ax.set(xlim=(0, df['Qo(bpd)'].max() + 10), ylim=(0, df['Pwf(psia)'][0] + 100))
ax.grid()
plt.show()
# IPR Curve
def IPR_curve(q_test, pwf_test, pr, pwf:list, pb):
# Creating Dataframe
df = pd.DataFrame()
df['Pwf(psia)'] = pwf
df['Qo(bpd)'] = df['Pwf(psia)'].apply(lambda x: qo_ipr_compuesto(q_test, pwf_test, pr, x, pb))
fig, ax = plt.subplots(figsize=(20, 10))
x = df['Qo(bpd)']
y = df['Pwf(psia)']
# The following steps are used to smooth the curve
X_Y_Spline = make_interp_spline(x, y)
X_ = np.linspace(x.min(), x.max(), 500)
Y_ = X_Y_Spline(X_)
#Build the curve
ax.plot(X_, Y_, c='g')
ax.set_xlabel('Qo(bpd)')
ax.set_ylabel('Pwf(psia)')
ax.set_title('IPR')
ax.set(xlim=(0, df['Qo(bpd)'].max() + 10), ylim=(0, df['Pwf(psia)'][0] + 100))
# Arrow and Annotations
plt.annotate(
'Bubble Point', xy=(Qb(q_test, pwf_test, pr, pb), pb),xytext=(Qb(q_test, pwf_test, pr, pb) + 100, pb + 100) ,
arrowprops=dict(arrowstyle='->',lw=1)
)
# Horizontal and Vertical lines at bubble point
plt.axhline(y=pb, color='r', linestyle='--')
plt.axvline(x=Qb(q_test, pwf_test, pr, pb), color='r', linestyle='--')
ax.grid()
plt.show()
# Datos
pr = 4000 #psi
pb = 3000 #psi
q_test = 600 #bpd
pwf_test = 2000 #bpd
pwf = np.array([4000, 3500, 3000, 2500, 1000, 0])
IPR_curve(q_test, pwf_test, pr, pwf, pb)
# IPR Curve
def IPR_curve_methods(q_test, pwf_test, pr, pwf:list, pb, method, ef=1, ef2=None):
# Creating Dataframe
fig, ax = plt.subplots(figsize=(20, 10))
df = pd.DataFrame()
df['Pwf(psia)'] = pwf
if method == 'Darcy':
df['Qo(bpd)'] = df['Pwf(psia)'].apply(lambda x: qo_darcy(q_test, pwf_test, pr, x, pb))
elif method == 'Vogel':
df['Qo(bpd)'] = df['Pwf(psia)'].apply(lambda x: qo_vogel(q_test, pwf_test, pr, x, pb))
elif method == 'IPR_compuesto':
df['Qo(bpd)'] = df['Pwf(psia)'].apply(lambda x: qo_ipr_compuesto(q_test, pwf_test, pr, x, pb))
# Stand the axis of the IPR plot
x = df['Qo(bpd)']
y = df['Pwf(psia)']
# The following steps are used to smooth the curve
X_Y_Spline = make_interp_spline(x, y)
X_ = np.linspace(x.min(), x.max(), 500)
Y_ = X_Y_Spline(X_)
#Build the curve
ax.plot(X_, Y_, c='g')
ax.set_xlabel('Qo(bpd)')
ax.set_ylabel('Pwf(psia)')
ax.set_title('IPR')
ax.set(xlim=(0, df['Qo(bpd)'].max() + 10), ylim=(0, df['Pwf(psia)'].max() + 100))
# Arrow and Annotations
plt.annotate(
'Bubble Point', xy=(Qb(q_test, pwf_test, pr, pb), pb),xytext=(Qb(q_test, pwf_test, pr, pb) + 100, pb + 100) ,
arrowprops=dict(arrowstyle='->',lw=1)
)
# Horizontal and Vertical lines at bubble point
plt.axhline(y=pb, color='r', linestyle='--')
plt.axvline(x=Qb(q_test, pwf_test, pr, pb), color='r', linestyle='--')
ax.grid()
plt.show()
IPR_curve_methods(q_test, pwf_test, pr, pwf, pb, 'Darcy')
IPR_curve_methods(q_test, pwf_test, pr, pwf, pb, 'IPR_compuesto')
# Datos de un yacimiento saturado
pr = 2400 #psia
pb = 2500 #psia
pwf = np.array([2400, 2000, 1500, 1000, 0]) #psia
q_test = 100 #stb/d
pwf_test = 1800 #psia
IPR_curve_methods(q_test, pwf_test, pr, pwf, pb, 'Vogel')
# Datos
pr = 120 #bar
pb = 65 #bar
q_test = 400 #m3/d
pwf_test = 100 #bar
# Datos
pr = 4000 #psi
pb = 3000 #psi
q_test = 600 #bpd
pwf_test = 2000 #bpd
pwf = np.array([4000, 3500, 3000, 2500, 1000, 0])
# Datos
pr = 3000 #psig
pb = 2130 #psig
q_test = 631 #stb/d
pwf_test = 1700 #psig
pwf = np.array([3000, 2500, 2000, 1000, 0])