Nodal Analysis

Definition: Nodal Analysis is a method used to evaluate a complete producing system by isolating a single point (the Node) and ensuring the pressure and flow rate are consistent across that point. In petroleum production, the most common node is the Bottom-hole, where the Inflow (IPR) meets the Outflow (VLP).

The Node Concept

For any given node, two conditions must be met:

Flow into the node equals flow out of the node.
Only one pressure can exist at the node at a given flow rate.

The Node Equations:

Inflow (Supply): $p_{node} = p_r - \Delta p_{reservoir}$
Outflow (Demand): $p_{node} = p_{surf} + \Delta p_{tubing} + \Delta p_{choke}$

Determining the Operating Point

The intersection of the IPR curve and the VLP curve represents the Operating Point. This is the only rate (:math:q_{actual}) at which the well will naturally flow for a given set of conditions.

Numerical Example:

Consider a well with:

Reservoir Pressure ($p_r$) = 3500 psi
Productivity Index ($J$) = 1.2 STB/day/psi
Surface Pressure ($p_{surf}$) = 250 psi
VLP is simplified as: $p_{wf} = p_{surf} + 0.00002 q^{1.8}; + \rho/144$

To find the operating point, we solve for $q$ where $p_{wf, IPR} = p_{wf, VLP}$.

//IPR Input
double q_max = 2000; // STB/day
double p_r = 2500; // psi
double p_wf = 1000; // psi

//IPR
double qfun(double pwf) => q_max * (1 - 0.2 * (pwf / p_r) - 0.8 * Pow(pwf / p_r, 2));
ColVec P_ipr = Linspace(0, p_r);
ColVec Q_ipr = Arrayfun(qfun, P_ipr);

// VLP Inputs
double p_surf = 200; // psi (Wellhead Pressure)
double depth = 2000; // ft


// VLP
double pressureGradient(double z, double p, double q)
{
    double density = 55.0; // lb/ft3 (Oil)
    double friction_grad = 0.00002 * Pow(q, 1.8); // Simplified friction term
    double hydro_grad = density / 144.0; // psi/ft
    return hydro_grad + friction_grad;
}
double pfun(double q_g)
{
    var (Z, P) = Ode45((z, p) => pressureGradient(z, p, q_g), p_surf, [0, depth]);
    double p_wf = P[^1]; // extract the pressure at the bottom
    return p_wf;
}

ColVec Q_vlp = Linspace(0, 800);
ColVec P_vlp = Arrayfun(pfun, Q_vlp);

(double operating_q, double operating_p) = Intersection(Q_vlp, P_vlp, Q_ipr, P_ipr);
// Assuming a simple bisection or search between 0 and AOF
Scatter(operating_q, operating_p, "for", 15); HoldOn();
Plot(Q_ipr, P_ipr, "b");
Plot(Q_vlp, P_vlp, "g"); HoldOff();
Legend(["Operating Condition", "IPR", "VLP"]);
SaveAs("Nodal_Analysis.png");

Sensitivity Analysis

Nodal analysis is most powerful when performing “What-If” scenarios. By shifting the curves, engineers can predict the impact of changes:

Stimulation Economics: A Nodal Analysis Case Study

Problem Statement: An oil company is evaluating two competing stimulation proposals to improve the productivity of a damaged well. The well is currently performing poorly due to a high skin factor ($s = +5$). The goal is to determine which intervention provides the best return on investment by calculating the production gain per million USD spent.

Reservoir and Well Data

Parameter	Value	Unit
Reservoir Pressure($p_r$)	2800	psia
Bubble Point($p_b$)	3000	psia
Current Skin Factor($s_{ current}$)	+5	dimensionless
Max Flow Rate(Ideal $q_{ max}$)	2000	STB/day
Reservoir Radius / Wellbore Radius($r_e/r_w$)	1000	dimensionless
Well Depth	4000	ft
Oil Density	55.0	$lb/ft^3$
Surface Pressure($p_{ surf}$)	200	psia

Stimulation Offers:

*Company A(Hydraulic Frac):**Reduces skin to :math:`-3`. Cost: **$10M. * Company B(Acid Wash):**Reduces skin to :math:`+1`. Cost: **$5M.

Solution - 1. Compute IPR and Flow Efficiency Since the reservoir pressure($p_r = 2800$) is below the bubble point($p_b = 3000$), the well follows Vogel’s non-linear behavior. We first determine the Flow Efficiency ($FE$) for each skin scenario. The relationship between skin and productivity adjustment is: $J_{ ratio} = \frac{\ln(r_e/r_w)}{\ln(r_e/r_w) + s}$

//IPR Input
double q_max = 2000; // STB/day
double p_r = 2800; // psi

//IPR
double qfun(double pwf) => q_max * (1 - 0.2 * (pwf / p_r) - 0.8 * Pow(pwf / p_r, 2));
ColVec P_ipr = Linspace(0, p_r);
ColVec Q_ipr = Arrayfun(qfun, P_ipr);

// Adjusting for Damage and Stimulation
double r_e_r_w = 1000;
double FE_Damaged = Log(r_e_r_w)/(Log(r_e_r_w) + 5);
double FE_Stimulation_A = Log(r_e_r_w)/(Log(r_e_r_w) + -3);
double FE_Stimulation_B = Log(r_e_r_w)/(Log(r_e_r_w) + 1);
ColVec Q_ipr_Damaged = Q_ipr * FE_Damaged;
ColVec Q_ipr_Stimulated_A = Q_ipr * FE_Stimulation_A;
ColVec Q_ipr_Stimulated_B = Q_ipr * FE_Stimulation_B;

Plot(Q_ipr_Damaged, P_ipr, "r", 2); HoldOn();
Plot(Q_ipr_Stimulated_A, P_ipr, "b", 2);
Plot(Q_ipr_Stimulated_B, P_ipr, "g", 2); HoldOff();
Legend(["Damaged", "Stimulated_A", "Stimulated_B"]);
SaveAs("IPR_Stimulation_Economics_Case_Study.png");

1. Compute the VLP

The Vertical Lift Performance is calculated by integrating the hydrostatic and frictional pressure drops from the surface to the bottom-hole.

// VLP Inputs
double p_surf = 200; // psi (Wellhead Pressure)
double depth = 4000; // ft


// VLP
double pressureGradient(double z, double p, double q)
{
    double density = 52.0; // lb/ft3 (Oil)
    double friction_grad = 1e-5 * Pow(q, 1.8); // Simplified friction term
    double hydro_grad = density / 144.0; // psi/ft
    return hydro_grad + friction_grad;
}
double pfun(double q_o)
{
    var (Z, P) = Ode45((z, p) => pressureGradient(z, p, q_o), p_surf, [0, depth]);
    double p_wf = P[^1]; // extract the pressure at the bottom
    return p_wf;
}

ColVec Q_vlp = Linspace(0, 1600);
ColVec P_vlp = Arrayfun(pfun, Q_vlp);

Plot(Q_vlp, P_vlp, "k", 2);
SaveAs("VLP_Stimulation_Economics_Case_Study.png");

1. Nodal Analysis and Operating Points

We find the intersection where: math:p_{ wf, IPR} = p_{ wf, VLP} for each case using ** SepalSolver**.

//IPR Input
double q_max = 2000; // STB/day
double p_r = 2800; // psi

//IPR
double qfun(double pwf) => q_max * (1 - 0.2 * (pwf / p_r) - 0.8 * Pow(pwf / p_r, 2));
ColVec P_ipr = Linspace(0, p_r);
ColVec Q_ipr = Arrayfun(qfun, P_ipr);

// Adjusting for Damage and Stimulation
double r_e_r_w = 1000;
double FE_Damaged = Log(r_e_r_w)/(Log(r_e_r_w) + 5);
double FE_Stimulation_A = Log(r_e_r_w)/(Log(r_e_r_w) + -3);
double FE_Stimulation_B = Log(r_e_r_w)/(Log(r_e_r_w) + 1);
ColVec Q_ipr_Damaged = Q_ipr * FE_Damaged;
ColVec Q_ipr_Stimulated_A = Q_ipr * FE_Stimulation_A;
ColVec Q_ipr_Stimulated_B = Q_ipr * FE_Stimulation_B;



// VLP Inputs
double p_surf = 200; // psi (Wellhead Pressure)
double depth = 3500; // ft


// VLP
double pressureGradient(double z, double p, double q)
{
    double density = 48.0; // lb/ft3 (Oil)
    double friction_grad = 5e-6 * Pow(q, 1.8); // Simplified friction term
    double hydro_grad = density / 144.0; // psi/ft
    return hydro_grad + friction_grad;
}
double pfun(double q_o)
{
    var (Z, P) = Ode45((z, p) => pressureGradient(z, p, q_o), p_surf, [0, depth]);
    double p_wf = P[^1]; // extract the pressure at the bottom
    return p_wf;
}

ColVec Q_vlp = Linspace(0, 1600);
ColVec P_vlp = Arrayfun(pfun, Q_vlp);


var (current_q, current_p) = Intersection(Q_ipr_Damaged, P_ipr, Q_vlp, P_vlp);
var (stimulatedA_q, stimulatedA_p) = Intersection(Q_ipr_Stimulated_A, P_ipr, Q_vlp, P_vlp);
var (stimulatedB_q, stimulatedB_p) = Intersection(Q_ipr_Stimulated_B, P_ipr, Q_vlp, P_vlp);

Plot(Q_ipr_Damaged, P_ipr, "r", 2); HoldOn();
Plot(Q_ipr_Stimulated_A, P_ipr, "b", 2);
Plot(Q_ipr_Stimulated_B, P_ipr, "g", 2);
Plot(Q_vlp, P_vlp, "y", 2);
Scatter([current_q, stimulatedA_q, stimulatedB_q],
    [current_p, stimulatedA_p, stimulatedB_p], "fok"); HoldOff();
Axis([0, 3500, 0, 3500]);
Legend(["Damaged", "Stimulated_A", "Stimulated_B", "VLP"]);
SaveAs("Nodal_Analysis_Stimulation_Economics_Case_Study.png");

1. Production Improvement and Economics

//IPR Input
double q_max = 2000; // STB/day
double p_r = 2800; // psi

//IPR
double qfun(double pwf) => q_max * (1 - 0.2 * (pwf / p_r) - 0.8 * Pow(pwf / p_r, 2));
ColVec P_ipr = Linspace(0, p_r);
ColVec Q_ipr = Arrayfun(qfun, P_ipr);

// Adjusting for Damage and Stimulation
double r_e_r_w = 1000;
double FE_Damaged = Log(r_e_r_w)/(Log(r_e_r_w) + 5);
double FE_Stimulation_A = Log(r_e_r_w)/(Log(r_e_r_w) + -3);
double FE_Stimulation_B = Log(r_e_r_w)/(Log(r_e_r_w) + 1);
ColVec Q_ipr_Damaged = Q_ipr * FE_Damaged;
ColVec Q_ipr_Stimulated_A = Q_ipr * FE_Stimulation_A;
ColVec Q_ipr_Stimulated_B = Q_ipr * FE_Stimulation_B;



// VLP Inputs
double p_surf = 200; // psi (Wellhead Pressure)
double depth = 4000; // ft


// VLP
double pressureGradient(double z, double p, double q)
{
    double density = 52.0; // lb/ft3 (Oil)
    double friction_grad = 1e-5 * Pow(q, 1.8); // Simplified friction term
    double hydro_grad = density / 144.0; // psi/ft
    return hydro_grad + friction_grad;
}
double pfun(double q_o)
{
    var (Z, P) = Ode45((z, p) => pressureGradient(z, p, q_o), p_surf, [0, depth]);
    double p_wf = P[^1]; // extract the pressure at the bottom
    return p_wf;
}

ColVec Q_vlp = Linspace(0, 1600);
ColVec P_vlp = Arrayfun(pfun, Q_vlp);

var (current_q, current_p) = Intersection(Q_ipr_Damaged, P_ipr, Q_vlp, P_vlp);
var (stimulatedA_q, stimulatedA_p) = Intersection(Q_ipr_Stimulated_A, P_ipr, Q_vlp, P_vlp);
var (stimulatedB_q, stimulatedB_p) = Intersection(Q_ipr_Stimulated_B, P_ipr, Q_vlp, P_vlp);

double BarrelPerDollar_A = (stimulatedA_q - current_q)/10;
double BarrelPerDollar_B = (stimulatedB_q - current_q)/5;

Console.WriteLine($"Barrel Per Dollar for Quote A: {BarrelPerDollar_A}");
Console.WriteLine($"Barrel Per Dollar for Quote B: {BarrelPerDollar_B}");

Ouput

Barrel Per Dollar for Quote A: 3.5733706107198175
Barrel Per Dollar for Quote B: 3.4741949196266146