Anthony Jimenez
21 August 2021
http://www.fa-jimenez.com/
import numpy as np
import pandas as pd
from scipy.optimize import shgo
from plotly.subplots import make_subplots
import plotly.graph_objects as go
Here we must recall that the Dranchuk and Abou-Kassem (DAK) correlation is an implicit equation that was created by fitting 1500 data points as an approximation of the Standing and Katz prior work.
Being that the equation is implicit, that means that we must solve the equation iteratively for any combination of $T_{pr}$ and $p_{pr}$. To do this, we will use the scipy.optimize function shgo(). This is an algorithm formally known as Simplical Homology Global Optimization. This global optimization function will seek to minimize the residual objective function that we define in the user-created function z_factor_solver(). An important key here is to make the objective function the abs() of the residual, because if not, the optimizer will find an erroneous answer by forcing the residual to a negative value.
# Inner z_factor_solver() function that computes coefficients
def z_factor_coefficients(z, T_pr, p_pr):
density = 0.27 * p_pr / (z * T_pr)
c1 = 0.3265 - (1.0700 / T_pr) - (0.5339 / T_pr**3) + (0.01569 / T_pr**4) - (0.05165 / T_pr**5)
c2 = 0.5475 - (0.7361 / T_pr) + (0.1844 / T_pr**2)
c3 = 0.1056 * ((-0.7361 / T_pr) + (0.1844 / T_pr**2))
c4 = 0.6134 * (1 + 0.7210 * density**2) * (density**2 / T_pr**3) * np.exp(-0.7210 * density**2)
out = [density, c1, c2, c3, c4]
return out
# Main z_factor_solver() that computes the residual
def z_factor_solver(z, T_pr, p_pr):
# Compute the z-factor coefficients from the inputs
z_coeff = z_factor_coefficients(z, T_pr, p_pr)
density, c1, c2, c3, c4 = z_coeff[0], z_coeff[1], z_coeff[2], z_coeff[3], z_coeff[4]
# Compute the residual of the z-factor equation that is implicit in z
res = abs(z - 1 - c1*density - c2*density**2 + c3*density**5 - c4)
return res
# Main function call that solves the z_factor_solver() function and returns the z-factor
def z_factor_func(T_pr, p_pr):
bounds = [(0.001, 3)]
return shgo(z_factor_solver, bounds=bounds, args=(T_pr, p_pr)).x[0]
The user-created function matches quite well with published Standing and Katz graphs. Other sets of pressure and temperature values can be tested and compared against published data.
# Example of the function call using user-defined pseudopressure and pseudotemperature values
p_pr = 4.87
T_pr = 1.8
z = z_factor_func(T_pr, p_pr)
z
# Main graph that generates all the z-factor traces
def generate_z_factor_graph():
# List values that will be used to generate graph
T_pr_graph = [1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3.0]
p_pr_graph = np.linspace(0, 15, 500)
# Main loop to populate the z-factors
for idx, T_pr in enumerate(T_pr_graph):
# Empty array to populate all p_pr for each T_pr
z = np.zeros(len(p_pr_graph))
T_pr_label = len(p_pr_graph) * [T_pr]
# Loop for z-factor calculations for each T_pr
for jdx, p_pr in enumerate(p_pr_graph):
z[jdx] = z_factor_func(T_pr, p_pr)
# Record z-factors to dataframe
if idx == 0:
zdf = pd.DataFrame(data=[T_pr_label, p_pr_graph, z]).T
zdf.columns = ['T_pr', 'p_pr', 'z']
else:
zdf_temp = pd.DataFrame(data=[T_pr_label, p_pr_graph, z]).T
zdf_temp.columns = ['T_pr', 'p_pr', 'z']
zdf = zdf.append(zdf_temp)
return zdf
z = generate_z_factor_graph()
# Illustration of z-factor dataframe
z.head()
# Generate subplots that show zoomed-in views and total views of the z-factor
fig = make_subplots(rows=3, cols=1)
# Sub-filter the dataframes to make different subplots
z_plot1 = z.loc[z['p_pr'] <= 8]
z_plot2 = z.loc[z['p_pr'] >= 7]
# Make subplot #1
for idx, T_pr in enumerate(z_plot1.T_pr.unique()):
fig.add_trace(
go.Scatter(x=z_plot1[z_plot1['T_pr'] == T_pr]['p_pr'], y=z_plot1[z_plot1['T_pr'] == T_pr]['z'],
name=f'T= {T_pr}', showlegend=False
), row=1, col=1
)
# Make subplot #2
for idx, T_pr in enumerate(z_plot2.T_pr.unique()):
fig.add_trace(
go.Scatter(x=z_plot2[z_plot2['T_pr'] == T_pr]['p_pr'], y=z_plot2[z_plot2['T_pr'] == T_pr]['z'],
name=f'T= {T_pr}', showlegend=False
), row=2, col=1
)
# Make subplot #3
for idx, T_pr in enumerate(z.T_pr.unique()):
fig.add_trace(
go.Scatter(x=z[z['T_pr'] == T_pr]['p_pr'], y=z[z['T_pr'] == T_pr]['z'],
name=f'T= {T_pr}', showlegend=False
), row=3, col=1
)
fig['layout']['xaxis']['title'] = r'$\text{Pseudoreduced Pressure, } p_{pr}$'
fig['layout']['yaxis']['title'] = r'$\text{Compressibility Factor, } z$'
fig['layout']['xaxis2']['title'] = r'$\text{Pseudoreduced Pressure, } p_{pr}$'
fig['layout']['yaxis2']['title'] = r'$\text{Compressibility Factor, } z$'
fig['layout']['xaxis3']['title'] = r'$\text{Pseudoreduced Pressure, } p_{pr}$'
fig['layout']['yaxis3']['title'] = r'$\text{Compressibility Factor, } z$'
fig.update_layout(height=1500, title_text='Z-Factor Graphs Using Dranchuk and Abou-Kassem Correlation')
fig.show()