Model complete production well networks — from reservoir sandface through tubing, chokes, flowlines, and manifolds to delivery — using the LoopedPipeNetwork Newton-Raphson Global Gradient Algorithm (NR-GGA) solver.

Table of Contents

Overview

Location: neqsim.process.equipment.network

The LoopedPipeNetwork class extends the basic pipeline network solver with production-specific resistance elements. Instead of modeling every element as a pipe segment, you can now represent each component of the production system with its own physics:

Element Physics Typical Use
PIPE Darcy-Weisbach friction Flowlines, export lines, risers
WELL_IPR Inflow Performance Relationship Reservoir-to-wellbore inflow
CHOKE Valve flow equation (IEC 60534-style) Wellhead/surface chokes
TUBING Gravity + friction (VLP) Vertical tubing strings
MULTIPHASE_PIPE Segmented multiphase correlations Subsea flowlines
COMPRESSOR Gas pressure rise, power, and chart constraints Gas compression stations
PUMP Liquid pressure rise, pump curve, power, and NPSH Oil export boosters
REGULATOR Downstream pressure control Delivery pressure control

All element types participate in the same Newton-Raphson global gradient solver. The Schur complement reduction handles any mix of elements in a single network.

Network Element Types

Enum: NetworkElementType

public enum NetworkElementType {
    PIPE,             // Standard Darcy-Weisbach pipe
    WELL_IPR,         // Inflow Performance Relationship
    CHOKE,            // Valve / choke with Kv rating
    TUBING,           // Vertical tubing with gravity + friction
    MULTIPHASE_PIPE,  // Segmented multiphase pipe flow
    COMPRESSOR,       // Gas pressure-rise edge
    PUMP,             // Liquid pressure-rise edge
    REGULATOR         // Pressure-control edge
}

Each element type has its own head-loss and derivative calculation, allowing the NR-GGA solver to converge efficiently even when mixing very different resistance characteristics (e.g., a low-resistance IPR next to a high-resistance choke).

IPR Models

The Inflow Performance Relationship describes how much fluid flows from the reservoir into the wellbore as a function of drawdown $\Delta P = P_{res} - P_{wf}$.

Enum: IPRType

public enum IPRType {
    PRODUCTIVITY_INDEX,  // Linear PI (oil or gas)
    VOGEL,               // Vogel's equation for solution-gas drive
    FETKOVICH            // Fetkovich empirical model
}

Productivity Index (PI)

For oil (incompressible), the flow rate is:

\[Q = PI \times (P_{res} - P_{wf})\]

For gas (compressible), the flow is proportional to the pressure-squared difference:

\[Q = PI_{gas} \times (P_{res}^2 - P_{wf}^2)\]

Parameters:

Parameter Method Unit Description
Reservoir boundary pressure network.setReservoirPressure(sourceNodeName, pressure) bara Updates the source node and all connected IPR elements
IPR stored pressure NetworkPipe.setReservoirPressure(double) Pa Low-level element parameter; does not update the source-node boundary
Productivity index (oil) setProductivityIndex(double) kg/s/Pa Oil PI
Gas IPR flag setGasIPR(boolean) — Use pressure-squared formulation
Productivity index (gas) setProductivityIndex(double) kg/s/Pa² Gas PI

Convenience method:

network.addSourceNode("reservoir", 350.0, 0.0); // pressure [bara], supply [kg/h]
network.addJunctionNode("wellhead");
network.addWellIPR("reservoir", "wellhead", "ipr1", 5e-7, false);

The IPR inherits its reservoir pressure from the source. Use the network-level setter for depletion and the sensitivity API for sweeps; changing only the stored IPR pressure leaves the source boundary unchanged.

Vogel’s Equation

For solution-gas-drive reservoirs below the bubble point:

\[\frac{Q}{Q_{max}} = 1 - 0.2 \left(\frac{P_{wf}}{P_{res}}\right) - 0.8 \left(\frac{P_{wf}}{P_{res}}\right)^2\]

The absolute open-flow potential $Q_{max}$ is specified instead of PI.

Convenience method:

network.addWellIPRVogel("reservoir", "wellhead", "vogel1",
    50.0);           // Qmax [kg/s]

Fetkovich Empirical Model

A generalized IPR with empirical exponent $n$ (0.5 to 1.0):

\[Q = C \times (P_{res}^2 - P_{wf}^2)^n\]

Convenience method:

network.addWellIPRFetkovich("reservoir", "wellhead", "fetk1",
    1e-12,           // C coefficient
    0.8);            // n exponent (0.5–1.0)

Choke Model

Two choke modes are available. Select the model before solving and check isProductionOptimizationApplicable() before using its result in optimization.

Mode Supported use Critical-flow behavior
Default simplified equation Subcritical screening at the template fluid density Reports UNSUPPORTED_CRITICAL_FLOW; optimization rejects the point
setChokeUseValveModel(true) Newton-Raphson, steady single-phase gas, linear opening, no attached-fitting correction Solves the forward IEC 60534 gas capacity relation, retaining Kv/opening sensitivity

The legacy screening equation is \(\Delta P_{bar} = \left(\frac{3600\dot m}{\rho K_v\theta}\right)^2\) where $\dot m$ is kg/s, $\rho$ is kg/m³, and $\theta$ is opening/100. This historical relation omits the liquid specific-gravity factor and gas expansion; its coefficient must not be interpreted as a qualified installed valve Kv. Its pressure-drop cap is retained for compatibility, but a capped solution cannot establish valid production or an optimum. A closed screening choke is also unsupported; use the forward gas mode for shutoff.

Gas Capacity and Critical Flow

The opt-in mode uses the same forward relation as ThrottlingValve with the IEC 60534 sizing method. It flashes the physical upstream composition at the node pressure and temperature, then converts actual volumetric flow to kg/s. Set the upstream node temperature explicitly for the intended case.

Define the pressure-drop fraction $x=(P_1-P_2)/P_1$, with absolute pressures, and $x_c=(\gamma/1.4)x_T$. The sizing relation uses \(x_s=\min(x,x_c),\qquad Y=1-\frac{x_s}{3x_c}\) so $Y=2/3$ at and beyond critical flow. At fixed upstream conditions, reducing downstream pressure further leaves capacity unchanged; increasing effective Kv increases capacity. The solver enforces this flow limit and includes its pressure derivatives in the nodal mass-balance equations. It does not cap the physical pressure difference across the choke.

This formulation follows Emerson’s Control Valve Handbook, chapter 5, sections 5.9.1–5.9.4. The network regression compares a synthetic IPR → choke → arrival fixture against a separately configured ThrottlingValve. This qualifies the network coupling to the existing valve model; it is not an independent experimental qualification of a production well or a multiphase choke.

Parameter Method Value convention
Kv setChokeKv(double) Finite, nonnegative m³/h/√bar coefficient
Opening setChokeOpening(double) Percent, 0–100, with linear characterization
Critical ratio setChokeCriticalPressureRatio(double) Default 0.5; screening drop fraction or gas-mode IEC $x_T$, strictly between 0 and 1
Applicability getChokeModelStatus() Last evaluated model/regime status
Gas capacity getChokeCapacityKgS() Signed kg/s at the last evaluated pressure pair 
network.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
network.getNode("wellhead").setTemperature(313.15); // K
LoopedPipeNetwork.NetworkPipe choke =
    network.addChoke("wellhead", "arrival", "choke", 10.0, 60.0);
choke.setChokeCriticalPressureRatio(0.5);
choke.setChokeUseValveModel(true);
network.run();
if (!network.isProductionOptimizationApplicable()) {
  throw new IllegalStateException(choke.getChokeModelStatus());
}

The gas mode supports zero/full opening and reverse flow, using the physical upstream node. A liquid or multiphase inlet, invalid pressure, or a solver other than Newton-Raphson raises an explicit error. It does not silently fall back to screening. Set actual trim Kv and $x_T$ from applicable valve data; constant $x_T$ and linear travel are assumptions, not universal valve characteristics. Availability scales effective Kv; zero availability closes the gas choke. Apply artificial lift to a separate network element.

optimizeFullField() returns chokeModelApplicable and a per-choke chokeModelStatus map. An unsupported screening point has converged = 0 and NaN objective/revenue even if isConverged() reports an algebraic solution. Start choke searches from an applicable point. Topside limits and other user-defined constraints must still be checked separately.

Tubing VLP Model

The tubing element models vertical lift performance (VLP) — the pressure drop in the production tubing due to gravity and friction:

\[\Delta P = \rho g L \sin\alpha + f \frac{L}{D} \frac{\rho v^2}{2}\]

Where:

For deep wells, the tubing is divided into segments, each evaluated at local conditions (pressure-averaged density).

Parameters:

Parameter Setter on NetworkPipe Unit Description
Inclination setTubingInclination(double) degrees From horizontal
Segments setTubingSegments(int) — Number of integration segments

Convenience method:

network.addTubing("tubing1", "bottomhole", "wellhead",
    3000.0,          // measured depth [m]
    0.1,             // tubing ID [m]
    90.0,            // inclination [degrees from horizontal]
    10);             // number of segments

Multiphase Pipe Flow

Multiphase pipe elements use a segmented approach similar to tubing but oriented horizontally (or at any angle). Each segment evaluates local friction and gravity based on bulk fluid properties.

Convenience method:

network.addMultiphasePipe("manifold", "platform", "mp1",
    15000.0,         // length [m]
    0.3);            // ID [m]

Pipe roughness is specified in metres throughout the network and PipeBeggsAndBrills APIs. For example, commercial-steel roughness can be set with network.getPipe("mp1").setRoughness(4.5e-5).

To use Beggs–Brill for ordinary pipes added afterwards, configure the network-wide default before calling addPipe:

network.setPipeModelType(LoopedPipeNetwork.PipeModelType.BEGGS_BRILL);
network.addPipe("manifold", "platform", "mp1", 15000.0, 0.3);

The setting applies only to later addPipe calls. It does not change existing pipes or elements created with specialized methods such as addChoke and addTubing. addMultiphasePipe always selects Beggs–Brill explicitly.

Building a Production Network

Step-by-Step Pattern

// 1. Create network and set fluid template
LoopedPipeNetwork net = new LoopedPipeNetwork("Gathering System");
SystemSrkEos gas = new SystemSrkEos(288.15, 100.0);
gas.addComponent("methane", 0.85);
gas.addComponent("ethane", 0.10);
gas.addComponent("propane", 0.05);
gas.setMixingRule("classic");
net.setFluidTemplate(gas);

// 2. Add reservoir nodes (SOURCE, fixed pressure)
net.addSourceNode("res1", 350.0, 0.0);  // P=350 bara, demand=0

// 3. Add delivery node (SINK, fixed pressure)
net.addFixedPressureSinkNode("separator", 50.0);

// 4. Add intermediate junction nodes
net.addJunctionNode("wellhead");
net.addJunctionNode("downstream_choke");

// 5. Add production elements
net.addWellIPR("ipr", "res1", "wellhead", 350.0, 5e-7, false);
net.addChoke("choke", "wellhead", "downstream_choke", 150.0, 0.8, 0.7);
net.addPipe("downstream_choke", "separator", "flowline", 5000.0, 0.2);

// 6. Configure solver
net.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
net.setTolerance(500.0);       // Pa
net.setMaxIterations(500);

// 7. Solve
net.run();

// 8. Inspect results
System.out.println("Converged: " + net.isConverged());
System.out.println("Well rate: " + net.getPipeFlowRate("ipr") + " kg/hr");
System.out.println("WHP: " + net.getNodePressure("wellhead") + " bara");

Multi-Well Gathering Network

     [Res 1]            [Res 2]
        │                   │
    (IPR 1)             (IPR 2)
        │                   │
     [WH 1]             [WH 2]
        │                   │
    (Choke 1)           (Choke 2)
        │                   │
     [DS 1]             [DS 2]
        │                   │
    (Pipe 1)            (Pipe 2)
        │                   │
        └───────┬───────────┘
                │
           [Manifold]
                │
           (Export Pipe)
                │
          [Separator]

LoopedPipeNetwork net = new LoopedPipeNetwork("Multi-Well");
net.setFluidTemplate(gas);

// Well 1: high-pressure, high-PI
net.addSourceNode("res1", 350.0, 0.0);
net.addJunctionNode("wh1");
net.addJunctionNode("ds1");
net.addWellIPR("ipr1", "res1", "wh1", 350.0, 8e-13, true);   // gas IPR
net.addChoke("choke1", "wh1", "ds1", 150.0, 0.8, 0.7);
net.addPipe("ds1", "manifold", "pipe1", 5000.0, 0.2);

// Well 2: lower-pressure, lower-PI
net.addSourceNode("res2", 280.0, 0.0);
net.addJunctionNode("wh2");
net.addJunctionNode("ds2");
net.addWellIPR("ipr2", "res2", "wh2", 280.0, 5e-13, true);
net.addChoke("choke2", "wh2", "ds2", 150.0, 0.6, 0.7);
net.addPipe("ds2", "manifold", "pipe2", 8000.0, 0.2);

// Gathering
net.addJunctionNode("manifold");
net.addFixedPressureSinkNode("separator", 50.0);
net.addPipe("manifold", "separator", "export", 10000.0, 0.3);

// Solve
net.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
net.setTolerance(500.0);
net.setMaxIterations(500);
net.run();

Transmission and Terminal Extensions

The production-network solver also supports heterogeneous transmission and terminal studies:

Legacy hydraulics use one fluidTemplate for all edges. Enable coupled mode only when source composition, condensation, compression temperature, or heat transfer can materially change allocation:

network.setCompositionalHydraulicsEnabled(true);
network.setThermalHydraulicsEnabled(true);
network.setCouplingMaxIterations(20);
network.setCouplingTolerances(1.0e-5, 100.0);
network.run();

NetworkCouplingReport report =
    network.getNetworkCouplingReport();

At a mixing node, component molar flows are conserved and aligned by component name. The mixed fluid is pressure/enthalpy flashed and initialized. Loops use a synchronous fixed-point iteration. Incompatible EOS classes, mixing rules, or pseudo-component definitions fail explicitly.

Use the dedicated guides for units, assumptions, API examples, JSON behavior, and limitations:

The NCS-style examples use synthetic compositions, capacities, quality limits, and nominations. Current operator data, approved assays, commercial agreements, and validated asset models remain controlling.

Solver Details

Newton-Raphson Global Gradient Algorithm

The production network solver uses the NR-GGA (Todini & Pilati, 1988) extended for generalized resistance elements. The system is formulated as:

\[\begin{bmatrix} D & -A^T \\ A & 0 \end{bmatrix} \begin{bmatrix} \delta Q \\ \delta H \end{bmatrix} = \begin{bmatrix} -f_1 \\ -f_2 \end{bmatrix}\]

Where:

For pipe elements, $f_1$ uses the Darcy-Weisbach resistance inline:

\[f_1 = D_i Q_i + (H_{from} - H_{to})\]

For production elements (IPR, choke, tubing), $f_1$ compares the element-specific head loss to the available pressure difference:

\[f_1 = h_{element}(Q_i) - (P_{from} - P_{to})\]

The Schur complement reduces the $(np + nn) \times (np + nn)$ system to an $nn \times nn$ system for node heads, then back-substitutes for flows.

Adaptive Relaxation

The solver uses adaptive under-relaxation for stability:

Iteration Range Residual Condition Relaxation Factor
1–5 $> 10^6$ Pa 0.3
6–15 $> 10^4$ Pa 0.6
After 15 Any 1.0 (full Newton)

Flow Initialization

Production elements get physics-based initial flow estimates:

Pressures at free nodes are initialized by BFS propagation from fixed-pressure nodes (sources and sinks).

Usage Examples

Example 1: Single Well with PI

LoopedPipeNetwork net = new LoopedPipeNetwork("Single Well PI");
net.setFluidTemplate(gas);

net.addSourceNode("reservoir", 350.0, 0.0);
net.addJunctionNode("wellhead");
net.addFixedPressureSinkNode("separator", 40.0);

net.addWellIPR("ipr", "reservoir", "wellhead", 350.0, 5e-7, false);
net.addPipe("wellhead", "separator", "flowline", 5000.0, 0.2);

net.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
net.setTolerance(100.0);
net.setMaxIterations(200);
net.run();

double wellRate = net.getPipeFlowRate("ipr");     // kg/hr
double whp = net.getNodePressure("wellhead");      // bara

Example 2: Vogel IPR

net.addWellIPRVogel("vogel", "reservoir", "wellhead", 350.0, 50.0);

Example 3: Complete Well System

IPR + tubing + choke + flowline in series:

LoopedPipeNetwork net = new LoopedPipeNetwork("Complete Well");
net.setFluidTemplate(gas);

net.addSourceNode("reservoir", 400.0, 0.0);
net.addJunctionNode("bottomhole");
net.addJunctionNode("wellhead");
net.addJunctionNode("ds_choke");
net.addFixedPressureSinkNode("separator", 40.0);

net.addWellIPR("ipr", "reservoir", "bottomhole", 400.0, 5e-7, false);
net.addTubing("tubing", "bottomhole", "wellhead", 3000.0, 0.1, 90.0, 10);
net.addChoke("choke", "wellhead", "ds_choke", 150.0, 0.7, 0.7);
net.addPipe("ds_choke", "separator", "flowline", 5000.0, 0.2);

net.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
net.setTolerance(500.0);
net.setMaxIterations(500);
net.run();

Example 4: Choke Sensitivity Study

Sweep choke opening from 10% to 100% and observe how production rate and wellhead pressure respond:

double[] openings = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
for (double opening : openings) {
    net.getNetworkPipe("choke").setChokeOpening(opening);
    net.run();
    System.out.printf("Opening=%.0f%%  Q=%.0f kg/hr  WHP=%.1f bara%n",
        opening * 100,
        net.getPipeFlowRate("ipr"),
        net.getNodePressure("wellhead"));
}

Results and Inspection

Solution Summary

// JSON-like summary of the network solution
String summary = net.getSolutionSummary();

Returns a map with keys:

Element Inspection

// Inspect each element's type and calculated values
for (String pipeName : net.getPipeNames()) {
    NetworkPipe pipe = net.getNetworkPipe(pipeName);
    System.out.printf("%-12s type=%-18s flow=%10.0f kg/hr%n",
        pipeName,
        pipe.getElementType(),
        net.getPipeFlowRate(pipeName));
}

Hydraulic Report

// Detailed pipe-by-pipe hydraulic report
String report = net.getHydraulicReport();

Mass Balance

// Verify mass balance at every node
String massBalance = net.getMassBalanceReport();

Artificial Lift

The LoopedPipeNetwork supports four artificial lift methods that inject additional energy into wells to increase production. Each method applies a pressure boost in the NR-GGA solver’s calculateHeadLoss() step.

Lift Types

Type Enum Value Method Typical Application
Gas lift GAS_LIFT setGasLift() High GOR wells, offshore
ESP ESP setESP() High-rate wells, deepwater
Jet pump JET_PUMP setJetPump() Low-rate, deviated wells
Rod pump ROD_PUMP setRodPump() Onshore, low-rate stripper wells

API

// Gas lift: inject gas to reduce hydrostatic gradient
net.setGasLift("W3", 500.0);              // 500 kg/hr gas lift rate

// ESP: electrical submersible pump
net.setESP("W4", 80.0, 0.55);             // 80 kW power, 55% efficiency

// Jet pump: Venturi-based downhole pump
net.setJetPump("W5", 600.0, 0.40);        // 600 kW nozzle power, 40% efficiency

// Rod pump: sucker-rod beam pump
net.setRodPump("W6", 30.0, 0.35);         // 30 kW motor, 35% efficiency

// Query totals
double totalGL = net.getTotalGasLiftRate();   // kg/hr
double totalESP = net.getTotalESPPower();     // kW

Effect on Solver

Each artificial lift type converts its input energy (gas injection rate or electrical power) to a pressure boost $\Delta P_{lift}$ that is subtracted from the well element’s head loss. For example, gas lift reduces the effective hydrostatic gradient, while ESP/jet/rod pump add a direct pressure rise to the well’s IPR resistance.

Water Handling

Track produced water, water injection, and water breakthrough per well.

API

// Set water cut on a well (fraction, 0-1)
net.setWaterCut("W1", 0.15);       // 15% water cut
net.setWaterCut("W2", 0.40);       // 40% water cut

// Add water injection support
// addWaterInjection(sourceNode, reservoirNode, elementName, rateKgHr)
net.addWaterInjection("WI_supply", "R2", "WI-W2", 2000.0);

// Water breakthrough tracking
// setWaterBreakthrough(elementName, breakthroughWC, finalWC, currentWC)
net.setWaterBreakthrough("W1", 0.10, 0.65, 0.10);

// Calculate water balance (after run)
Map<String, double[]> wb = net.calculateWaterBalance();
// Per node: [0] = water production (kg/hr), [1] = injection (kg/hr), [2] = net

double totalProd = net.getTotalWaterProduction();   // kg/hr
double totalInj  = net.getTotalWaterInjection();    // kg/hr

Sand and Solids Tracking

Per-element sand production, erosion rate (DNV RP O501), and deposition tracking with configurable limits and violation reporting.

API

net.setSandRate("W1", 3.0);             // 3 kg/hr sand production
net.setSandRate("W2", 12.0);            // 12 kg/hr

net.setMaxAllowableSandRate(10.0);      // kg/hr limit
net.setMaxAllowableErosionRate(5.0);    // mm/yr limit

net.run();
Map<String, double[]> sand = net.calculateSandTransport();
// Per element: [0] = sand rate (kg/hr), [1] = concentration (kg/m3),
//              [2] = erosion rate (mm/yr), [3] = deposition rate (mm/yr)

List<String> violations = net.getSandViolations();

Erosion Model (DNV RP O501)

\[E = K \cdot C_{sand} \cdot v^{2.6} \cdot d_p^{0.2} \cdot 3600 \cdot 8760\]

where $E$ is erosion rate (mm/yr), $K$ is the material factor, $C_{sand}$ is sand concentration (kg/m3), $v$ is flow velocity (m/s), and $d_p$ is particle diameter (m). Deposition is flagged when velocity falls below 1 m/s.

Corrosion and Integrity

Inline CO2/H2S corrosion rate estimation with two industry-standard models and remaining wall life calculation.

API

// Set corrosive gas composition on a pipe or well
net.setCorrosiveGas("trunk", 0.035, 0.002);  // 3.5 mol% CO2, 0.2 mol% H2S
net.setCorrosionModel("trunk", "NORSOK");     // "DEWAARD" (default) or "NORSOK"

net.setMinAllowableWallLife(20.0);  // years

net.run();
Map<String, double[]> corr = net.calculateCorrosion();
// Per element: [0] = corrosion rate (mm/yr), [1] = pCO2 (bar),
//              [2] = remaining wall life (yr)

List<String> violations = net.getCorrosionViolations();

de Waard-Milliams Model

\[\log_{10}(V_{corr}) = 5.8 - \frac{1710}{T} + 0.67 \cdot \log_{10}(p_{CO_2})\]

where $T$ is temperature (K) and $p_{CO_2}$ is CO2 partial pressure (bar). An H2S multiplier $(1 + 2 \cdot y_{H_2S}/y_{CO_2})$ is applied when H2S is present.

NORSOK M-506 Model

\[V_{corr} = K_t \cdot f_{CO_2}^{0.62} \cdot \left(\frac{S}{19}\right)^{0.146}\]

where $K_t$ is a temperature-dependent factor, $f_{CO_2}$ is CO2 fugacity, and $S$ is the wall shear stress (Pa). The NORSOK model generally predicts lower rates than de Waard-Milliams and accounts for protective scale formation.

GHG Emissions Tracking

Calculate CO2 and methane emissions per compressor station with field-level totals, annual figures, and emissions intensity.

API

net.setCO2EmissionFactor(2.75);     // kg CO2 per kg fuel gas (default)
net.setMethaneSlipFactor(0.02);     // 2% methane slip (default)

net.run();

// Set compressor power (if not computed by the solver)
net.getPipe("booster").setCompressorPower(1500.0);  // kW

Map<String, double[]> em = net.calculateEmissions();
// Per compressor: [0] = CO2 (kg/hr), [1] = CH4 slip (kg/hr),
//                 [2] = CO2-eq (kg/hr), [3] = power (kW), [4] = fuel gas (kg/hr)

double totalCO2eq = net.getTotalCO2Emissions();        // kg/hr
double annual     = net.getAnnualCO2EmissionsTonnes();  // tonnes/yr
double intensity  = net.getEmissionsIntensity();        // kgCO2-eq/tonne product

Emission Factors

Parameter Default Units Source
CO2 emission factor 2.75 kgCO2/kg fuel Natural gas combustion
Methane slip 2% fraction Gas turbine typical
Methane GWP 28 - IPCC AR5 (100-yr)
Fuel gas heat rate 10,000 kJ/kWh Gas turbine typical
LHV natural gas 48,000 kJ/kg -

Calculations

\[\dot{m}_{fuel} = \frac{P \cdot HR}{LHV \cdot 3600}\] \[CO_{2,eq} = \dot{m}_{fuel} \cdot EF_{CO_2} + \dot{m}_{fuel} \cdot f_{slip} \cdot GWP_{CH_4}\]

where $P$ is compressor power (kW), $HR$ is heat rate (kJ/kWh), $LHV$ is LHV (kJ/kg), $EF_{CO_2}$ is the CO2 emission factor, $f_{slip}$ is methane slip fraction, and $GWP_{CH_4}$ is the methane global warming potential.

Large-Scale Networks

The NR-GGA solver handles networks with 100+ wells efficiently. The Schur complement reduction keeps the system matrix size proportional to the number of independent loops, not the total number of elements.

Scaling Example

LoopedPipeNetwork net = new LoopedPipeNetwork("large");
net.setFluidTemplate(gas);
net.setSolverType(LoopedPipeNetwork.SolverType.NEWTON_RAPHSON);
net.setMaxIterations(500);
net.setTolerance(500.0);

// 6 manifolds, 20 wells each = 120 wells
for (int m = 0; m < 6; m++) {
    net.addJunctionNode("MF" + (m + 1));
    for (int w = 0; w < 20; w++) {
        double resP = 230.0 + random(-30, 30);
        double pi = 5e-13 * random(0.5, 2.0);
        net.addSourceNode("R" + m + "_" + w, resP, 0.0);
        net.addWellIPR("R" + m + "_" + w, "MF" + (m + 1),
                        "W" + m + "_" + w, pi, true);
    }
}

// Star topology: manifolds to hub to export
net.addJunctionNode("hub");
for (int m = 0; m < 6; m++) {
    net.addPipe("MF" + (m + 1), "hub", "trunk" + (m + 1),
                5000.0 + m * 3000.0, 0.3);
}
net.addFixedPressureSinkNode("export", 70.0);
net.addPipe("hub", "export", "export_line", 30000.0, 0.4);

net.run();
// Typically converges in 15-20 iterations, < 0.1 s

Performance Benchmarks

Network Size Wells Elements Iterations Time
Small 4 10 8-12 < 0.01 s
Medium 20 50 12-15 < 0.03 s
Large 120 300+ 15-20 < 0.1 s

Troubleshooting

Common Issues

Symptom Likely Cause Fix
Non-convergence with choke Kv too small for the flow Increase Kv or reduce flow rate
Unsupported critical screening point Legacy pressure-drop cap has no flow-capacity relation Use the supported single-phase gas valve mode or retain a subcritical screening case
Zero flow from IPR Separator pressure > reservoir pressure Check node pressures
Negative flow in tubing Wrong node ordering Ensure from-node is bottomhole
Slow convergence Large pressure span in network Increase maxIterations to 500

Choosing Tolerances

For pipe-only networks (water distribution, gas transmission):

For production networks with chokes and IPR elements:

The larger tolerance is needed because choke elements have steep resistance curves that amplify small flow changes into large pressure residuals.

API Reference

Convenience Methods on LoopedPipeNetwork

// Default model for ordinary pipes added afterwards
void setPipeModelType(PipeModelType type)

// Ordinary pipe using the configured default model
NetworkPipe addPipe(String fromNode, String toNode, String pipeName,
    double length_m, double diameter_m)

// Well IPR (oil or gas PI model)
void addWellIPR(String name, String fromNode, String toNode,
    double reservoirPressure_bara, double pi, boolean gasIPR)

// Vogel IPR (solution-gas drive)
void addWellIPRVogel(String name, String fromNode, String toNode,
    double reservoirPressure_bara, double qmax_kgs)

// Fetkovich empirical IPR
void addWellIPRFetkovich(String name, String fromNode, String toNode,
    double reservoirPressure_bara, double C, double n)

// Choke / valve
void addChoke(String name, String fromNode, String toNode,
    double kv, double opening, double criticalPressureRatio)

// Tubing (vertical lift performance)
void addTubing(String name, String fromNode, String toNode,
    double length_m, double diameter_m,
    double inclination_deg, int segments)

// Multiphase pipe using Beggs-Brill explicitly
NetworkPipe addMultiphasePipe(String fromNode, String toNode,
    String elementName, double length_m, double diameter_m)

// Liquid pump / booster
NetworkPipe addPump(String fromNode, String toNode, String name,
    double outletPressure_bara, double efficiency)
NetworkPipe addPumpDifferentialPressure(
    String fromNode, String toNode, String name,
    double differentialPressure_bar, double efficiency)
NetworkPipe addPumpWithCurve(
    String fromNode, String toNode, String name, Pump pump)

// Fixed-pressure sink (separator, platform)
void addFixedPressureSinkNode(String name, double pressure_bara)

NetworkPipe Production Fields

// Element type
NetworkElementType getElementType()
void setElementType(NetworkElementType type)

// IPR fields
IPRType getIprType()
void setIprType(IPRType type)
double getReservoirPressure()
void setReservoirPressure(double p)
double getProductivityIndex()
void setProductivityIndex(double pi)
boolean isGasIPR()
void setGasIPR(boolean gas)
double getVogelQmax()
void setVogelQmax(double qmax)
double getFetkovichC()
void setFetkovichC(double c)
double getFetkovichN()
void setFetkovichN(double n)

// Choke fields
double getChokeKv()
void setChokeKv(double kv)
double getChokeOpening()
void setChokeOpening(double opening)
double getChokeCriticalPressureRatio()
void setChokeCriticalPressureRatio(double ratio)

// Tubing / multiphase fields
double getTubingInclination()
void setTubingInclination(double degrees)
int getTubingSegments()
void setTubingSegments(int n)
int getMultiphaseSegments()
void setMultiphaseSegments(int n)

// Pump fields and constraint residuals
PumpOperatingMode getPumpOperatingMode()
void setPumpOperatingMode(PumpOperatingMode mode)
double getPumpSpeed()
void setPumpSpeed(double rpm)
double getPumpPowerKW()
double getPumpHeadM()
double getPumpPowerResidualKW()
double getPumpMinimumFlowResidualKgS()
double getPumpNpshResidualM()
String getPumpOperatingStatus()

NetworkElementType Enum

public enum NetworkElementType {
    PIPE,
    WELL_IPR,
    CHOKE,
    TUBING,
    MULTIPHASE_PIPE,
    COMPRESSOR,
    PUMP,
    REGULATOR
}

IPRType Enum

public enum IPRType {
    PRODUCTIVITY_INDEX,
    VOGEL,
    FETKOVICH
}

References

  1. Todini, E. & Pilati, S. (1988). “A gradient algorithm for the analysis of pipe networks.” Computer Applications in Water Supply, Vol. 1, pp. 1–20.
  2. Vogel, J.V. (1968). “Inflow Performance Relationships for Solution-Gas Drive Wells.” JPT, 83–92.
  3. Fetkovich, M.J. (1973). “The Isochronal Testing of Oil Wells.” SPE 4529.
  4. IEC 60534-2-1. “Industrial-process control valves — Flow capacity.”