Part III: Reservoir & Production Engineering

Chapter 12

Nodal Analysis and Well Performance

schedule15 min readfitness_center10 exercises

Why This Chapter Exists

A well that has been drilled, cased, cemented, perforated, and completed is not automatically productive. The reservoir may contain millions of barrels of oil, but how much of it actually reaches the surface depends on the entire production system working together --- the reservoir rock, the perforations, the wellbore, the tubing, the choke, and the surface facilities. A restriction at any point in that chain limits the whole system.

Nodal analysis is the method petroleum engineers use to diagnose where the restriction is and what to do about it. It works by dividing the production system into components, modelling each one mathematically, and finding the operating point --- the flow rate and pressure at which the system naturally settles. If that operating rate is lower than expected, nodal analysis tells you why and points you toward the fix: a larger tubing size, a different choke setting, artificial lift, or a stimulation treatment.

This chapter builds the complete nodal analysis toolkit in Python. By the end, you will be able to construct inflow performance curves from reservoir data, model tubing performance under multiphase flow conditions, find the operating point for any well configuration, and run sensitivity analyses that quantify how each parameter affects production.

infoWhat You Will Learn

Construct IPR curves using the Vogel, Fetkovich, and composite models
Calculate tubing performance relationships for vertical multiphase flow
Find the operating point where inflow meets outflow
Evaluate artificial lift options quantitatively
Run sensitivity analyses and generate tornado charts
Diagnose underperforming wells systematically

The Production System

Before writing any code, it is important to understand what the production system actually consists of and why each component matters.

Fluid in the reservoir exists at high pressure --- typically thousands of psi --- held in place by the weight of the overlying rock and the surrounding formation. When a well is drilled and completed, it creates a pressure sink. Reservoir pressure pushes fluid through the porous rock toward the wellbore. This is inflow.

Once fluid enters the wellbore, it must travel up through the tubing to the surface. Along the way, it loses pressure to gravity (the weight of the fluid column), friction (the fluid rubbing against the pipe walls), and acceleration (changes in velocity, especially as gas comes out of solution and the fluid expands). This is outflow.

At the surface, fluid passes through a choke and into the processing facilities, where it is separated into oil, gas, and water. Back-pressure from the surface equipment further reduces the driving force available to move fluid from the reservoir to the surface.

The key insight of nodal analysis is this: the well produces at the rate where the pressure supplied by the reservoir exactly matches the pressure required to lift the fluid to the surface. If you plot inflow pressure against rate on one curve and outflow pressure against rate on another, the intersection is your operating point. That single point tells you your well's flow rate and flowing bottomhole pressure.

Inflow Performance Relationships

The inflow performance relationship (IPR) describes how much fluid the reservoir can deliver to the wellbore at a given flowing bottomhole pressure. It is the supply side of the system.

Productivity Index --- The Simplest Case

For an undersaturated oil reservoir (where the flowing bottomhole pressure stays above the bubble point), the relationship between rate and pressure is linear:

q_o = J \times (P_e - P_{wf})

where:

$q_o$ = oil flow rate (STB/day)
$J$ = productivity index (STB/day/psi)
$P_e$ = average reservoir pressure (psi)
$P_{wf}$ = flowing bottomhole pressure (psi)

The productivity index $J$ is a single number that captures everything about the reservoir's ability to deliver fluid: permeability, thickness, viscosity, wellbore radius, drainage radius, and skin factor. A higher $J$ means a more productive well. This linear model works well when the reservoir pressure is high and no gas comes out of solution in the formation.

main.py

import numpy as np
import matplotlib.pyplot as plt

# Reservoir parameters
Pe = 4500           # average reservoir pressure, psi
J = 8.5             # productivity index, STB/day/psi
Pb = 3200           # bubble point pressure, psi

# Generate IPR curve (linear model, valid above bubble point)
Pwf = np.linspace(0, Pe, 200)
qo = J * (Pe - Pwf)

# Key points
q_max = J * Pe                    # absolute open flow (Pwf = 0)
q_at_bubble = J * (Pe - Pb)      # rate at bubble point pressure

print(f"Reservoir pressure:      {Pe:,} psi")
print(f"Productivity index:      {J} STB/day/psi")
print(f"Bubble point pressure:   {Pb:,} psi")
print(f"Rate at bubble point:    {q_at_bubble:,.0f} STB/day")
print(f"Absolute open flow:      {q_max:,.0f} STB/day")

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(qo, Pwf, 'k-', linewidth=2)
ax.axhline(y=Pb, color='gray', linestyle='--', linewidth=0.8, label=f'Bubble Point ({Pb:,} psi)')
ax.plot(q_at_bubble, Pb, 'ro', markersize=8, label=f'q at Pb = {q_at_bubble:,.0f} STB/day')
ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Flowing Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Linear IPR — Undersaturated Reservoir', fontsize=12)
ax.legend(fontsize=9)
ax.set_xlim(0)
ax.set_ylim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The linear model is a starting point, but most wells do not operate entirely above the bubble point. Once the flowing bottomhole pressure drops below the bubble point, gas comes out of solution in the reservoir rock around the wellbore. This gas reduces the effective permeability to oil --- a phenomenon called the relative permeability effect --- and the IPR curve bends away from the straight line.

Vogel's IPR --- Below the Bubble Point

In 1968, Vogel published a dimensionless IPR equation for solution-gas drive reservoirs that accounts for the two-phase flow effect:

\frac{q_o}{q_{o,max}} = 1 - 0.2\left(\frac{P_{wf}}{P_e}\right) - 0.8\left(\frac{P_{wf}}{P_e}\right)^2

This equation is empirical --- it was derived from simulation studies, not from first principles --- but it has been validated extensively against field data and is one of the most widely used IPR models in the industry.

The parameter $q_{o,max}$ is the maximum theoretical flow rate if the bottomhole pressure were reduced to zero (absolute open flow potential). In practice, wells never operate at $P_{wf} = 0$ , but $q_{o,max}$ is a useful reference point for comparing wells.

main.py

import numpy as np
import matplotlib.pyplot as plt

def vogel_ipr(Pe, qo_max, Pwf_array):
    """
    Vogel's IPR for solution-gas drive reservoirs.

    Parameters
    ----------
    Pe : float
        Average reservoir pressure (psi).
    qo_max : float
        Absolute open flow potential (STB/day).
    Pwf_array : array-like
        Flowing bottomhole pressures (psi).

    Returns
    -------
    numpy.ndarray
        Oil flow rates (STB/day).
    """
    Pwf = np.asarray(Pwf_array)
    ratio = Pwf / Pe
    return qo_max * (1 - 0.2 * ratio - 0.8 * ratio**2)


def composite_ipr(Pe, Pb, J, Pwf_array):
    """
    Composite IPR: linear above bubble point, Vogel below.

    This is the standard approach when reservoir pressure
    exceeds the bubble point but flowing pressure may drop below it.
    """
    Pwf = np.asarray(Pwf_array)
    q = np.zeros_like(Pwf, dtype=float)

    # Rate at the bubble point (transition point)
    qb = J * (Pe - Pb)

    # Maximum rate using Vogel below bubble point
    qo_max = qb + J * Pb / 1.8

    for i, pw in enumerate(Pwf):
        if pw >= Pb:
            # Linear region (undersaturated)
            q[i] = J * (Pe - pw)
        else:
            # Vogel region (two-phase near wellbore)
            q[i] = qb + (qo_max - qb) * (1 - 0.2*(pw/Pb) - 0.8*(pw/Pb)**2)

    return q


# Well parameters
Pe = 4500       # reservoir pressure, psi
Pb = 3200       # bubble point, psi
J = 8.5         # productivity index, STB/day/psi

Pwf_range = np.linspace(0, Pe, 300)

# Calculate both models
q_linear = J * (Pe - Pwf_range)
q_composite = composite_ipr(Pe, Pb, J, Pwf_range)

# Absolute open flow for Vogel
qb = J * (Pe - Pb)
qo_max = qb + J * Pb / 1.8

print(f"Rate at bubble point:    {qb:,.0f} STB/day")
print(f"Absolute open flow:      {qo_max:,.0f} STB/day")
print(f"Difference from linear:  {J*Pe - qo_max:,.0f} STB/day less")
print()
print("The Vogel correction reduces the predicted flow rate")
print("at low bottomhole pressures because gas liberation in")
print("the reservoir impedes oil flow.")

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_linear, Pwf_range, 'b--', linewidth=1.5, alpha=0.5, label='Linear (no gas effect)')
ax.plot(q_composite, Pwf_range, 'k-', linewidth=2, label='Composite IPR')
ax.axhline(y=Pb, color='gray', linestyle=':', linewidth=0.8, label=f'Bubble Point ({Pb:,} psi)')
ax.plot(qb, Pb, 'ro', markersize=8, zorder=5)
ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Flowing Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Composite IPR — Linear Above Pb, Vogel Below', fontsize=12)
ax.legend(fontsize=9)
ax.set_xlim(0)
ax.set_ylim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The difference between the linear model and the composite model is significant. At low flowing pressures, the linear model overestimates the well's deliverability because it ignores the gas interference effect. Using the linear model alone to size equipment or plan production would lead to oversized facilities and disappointed production targets.

Fetkovich's IPR

Fetkovich proposed an alternative empirical model based on isochronal well test data:

q_o = C \times (P_e^2 - P_{wf}^2)^n

where $C$ is the deliverability coefficient and $n$ is the deliverability exponent (ranging from 0.5 to 1.0). When $n = 1$ , the equation reduces to a form similar to the linear PI. When $n = 0.5$ , it represents fully turbulent flow near the wellbore.

The Fetkovich model is particularly useful for gas wells and for oil wells where you have multi-rate test data to determine $C$ and $n$ .

main.py

import numpy as np
import matplotlib.pyplot as plt

def fetkovich_ipr(Pe, C, n, Pwf_array):
    """
    Fetkovich IPR model.

    Parameters
    ----------
    Pe : float
        Average reservoir pressure (psi).
    C : float
        Deliverability coefficient.
    n : float
        Deliverability exponent (0.5 to 1.0).
    Pwf_array : array-like
        Flowing bottomhole pressures (psi).

    Returns
    -------
    numpy.ndarray
        Flow rates (STB/day).
    """
    Pwf = np.asarray(Pwf_array)
    return C * (Pe**2 - Pwf**2)**n


Pe = 4500
Pwf_range = np.linspace(0, Pe, 300)

# Compare different n values
fig, ax = plt.subplots(figsize=(8, 5))
for n, C, style in [(1.0, 0.00042, '--'), (0.8, 0.0025, '-'), (0.5, 0.85, ':')]:
    q = fetkovich_ipr(Pe, C, n, Pwf_range)
    ax.plot(q, Pwf_range, style, linewidth=2, label=f'n = {n} (C = {C})')

ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Flowing Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Fetkovich IPR — Effect of Deliverability Exponent', fontsize=12)
ax.legend(fontsize=9)
ax.set_xlim(0)
ax.set_ylim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Tubing Performance Relationships

The outflow side of nodal analysis models what happens to the fluid as it travels from the bottom of the well to the surface through the production tubing. This is called the tubing performance relationship (TPR) or the vertical lift performance (VLP).

The pressure at the bottom of the tubing must overcome three components to deliver fluid to the surface:

P_{wf} = P_{wh} + \Delta P_{gravity} + \Delta P_{friction} - \Delta P_{acceleration}

where:

$P_{wh}$ = wellhead pressure (what the surface facilities impose)
$\Delta P_{gravity}$ = the hydrostatic weight of the fluid column in the tubing
$\Delta P_{friction}$ = frictional pressure losses
$\Delta P_{acceleration}$ = pressure changes due to velocity changes (usually small)

In single-phase liquid flow, this calculation is straightforward. In multiphase flow --- which is the norm in oil production, where oil, gas, and often water flow simultaneously --- the calculation becomes complex because the fluid density, velocity, and flow regime all change with depth as pressure decreases and gas comes out of solution.

Several correlations exist for multiphase vertical flow: Hagedorn and Brown, Beggs and Brill, Duns and Ros, and others. Each has strengths and limitations depending on the flow conditions. For this chapter, we implement a simplified gradient approach that captures the essential physics.

main.py

import numpy as np
import matplotlib.pyplot as plt

def tubing_performance(Pwh, depth_ft, tubing_id_in, oil_rate_stbd,
                       oil_sg=0.85, gas_sg=0.65, gor_scf_stb=500,
                       wc=0.1, oil_visc_cp=2.0):
    """
    Simplified tubing performance calculation.

    Estimates the required bottomhole flowing pressure
    for a given surface rate, accounting for gravity
    and friction in the tubing.

    This is a teaching implementation. Production environments
    use rigorous multiphase correlations (Beggs-Brill, etc.).

    Parameters
    ----------
    Pwh : float
        Wellhead pressure (psi).
    depth_ft : float
        True vertical depth (ft).
    tubing_id_in : float
        Tubing inner diameter (inches).
    oil_rate_stbd : float
        Oil production rate (STB/day).
    oil_sg : float
        Oil specific gravity (water = 1.0).
    gas_sg : float
        Gas specific gravity (air = 1.0).
    gor_scf_stb : float
        Gas-oil ratio (scf/STB).
    wc : float
        Water cut (fraction).
    oil_visc_cp : float
        Oil viscosity (cp).

    Returns
    -------
    float
        Required flowing bottomhole pressure (psi).
    """
    # Mixed fluid gradient (simplified)
    water_sg = 1.07
    liquid_sg = oil_sg * (1 - wc) + water_sg * wc

    # Gravity component
    # 0.433 psi/ft is the gradient for fresh water
    gravity_dp = 0.433 * liquid_sg * depth_ft

    # Friction component (Fanning equation, simplified)
    # Convert rate to velocity
    tubing_area_ft2 = np.pi * (tubing_id_in / 2 / 12)**2  # ft²
    liquid_rate_ft3_s = oil_rate_stbd * 5.615 / 86400       # ft³/s
    velocity_ft_s = liquid_rate_ft3_s / tubing_area_ft2

    # Friction factor (smooth pipe approximation)
    Re = 1488 * liquid_sg * 62.4 * velocity_ft_s * (tubing_id_in/12) / oil_visc_cp
    if Re < 2100:
        f = 64 / max(Re, 1)
    else:
        f = 0.316 / Re**0.25  # Blasius correlation

    friction_dp = f * depth_ft * liquid_sg * 62.4 * velocity_ft_s**2 / (2 * 32.174 * (tubing_id_in/12)) / 144

    Pwf = Pwh + gravity_dp + friction_dp
    return Pwf


# Generate TPR curve for a range of flow rates
Pwh = 200               # wellhead pressure, psi
depth = 9200             # TVD, ft
tubing_id = 2.441        # 2-7/8" tubing ID, inches

rates = np.linspace(100, 6000, 100)
Pwf_tpr = [tubing_performance(Pwh, depth, tubing_id, q) for q in rates]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(rates, Pwf_tpr, 'b-', linewidth=2, label=f'TPR ({tubing_id}" ID, Pwh={Pwh} psi)')
ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Required Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Tubing Performance Relationship', fontsize=12)
ax.legend(fontsize=9)
ax.set_xlim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The TPR curve shows that higher flow rates require higher bottomhole pressures --- the well has to push harder to move more fluid up the tubing. The curve rises steeply at high rates because friction losses increase with the square of velocity.

Finding the Operating Point

The operating point is where the IPR and TPR curves intersect. At this point, the pressure the reservoir can deliver at a given rate exactly matches the pressure the tubing system requires to lift that rate to the surface.

main.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq

# Reservoir parameters
Pe = 4500
Pb = 3200
J = 8.5

# Tubing parameters
Pwh = 200
depth = 9200
tubing_id = 2.441

def ipr_Pwf(q):
    """Return Pwf from IPR for a given rate (inverse of composite IPR)."""
    qb = J * (Pe - Pb)
    qo_max = qb + J * Pb / 1.8

    if q <= qb:
        return Pe - q / J
    else:
        # Solve Vogel equation for Pwf given q
        # q = qb + (qo_max - qb) * (1 - 0.2*(Pwf/Pb) - 0.8*(Pwf/Pb)**2)
        # Rearrange to find Pwf
        frac = (q - qb) / (qo_max - qb)
        # 0.8*x² + 0.2*x + (frac - 1) = 0 where x = Pwf/Pb
        a, b, c = 0.8, 0.2, (frac - 1)
        discriminant = b**2 - 4*a*c
        if discriminant < 0:
            return 0.0
        x = (-b + np.sqrt(discriminant)) / (2*a)
        return max(x * Pb, 0)


def tpr_Pwf(q):
    """Return Pwf from TPR for a given rate."""
    return tubing_performance(Pwh, depth, tubing_id, max(q, 1))


def residual(q):
    """Difference between IPR and TPR pressures at rate q."""
    return ipr_Pwf(q) - tpr_Pwf(q)


# Find operating point
qb = J * (Pe - Pb)
qo_max = qb + J * Pb / 1.8

# Search for intersection
q_op = brentq(residual, 100, qo_max * 0.99)
Pwf_op = ipr_Pwf(q_op)

print("=" * 45)
print("  NODAL ANALYSIS — OPERATING POINT")
print("=" * 45)
print(f"  Operating rate:     {q_op:,.0f} STB/day")
print(f"  Flowing BHP:        {Pwf_op:,.0f} psi")
print(f"  Drawdown:           {Pe - Pwf_op:,.0f} psi")
print(f"  Wellhead pressure:  {Pwh} psi")
print(f"  Tubing ID:          {tubing_id}\"")
print("=" * 45)

# Plot both curves with operating point
rates_plot = np.linspace(100, qo_max * 0.95, 300)
Pwf_ipr = [ipr_Pwf(q) for q in rates_plot]
Pwf_tpr = [tpr_Pwf(q) for q in rates_plot]

fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(rates_plot, Pwf_ipr, 'k-', linewidth=2, label='IPR (Reservoir Inflow)')
ax.plot(rates_plot, Pwf_tpr, 'b-', linewidth=2, label='TPR (Tubing Outflow)')
ax.plot(q_op, Pwf_op, 'ro', markersize=10, zorder=5, label=f'Operating Point ({q_op:,.0f} STB/d, {Pwf_op:,.0f} psi)')
ax.axhline(y=Pb, color='gray', linestyle=':', linewidth=0.8, alpha=0.5, label=f'Bubble Point ({Pb:,} psi)')
ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Flowing Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Nodal Analysis — IPR vs TPR', fontsize=13)
ax.legend(fontsize=9, loc='upper right')
ax.set_xlim(0)
ax.set_ylim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Sensitivity Analysis

The operating point depends on multiple parameters, and in practice, engineers need to know which parameters have the greatest impact on production rate. A sensitivity analysis systematically varies each parameter while holding the others constant, then measures the effect on the operating rate.

This is one of the most practical applications of nodal analysis. It answers questions like:

Would a larger tubing size increase production enough to justify the cost?
How much does wellhead pressure affect the rate?
What happens to production as reservoir pressure declines over time?
Would reducing skin (through stimulation) be more effective than installing artificial lift?

main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import brentq

def find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id):
    """Find the operating rate for given system parameters."""

    def ipr_pw(q):
        qb = J * (Pe - Pb)
        qo_max = qb + J * Pb / 1.8
        if q <= qb:
            return Pe - q / J
        else:
            frac = (q - qb) / (qo_max - qb)
            a, b, c = 0.8, 0.2, (frac - 1)
            disc = b**2 - 4*a*c
            if disc < 0:
                return 0.0
            x = (-b + np.sqrt(disc)) / (2*a)
            return max(x * Pb, 0)

    def tpr_pw(q):
        return tubing_performance(Pwh, depth, tubing_id, max(q, 1))

    def residual(q):
        return ipr_pw(q) - tpr_pw(q)

    qb = J * (Pe - Pb)
    qo_max = qb + J * Pb / 1.8

    try:
        return brentq(residual, 50, qo_max * 0.99)
    except ValueError:
        return 0.0


# Base case
base = {"Pe": 4500, "Pb": 3200, "J": 8.5, "Pwh": 200, "depth": 9200, "tubing_id": 2.441}
q_base = find_operating_rate(**base)

# Parameter variations
sensitivities = {
    "Reservoir Pressure (psi)":  ("Pe", [3500, 4000, 4500, 5000, 5500]),
    "Productivity Index (STB/d/psi)": ("J", [4.0, 6.0, 8.5, 12.0, 15.0]),
    "Wellhead Pressure (psi)":   ("Pwh", [100, 200, 300, 400, 500]),
    "Tubing ID (inches)":        ("tubing_id", [1.995, 2.441, 2.992, 3.476, 3.958]),
}

print(f"Base case operating rate: {q_base:,.0f} STB/day\n")
print(f"{'Parameter':<35} {'Value':>10} {'Rate (STB/d)':>14} {'Change':>10}")
print("-" * 75)

results = {}
for label, (param, values) in sensitivities.items():
    rates = []
    for v in values:
        params = base.copy()
        params[param] = v
        q = find_operating_rate(**params)
        rates.append(q)
        marker = "<-- base" if v == base[param] else ""
        print(f"  {label:<33} {v:>10.3f} {q:>14,.0f} {q - q_base:>+10,.0f}  {marker}")
    results[label] = (values, rates)
    print()

# Tornado chart — impact of each parameter
fig, ax = plt.subplots(figsize=(9, 5))
labels = []
low_vals = []
high_vals = []
for label, (values, rates) in results.items():
    labels.append(label)
    low_vals.append(min(rates) - q_base)
    high_vals.append(max(rates) - q_base)

y_pos = range(len(labels))
ax.barh(y_pos, high_vals, align='center', color='#2E86AB', alpha=0.8, label='Upside')
ax.barh(y_pos, low_vals, align='center', color='#D7263D', alpha=0.8, label='Downside')
ax.axvline(x=0, color='k', linewidth=0.8)
ax.set_yticks(y_pos)
ax.set_yticklabels(labels, fontsize=10)
ax.set_xlabel('Change in Operating Rate (STB/day)', fontsize=11)
ax.set_title('Tornado Chart — Sensitivity of Operating Rate', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.show()

The tornado chart shows at a glance which parameters most affect production. This is the tool that supports capital allocation decisions: if tubing size has the biggest impact, the investment goes toward a workover to upsize the tubing. If productivity index dominates, the money goes toward stimulation or re-perforation.

Artificial Lift Selection

When a well cannot flow naturally --- either because reservoir pressure has declined or because the fluid is too heavy for the available pressure to lift --- artificial lift is required. The main options are:

Rod pump (sucker rod) --- a surface-mounted beam pump that mechanically lifts fluid. Common for low-rate onshore wells.
Electrical submersible pump (ESP) --- a downhole centrifugal pump driven by an electric motor. Handles high rates but is expensive to replace.
Gas lift --- gas is injected into the tubing to reduce the fluid density, allowing reservoir pressure to push the lightened column to the surface.
Progressive cavity pump (PCP) --- a positive displacement pump suited to viscous oils and sandy conditions.

Each has a range of applicability defined by rate, depth, fluid properties, and economics. Nodal analysis helps quantify the benefit of each option by modifying the TPR curve and recalculating the operating point.

main.py

import numpy as np
import matplotlib.pyplot as plt

# Uses ipr_Pwf() and tubing_performance() from earlier code blocks
# Demonstrate gas lift effect on TPR
# Gas injection reduces the average fluid density in the tubing,
# lowering the gravity component of pressure drop.

def tpr_with_gas_lift(Pwh, depth, tubing_id, rate, gas_injection_rate_mscfd=0):
    """
    Simplified TPR with gas lift effect.

    Gas injection reduces effective fluid density in the tubing,
    lowering the required bottomhole pressure for a given rate.
    """
    base_Pwf = tubing_performance(Pwh, depth, tubing_id, rate)

    # Gas lightening effect (simplified)
    # Each Mscf/d of injected gas reduces the effective liquid gradient
    if gas_injection_rate_mscfd > 0:
        lightening_factor = 1 / (1 + 0.0003 * gas_injection_rate_mscfd)
        gravity_reduction = 0.433 * 0.85 * depth * (1 - lightening_factor)
        return base_Pwf - gravity_reduction
    return base_Pwf


Pwh = 200
depth = 9200
tubing_id = 2.441
rates = np.linspace(100, 6000, 200)

fig, ax = plt.subplots(figsize=(9, 6))

# Natural flow TPR
tpr_natural = [tubing_performance(Pwh, depth, tubing_id, q) for q in rates]
ax.plot(rates, tpr_natural, 'b-', linewidth=2, label='Natural Flow')

# Gas lift at different injection rates
for gl_rate, color, ls in [(500, '#2E86AB', '--'), (1000, '#27AE60', '--'), (2000, '#E67E22', '--')]:
    tpr_gl = [tpr_with_gas_lift(Pwh, depth, tubing_id, q, gl_rate) for q in rates]
    ax.plot(rates, tpr_gl, color=color, linestyle=ls, linewidth=1.5,
            label=f'Gas Lift @ {gl_rate} Mscf/d')

# IPR for reference
Pwf_ipr = [ipr_Pwf(q) for q in rates if q < qo_max * 0.95]
ax.plot(rates[:len(Pwf_ipr)], Pwf_ipr, 'k-', linewidth=2, label='IPR')

ax.set_xlabel('Oil Rate (STB/day)', fontsize=11)
ax.set_ylabel('Flowing Bottomhole Pressure (psi)', fontsize=11)
ax.set_title('Effect of Gas Lift on Tubing Performance', fontsize=12)
ax.legend(fontsize=9)
ax.set_xlim(0)
ax.set_ylim(0, 5000)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print("Gas lift shifts the TPR curve downward, reducing the")
print("required bottomhole pressure at every rate. The operating")
print("point moves to a higher rate and lower flowing pressure.")

Exercises

fitness_center

Exercise 12.1Practice

-- Linear IPR

A well has a productivity index of 12 STB/day/psi and a reservoir pressure of 3,800 psi. Calculate and plot the IPR curve. Determine the flow rate at ...

arrow_forward