In the oil and gas industry, maximizing the efficiency and productivity of oil wells is crucial. One of the most effective techniques for achieving this is NODAL analysis. This method allows engineers to evaluate and optimize the performance of oil wells by analyzing the inflow performance relationship (IPR) and vertical lift performance (VLP) curves.
In today’s tutorial, you’ll explore how to perform NODAL analysis using Python. We’ll start by calculating production rates with the Vogel IPR function and then plot the VLP curves. By combining these two analyses, we can determine the operating point where the well performs optimally. To make this process interactive and user-friendly, we’ll create a graphical interface using the tkinter library.
Whether you are a petroleum engineer looking to refine your reservoir management skills or a Python enthusiast eager to apply your programming knowledge to real-world problems, this article will provide you with valuable insights and practical tools for optimizing oil well production. Let’s dive into the fascinating world of NODAL analysis and see how Python can help us unlock the full potential of our oil wells.
Table of Contents
- Getting Started
- Imports
- Function to Calculate Vogel IPR Curve
- Function to Calculate VLP Example Curve
- Function-for-NODAL-Analysis
- Plotting IPR and VLP Curves
- Updating the Plot Based on User Inputs
- Main Application Setup
- Example
- Full Code
Getting Started
To get this code working, you’ll need to install a few libraries. Just open your terminal or command prompt and run the following commands:
$ pip install tk
$ pip install numpy
$ pip install pandas
$ pip install matplotlib
Imports
import tkinter as tk
from tkinter import ttk
import numpy as np
import pandas as pd
from matplotlib.figure import Figure
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
Well, as usual, before heading out on any adventure, you need to gear up with the right tools. Here’s what we’ll need:
- First, we have tkinter to create our graphical user interface.
- Next,
ttkhelps us use those sleek-themed widgets. - numpy is our go-to for handling arrays and all sorts of mathematical functions.
- pandas steps in for data analysis and manipulation, especially with tables.
matplotlib.figurelets us create stunning figures.- Finally,
FigureCanvasTkAggallows us to embed those Matplotlib figures right into our tkinter application.
Function to Calculate Vogel IPR Curve
# Function to calculate Vogel IPR curve
def vogel_ipr(q_max, p_wf, p_res):
return q_max * (1 - 0.2 * (p_wf / p_res) - 0.8 * (p_wf / p_res) ** 2)
After assembling our gear, the next step in our journey is to reveal the production capabilities of an oil well. This is where the vogel_ipr() function comes in. It calculates the production rate of a well based on the maximum production rate, the bottom-hole flowing pressure, and the reservoir pressure.
The formula for the Vogel IPR curve is:
Where:
- q: Production rate (STB/day)
- qmax: Maximum production rate (STB/day)
- pwf: Bottom-hole flowing pressure (psi)
- pres: Reservoir pressure (psi)
Function to Calculate VLP Example Curve
# Function to calculate VLP example curve
def vlp_example(depth, rate, fluid_density=50, gas_fraction=0.1):
gas_density = 0.0764 # lb/ft³, typical value for natural gas
liquid_density = fluid_density * (1 - gas_fraction) + gas_density * gas_fraction
pressure_drop = 0.433 * liquid_density * depth / 1440 * rate / 1000 # Pressure drop in psi
return pressure_drop
The next step is to understand the pressure drop that occurs as fluid flows up the well. This is the goal of the vlp_example() function, which calculates the pressure drop using a formula based on the well depth, production rate, fluid density, fraction of gas in the fluid, natural gas density, and the combined density of the fluid and gas.
The formula for the pressure drop is:
Where:
- Pdrop: Pressure drop (psi)
- depth is the depth of the well (in feet)
- rate is the production rate (in STB/day)
- fluid_density is the density of the fluid (in lb/ft³), with a default value of 50 lb/ft³
- gas_fraction is the fraction of gas in the fluid, with a default value of 0.1 (or 10%)
- gas_density is the density of natural gas, set as a constant 0.0764 lb/ft³
By applying this formula, the vlp_example() function helps us calculate the pressure drop along the well, which is crucial for determining the well’s operating point.
Function for NODAL Analysis
Now that we have calculated the formulas for our two curves, it’s time to find the operating point where these curves intersect. This is the goal of the nodal_analysis() function. It works by taking two sets of data: one for the Inflow Performance Relationship (IPR) and another for the Vertical Lift Performance (VLP).
The function first estimates the VLP pressures at the flow rates given in the IPR data. Then, it calculates the difference between the IPR pressures and these estimated VLP pressures. By finding the point where this difference is the smallest, we identify the operating point. This operating point is where the two curves intersect, representing the optimal flow rate and pressure for the well.
# Function for NODAL analysis to find the operating point
def nodal_analysis(ipr_df, vlp_df):
vlp_interp = np.interp(ipr_df['q'], vlp_df['rate'], vlp_df['pressure_drop'])
differences = np.abs(ipr_df['p_wf'] - vlp_interp)
min_index = differences.idxmin()
operating_point = (ipr_df.loc[min_index, 'q'], ipr_df.loc[min_index, 'p_wf'])
return operating_point
Plotting IPR and VLP Curves
With all the theoretical groundwork laid, we can now use the plot_ipr_vlp() function. This function takes the DataFrames for both IPR and VLP data, along with the operating point, and creates a figure. It plots both the IPR and VLP curves and marks their intersection, which is the operating point. The bottom-hole pressure is used for the x-axis, and the production rate is used for the y-axis.
# Function to plot IPR and VLP curves and the operating point
def plot_ipr_vlp(ipr_df, vlp_df, operating_point):
fig = Figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(ipr_df['p_wf'], ipr_df['q'], label='IPR Curve')
ax.plot(vlp_df['pressure_drop'], vlp_df['rate'], label='VLP Curve')
ax.plot(operating_point[1], operating_point[0], 'ro', label='Operating Point')
ax.set_title('NODAL Analysis')
ax.set_xlabel('Bottom Hole Pressure (psi)')
ax.set_ylabel('Production Rate (STB/day)')
ax.legend()
return fig
Updating the Plot Based on User Inputs
Since we want our programs to be as user-friendly as possible, this one is no exception. That’s why we created the update_plot() function. It takes the user’s inputs and uses the vogel_ipr() function to generate IPR data, storing it in a DataFrame. Then, it calls the vlp_example() function to generate VLP data, which it stores in another DataFrame.
Next, it uses the nodal_analysis() function to find the operating point and updates the label with the operating point details using op_label.config(). Finally, it clears the old plots from the figure and adds the new ones using FigureCanvasTkAgg(). This way, our program remains interactive and responsive to user input.
# Function to update the plot based on user inputs
def update_plot():
q_max = float(q_max_entry.get())
p_res = float(p_res_entry.get())
depth = float(depth_entry.get())
fluid_density = float(fluid_density_entry.get())
p_wf = np.linspace(0, p_res, 50)
q = vogel_ipr(q_max, p_wf, p_res)
ipr_df = pd.DataFrame({'p_wf': p_wf, 'q': q})
rates = np.linspace(0, q_max, 50)
pressure_drops = vlp_example(depth, rates, fluid_density)
vlp_df = pd.DataFrame({'rate': rates, 'pressure_drop': pressure_drops})
operating_point = nodal_analysis(ipr_df, vlp_df)
op_label.config(
text=f"Operating Point: Flow Rate = {operating_point[0]:.2f} STB/day, Bottom Hole Pressure = {operating_point[1]:.2f} psi"
)
# Clear the previous plot
for widget in plot_frame.winfo_children():
widget.destroy()
# Create a new plot
fig = plot_ipr_vlp(ipr_df, vlp_df, operating_point)
canvas = FigureCanvasTkAgg(fig, master=plot_frame)
canvas.draw()
canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
Main Application Setup
At last, we have reached the grand finale, where we put everything together in a graphical interface:
First, we create the main window, set its title, and define its geometry. Next, we set up a frame that will contain our inputs. Within that frame, we create input fields along with labels to clarify each field. We also add an “Update Plot” button that triggers the update_plot() function, and don’t forget the op_label that displays the operating point information whenever it’s updated.
Now that we are done with the inputs, it’s time to move on to the plots. We create a new frame dedicated to displaying the plots front and center. Finally, we start the main event loop to keep the main window running and responsive to the user with the mainloop() method.
# Main application
app = tk.Tk()
app.title("NODAL Analysis Optimizer - The Pycodes")
app.geometry("900x700")
input_frame = ttk.LabelFrame(app, text="Input Parameters")
input_frame.pack(side=tk.TOP, fill=tk.X, padx=10, pady=10)
# Input fields
ttk.Label(input_frame, text="Maximum Production Rate (STB/day):").grid(row=0, column=0, padx=5, pady=5, sticky=tk.W)
q_max_entry = ttk.Entry(input_frame)
q_max_entry.grid(row=0, column=1, padx=5, pady=5)
q_max_entry.insert(0, "2000")
ttk.Label(input_frame, text="Reservoir Pressure (psi):").grid(row=1, column=0, padx=5, pady=5, sticky=tk.W)
p_res_entry = ttk.Entry(input_frame)
p_res_entry.grid(row=1, column=1, padx=5, pady=5)
p_res_entry.insert(0, "3000")
ttk.Label(input_frame, text="Well Depth (feet):").grid(row=2, column=0, padx=5, pady=5, sticky=tk.W)
depth_entry = ttk.Entry(input_frame)
depth_entry.grid(row=2, column=1, padx=5, pady=5)
depth_entry.insert(0, "10000")
ttk.Label(input_frame, text="Fluid Density (lb/ft³):").grid(row=3, column=0, padx=5, pady=5, sticky=tk.W)
fluid_density_entry = ttk.Entry(input_frame)
fluid_density_entry.grid(row=3, column=1, padx=5, pady=5)
fluid_density_entry.insert(0, "50")
# Update button
update_button = ttk.Button(input_frame, text="Update Plot", command=update_plot)
update_button.grid(row=4, column=0, columnspan=2, pady=10)
# Operating point label
op_label = ttk.Label(app, text="Operating Point: Flow Rate = N/A, Bottom Hole Pressure = N/A")
op_label.pack(side=tk.TOP, fill=tk.X, padx=10, pady=10)
# Plot frame
plot_frame = ttk.Frame(app)
plot_frame.pack(side=tk.TOP, fill=tk.BOTH, expand=1)
# Run the application
app.mainloop()
Example
As you can see below, I ran this script on a Windows system:
Also on a Linux system:
Full Code
import tkinter as tk
from tkinter import ttk
import numpy as np
import pandas as pd
from matplotlib.figure import Figure
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
# Function to calculate Vogel IPR curve
def vogel_ipr(q_max, p_wf, p_res):
return q_max * (1 - 0.2 * (p_wf / p_res) - 0.8 * (p_wf / p_res) ** 2)
# Function to calculate VLP example curve
def vlp_example(depth, rate, fluid_density=50, gas_fraction=0.1):
gas_density = 0.0764 # lb/ft³, typical value for natural gas
liquid_density = fluid_density * (1 - gas_fraction) + gas_density * gas_fraction
pressure_drop = 0.433 * liquid_density * depth / 1440 * rate / 1000 # Pressure drop in psi
return pressure_drop
# Function for NODAL analysis to find the operating point
def nodal_analysis(ipr_df, vlp_df):
vlp_interp = np.interp(ipr_df['q'], vlp_df['rate'], vlp_df['pressure_drop'])
differences = np.abs(ipr_df['p_wf'] - vlp_interp)
min_index = differences.idxmin()
operating_point = (ipr_df.loc[min_index, 'q'], ipr_df.loc[min_index, 'p_wf'])
return operating_point
# Function to plot IPR and VLP curves and the operating point
def plot_ipr_vlp(ipr_df, vlp_df, operating_point):
fig = Figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(ipr_df['p_wf'], ipr_df['q'], label='IPR Curve')
ax.plot(vlp_df['pressure_drop'], vlp_df['rate'], label='VLP Curve')
ax.plot(operating_point[1], operating_point[0], 'ro', label='Operating Point')
ax.set_title('NODAL Analysis')
ax.set_xlabel('Bottom Hole Pressure (psi)')
ax.set_ylabel('Production Rate (STB/day)')
ax.legend()
return fig
# Function to update the plot based on user inputs
def update_plot():
q_max = float(q_max_entry.get())
p_res = float(p_res_entry.get())
depth = float(depth_entry.get())
fluid_density = float(fluid_density_entry.get())
p_wf = np.linspace(0, p_res, 50)
q = vogel_ipr(q_max, p_wf, p_res)
ipr_df = pd.DataFrame({'p_wf': p_wf, 'q': q})
rates = np.linspace(0, q_max, 50)
pressure_drops = vlp_example(depth, rates, fluid_density)
vlp_df = pd.DataFrame({'rate': rates, 'pressure_drop': pressure_drops})
operating_point = nodal_analysis(ipr_df, vlp_df)
op_label.config(
text=f"Operating Point: Flow Rate = {operating_point[0]:.2f} STB/day, Bottom Hole Pressure = {operating_point[1]:.2f} psi"
)
# Clear the previous plot
for widget in plot_frame.winfo_children():
widget.destroy()
# Create a new plot
fig = plot_ipr_vlp(ipr_df, vlp_df, operating_point)
canvas = FigureCanvasTkAgg(fig, master=plot_frame)
canvas.draw()
canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
# Main application
app = tk.Tk()
app.title("NODAL Analysis Optimizer - The Pycodes")
app.geometry("900x700")
input_frame = ttk.LabelFrame(app, text="Input Parameters")
input_frame.pack(side=tk.TOP, fill=tk.X, padx=10, pady=10)
# Input fields
ttk.Label(input_frame, text="Maximum Production Rate (STB/day):").grid(row=0, column=0, padx=5, pady=5, sticky=tk.W)
q_max_entry = ttk.Entry(input_frame)
q_max_entry.grid(row=0, column=1, padx=5, pady=5)
q_max_entry.insert(0, "2000")
ttk.Label(input_frame, text="Reservoir Pressure (psi):").grid(row=1, column=0, padx=5, pady=5, sticky=tk.W)
p_res_entry = ttk.Entry(input_frame)
p_res_entry.grid(row=1, column=1, padx=5, pady=5)
p_res_entry.insert(0, "3000")
ttk.Label(input_frame, text="Well Depth (feet):").grid(row=2, column=0, padx=5, pady=5, sticky=tk.W)
depth_entry = ttk.Entry(input_frame)
depth_entry.grid(row=2, column=1, padx=5, pady=5)
depth_entry.insert(0, "10000")
ttk.Label(input_frame, text="Fluid Density (lb/ft³):").grid(row=3, column=0, padx=5, pady=5, sticky=tk.W)
fluid_density_entry = ttk.Entry(input_frame)
fluid_density_entry.grid(row=3, column=1, padx=5, pady=5)
fluid_density_entry.insert(0, "50")
# Update button
update_button = ttk.Button(input_frame, text="Update Plot", command=update_plot)
update_button.grid(row=4, column=0, columnspan=2, pady=10)
# Operating point label
op_label = ttk.Label(app, text="Operating Point: Flow Rate = N/A, Bottom Hole Pressure = N/A")
op_label.pack(side=tk.TOP, fill=tk.X, padx=10, pady=10)
# Plot frame
plot_frame = ttk.Frame(app)
plot_frame.pack(side=tk.TOP, fill=tk.BOTH, expand=1)
# Run the application
app.mainloop()
Happy Coding!
