TwoFluidPipe Model Documentation

Overview

The NeqSim TwoFluidPipe model implements a transient two-fluid multiphase flow solver for pipeline and riser simulations. It solves phase-resolved gas, hydrocarbon-liquid, and aqueous-liquid conservation equations, enabling prediction of:

This document provides comprehensive documentation of the model’s capabilities, governing equations, and usage.

The selectable closure sets are literature-inspired NeqSim implementations. Historical API names containing OLGA are retained for compatibility and do not claim numerical equivalence with OLGA, LedaFlow, or another commercial simulator.

Conservation Equations

Mass Conservation

Separate mass conservation equations are solved for gas, hydrocarbon liquid, and aqueous liquid:

Equation Mathematical Form Description
Gas mass ∂(A αG ρG)/∂t + ∂(A αG ρG vG)/∂x = ΓG Gas phase continuity with mass transfer
Oil mass ∂(A αO ρO)/∂t + ∂(A αO ρO vO)/∂x = ΓO Hydrocarbon-liquid continuity
Water mass ∂(A αW ρW)/∂t + ∂(A αW ρW vW)/∂x = ΓW Aqueous-liquid continuity
Phase transfer ΓG + ΓO + ΓW = 0 Flash-based transfer with inventory limits

Where:

Flash-driven phase identity

ThermodynamicCoupling identifies phases by PhaseType; it does not assume gas, oil, or aqueous phases occupy fixed array positions. A gas + oil + aqueous PT flash aggregates both liquid phases in the equilibrium-liquid target. For condensation, the new liquid is split using equilibrium oil and aqueous mass contributions. The current hydrodynamic water cut is intentionally not used at phase appearance because a gas-only cell contains no information about the identity of its first liquid.

For evaporation, oil and water withdrawals are distributed from the actual conservative phase inventories. Each withdrawal is bounded by phase mass / relaxation time, so an absent phase cannot evaporate and no phase can be removed faster than the relaxation step permits. The immutable PhaseMassTransfer result reports gas, oil, and water sources in kg/(m s), together with flash convergence and applicability metadata.

Transferred momentum uses donor velocity: condensing gas gives each receiving liquid gas momentum, while evaporating oil and water give the gas their respective liquid momenta. Transfer-only gas, oil, and water momentum sources therefore sum to zero. This closure preserves phase and total mass and mixture momentum. Full named-component transport is available through the opt-in conservative component mode described below.

FlashTable stores the same aggregate liquid fraction, oil/aqueous liquid mass split, and gas/liquid molar masses as the rigorous flash path. Interpolated liquid identity fractions are clamped and renormalized; an identity that is absent at all surrounding grid points remains exactly zero.

Conservative named-component transport

Call setComponentTransportEnabled(true) before run() to initialize independent component inventories in every physical cell and each gas, oil, and aqueous phase. The default remains false for backward compatibility.

For cell $i$, phase $k$, and named component $c$, the conserved inventory is $M_{i,k,c}$ [kg]. For an accepted hydrodynamic substep, the component face flux is

\[F_{f,k,c} = \dot m_{f,k}Y_{upwind(f,k),k,c},\]

where $\dot m_{f,k}$ is the integration-weighted phase mass flux already used by the phase continuity equation. The upwind state follows the sign of the internal face flux. The inlet stream supplies positive-flow inlet composition. Integrated boundary component masses use the same Euler, Runge-Kutta, or IMEX stage weights as the accepted phase update.

Flash-driven phase transfer is mapped by component name. Evaporation withdraws the donor oil or water composition. Condensation uses the receiving equilibrium phase composition. In every cell,

\[\sum_{k\in\{G,O,W\}} \Delta M_{i,k,c}^{transfer}=0\]

for every component. After transport, the component sums must match the accepted hydrodynamic phase inventories within componentConservationTolerance; the implementation only removes floating-point round-off and throws if a material mismatch would require a projection. Cell PT flashes are then built from total conservative named-component inventories. These flashes update density, viscosity, sound speed, phase identity, and phase enthalpy but never overwrite the conserved inventories.

The interphase energy closure evaluates phase-specific partial component enthalpies at the cell pressure and temperature. For one accepted transfer,

\[Q_{latent,i}=-\sum_k\sum_c \Delta M_{i,k,c}^{transfer}\bar h_{i,k,c}.\]

Positive $Q_{latent}$ is heat released into the fluid sensible-energy equation; negative values consume sensible energy. This term is exposed in both TwoFluidComponentConservationReport.getInterphaseLatentHeatEnergyJ() and TwoFluidThermalEnergyBalanceReport.getLatentHeatEnergyJ(). It is included exactly once in the thermal residual.

pipe.setComponentTransportEnabled(true);
pipe.setComponentConservationTolerance(1.0e-8);
pipe.setStoreComponentConservationHistory(true);
pipe.run();

pipe.runTransient(0.1, UUID.randomUUID());
TwoFluidComponentConservationReport components =
    pipe.getLastComponentConservationReport();
double[] co2Gas = pipe.getComponentMassFractionProfile(
    TwoFluidComponentConservationReport.Phase.GAS, "CO2");
double outletCo2 = pipe.getOutletComponentMassFraction(
    TwoFluidComponentConservationReport.Phase.GAS, "CO2");
String json = components.toJson();

TwoFluidComponentConservationHistory retains time-aligned immutable reports when history storage is enabled. All array getters return defensive copies. The same APIs are directly accessible from Python through JPype; devtools/neqsim_dev_setup.py imports the pipe, reports, history, and phase enum into the standard notebook namespace. The public report constructor rejects missing or duplicate component names, inconsistent phase/component/cell dimensions, and non-finite diagnostic values before an invalid report can cross the Java or JPype boundary.

Validated scope and fail-loud boundaries

The one-phase gas limit is regression-tested against the validated PipeFlowSystem conservative species pulse. On a compact 2 m case, coarse and jointly refined grids/timesteps must keep the outlet nitrogen mass fraction within 0.08 absolute mass fraction of the reference; refinement may not move away by more than 0.01. This stated margin accounts for different hydraulic grids and explicit versus implicit first-order transport. A closed SRK-CPA water-dew-point transition separately checks every component, cell-wise equal/opposite phase transfer, phase mass, boundedness, deterministic history, and the latent-inclusive thermal residual. Engineering applications should repeat mesh and time-step refinement at their own length, velocity, phase split, and event duration.

Momentum Conservation

Separate momentum equations for each phase:

Component Implementation
Gas momentum Full 1D momentum with wall shear, interfacial shear, pressure gradient
Liquid momentum Full 1D momentum with wall shear, interfacial shear, pressure gradient
Wall friction Pipe roughness-based (Colebrook/Blasius correlations)
Interfacial friction Flow-regime dependent correlations

Optional stiff dispersed-bubble drag

setEnableStiffBubbleDrag(true) opts into the dimensionally correct Schiller-Naumann dispersed-bubble force and a conservative local implicit source solve. For a spherical bubble population,

\[F_i=\frac{3}{4} C_D \rho_L \alpha_G \frac{A}{d_b} (v_G-v_L)|v_G-v_L|, \qquad a_i=\frac{6\alpha_G}{d_b}, \qquad f_i=\frac{C_D}{4}.\]

Bubble and dispersed-bubble regimes use the explicit algebraic diameter closure

d_b = min(2 sqrt(0.725 sigma_b / (g |rho_L-rho_G|)), f_D D).

The defaults $\sigma_b=0.02$ N/m and $f_D=0.20$ preserve the historical calculation. Change the fixed values, or explicitly use each section’s thermodynamic phase-property surface tension, through the public pipe API:

pipe.setBubbleSurfaceTension(0.025);
pipe.setMaximumBubbleDiameterFraction(0.15);
pipe.setUseLocalBubbleSurfaceTension(true);

Local mode uses the surface tension already stored for each section by the thermodynamic coupling; the default remains fixed for compatibility. This single algebraic scale does not model a bubble-size distribution, deformation, coalescence, breakup, or turbulent-dissipation dependence.

The source operator solves the active gas and combined-liquid momenta by backward Euler and applies one half-step on each side of the transport update. It conserves total active-phase momentum to roundoff, decreases slip and kinetic energy, removes exactly absent phases instead of applying a mass floor, and partitions the liquid impulse by active oil/water mass so existing oil-water slip is preserved. The source evaluation is local and retains no stage history.

This mode is opt-in for migration compatibility. The legacy force scaling remains the default because the corrected closure, although numerically stable, is not yet quantitatively validated by the public Tengesdal severe-slugging benchmark. That comparison was made when the benchmark still asserted a riser-head-scaled pressure swing, which has since been shown to come from a saturated minimum-slip bound rather than the momentum balance, so the recorded pass counts predate the rebased acceptance bounds and the comparison has to be repeated before it means anything. These results indicate a remaining closure or regime-transition limitation rather than a stiff-source instability.

Optional virtual-mass coupling

setEnableVirtualMassForce(true) enables a local added-inertia coupling between gas and combined liquid momentum. After the complete uncoupled finite-volume right-hand side is assembled, the model computes

\[K=C_{vm}\alpha_G\rho_L A, \qquad F_{vm,G}=\frac{-K(a_{G,0}-a_{L,0})} {1+K(1/m_G+1/m_L)}, \qquad F_{vm,L}=-F_{vm,G}.\]

Here $a_{k,0}=(d(m_kv_k)/dt-v_kdm_k/dt)/m_k$ is obtained from the current integration-stage state and its complete uncoupled rate. The operator retains no velocity history, so repeated RHS evaluations and rejected integrator stages cannot change a later evaluation. For gas-oil-water flow, the combined liquid correction is partitioned between oil and water by conservative liquid mass, preserving mixture momentum without creating oil-water transfer. Coupling tends continuously to zero when either phase inventory is absent.

The default $C_{vm}=0.5$ is the spherical-bubble value. It is not a universal calibration for slug, churn, or annular flow, and enabling it does not establish accuracy parity with a commercial multiphase simulator. Validate the coefficient and transient response against applicable public or project data before engineering use.

Energy Conservation

Feature Description
Mixture energy equation Full energy balance including kinetic and potential terms
Joule-Thomson effect Enabled by default for accurate temperature prediction
Multi-layer heat transfer RadialThermalLayer and MultilayerThermalCalculator classes

Flow Regime Detection

Gas-Liquid Flow Regimes

The gas-liquid flow regime detector uses Taitel-Dukler transitions:

Regime Detection Criteria Status
STRATIFIED_SMOOTH Low gas velocity, stable interface ✅
STRATIFIED_WAVY Kelvin-Helmholtz instability criterion ✅
SLUG Liquid bridging criterion ✅
ANNULAR Weber number > 30 ✅
CHURN Transition between slug and annular ✅
BUBBLE High liquid fraction, low gas velocity ✅

Oil-Water Flow Regime Detection

For three-phase (gas-oil-water) simulations the OilWaterFlowRegimeDetector classifies the liquid-phase configuration at every pipe section. This is critical for corrosion prediction (water wetting), effective viscosity calculation, and water dropout risk assessment.

Based on Trallero (1995), Brauner (2003), and Angeli & Hewitt (2000):

Regime Condition Description
STRATIFIED $v_m < 0.1\,v_{crit}$ Separate oil and water layers
STRATIFIED_WITH_MIXING $0.1\,v_{crit} < v_m < 0.5\,v_{crit}$ Stratified with interfacial mixing zone
DISPERSED_OIL_IN_WATER $v_m > v_{crit}$ and $w_c > w_{inv}$ Oil droplets in continuous water
DISPERSED_WATER_IN_OIL $v_m > v_{crit}$ and $w_c < w_{inv}$ Water droplets in continuous oil
DUAL_DISPERSION $v_m \approx v_{crit}$ and $w_c \approx w_{inv}$ Both O/W and W/O regions coexist
ANNULAR High velocity, large density difference Oil core with water annulus or vice versa
SINGLE_PHASE $w_c < 0.005$ or $w_c > 0.995$ Only oil or only water present

Key calculations:

Per-Section Access

Each TwoFluidSection exposes the oil-water results:

Method Returns Description
getOilWaterFlowRegime() OilWaterFlowRegime Detected regime for this section
getOilWaterResult() OilWaterResult Full result (regime, viscosity, inversion, droplet size, etc.)
isWaterWetting() boolean True if water wets the pipe wall (corrosion risk)
isWaterDropoutRisk() boolean True if water may separate and accumulate
getOilWaterInterfacialTension() double Oil-water IFT (N/m)
setOilWaterInterfacialTension(double) — Override IFT (default: 0.03 N/m)
getOilWaterDetector() OilWaterFlowRegimeDetector Access the detector for tuning
setOilWaterDetector(...) — Set custom detector instance

Tuning the Detector

OilWaterFlowRegimeDetector detector = section.getOilWaterDetector();
detector.setCriticalWeber(1.17);   // Hinze criterion (default 1.17)
detector.setInversionConstant(0.5); // Decarre-Fabre constant (default 0.5)

Holdup Correlations

Minimum Holdup Configuration

The default adaptive minimum is a closure relation that scales with no-slip holdup and tends continuously to zero as liquid input vanishes. It is not a phase-presence threshold. Exact phase presence comes from the conservative gas, oil, and water masses: an absent phase has zero mass, holdup, and velocity.

Configuration Methods

Method Default Description
setUseAdaptiveMinimumOnly(boolean) true Use correlation-based minimum only
setMinimumLiquidHoldup(double) 0.001 Optional absolute floor in fixed-floor mode; zero disables it
setMinimumSlipFactor(double) 2.0 Multiplier for no-slip holdup
setEnforceMinimumSlip(boolean) true Enable/disable minimum constraint

Where the minimum applies. The bound alphaL >= lambdaL * minimumSlipFactor states that the gas outruns the liquid by at least that factor, which is a property of gas-driven transport. It is applied only on level and uphill sections. On a downhill section gravity moves the liquid, the slip ratio legitimately falls, and the bound has no basis; applying it there overwrote the momentum balance with a constant. On a 5 km, 200 mm fixture undulating by +/-30 m it was binding on 39 of 42 downhill sections and on none of the uphill or level ones.

Horizontal Annular Criterion

Method Default Description
setUseEquilibriumLevelAnnularTransition(boolean) true Branch on the equilibrium liquid level instead of the droplet-entrainment criterion

The horizontal branch of the flow map decides annular flow from the equilibrium liquid level, following Taitel and Dukler (1976). Disabling it restores the earlier path, which used the vertical droplet-entrainment criterion U_SG > 3.1 * (sigma * g * drho / rhoG^2)^0.25 ahead of the stratified/slug transition. That threshold is around 0.75 m/s for a 14-inch high-pressure export line, so it classified a horizontal gas pipeline as annular on gas velocity alone and solved a shallow stratified layer with a thin-film closure.

The two paths differ only where the gas velocity clears the droplet threshold while the Kelvin-Helmholtz margin is still below one. On a 73.8 km export line they are identical at 10 MSm3/d; at 4 MSm3/d the equilibrium-level branch reclassifies 272 of 320 sections as stratified-wavy and moves the maximum holdup error from -25.5 to -2.4 per cent.

Friction Model

Method Default Description
setSeparatedFrictionModel(boolean) true Charge each phase its own wall shear where the phases are separated

The friction gradient uses per-phase wall shear, -dP/dx = (tau_wG*S_G + tau_wL*S_L)/A, in stratified flow, and the mixture correlation elsewhere. The mixture form charges the whole perimeter with a holdup-weighted density; on a stratified line that over-predicts the pressure drop by a factor of about 2.3 at 41 per cent holdup, and because it scales as G^2 / rho_mix it also makes extra liquid reduce the gradient, which inverts the terrain response.

The separated form is scoped to stratified flow because its wetted perimeters come from a circular-segment layer at the bottom of the bore. Annular flow, whose film wets the whole perimeter, is not described by that geometry: including it moved the export-line error from +1.4 to +14.7 per cent at 10 MSm3/d and pushed 12 MSm3/d into the pressure floor.

Lean Gas Systems

For lean wet gas (< 1% liquid loading), use adaptive-only mode:

pipe.setUseAdaptiveMinimumOnly(true);  // Default
pipe.setMinimumSlipFactor(2.0);
// Minimum holdup = lambdaL × 2.0 = 0.6% for 0.3% liquid loading

Rich Condensate Systems

For rich gas condensate (> 5% liquid loading), either mode works:

// Option 1: Adaptive (recommended)
pipe.setUseAdaptiveMinimumOnly(true);

// Option 2: Explicit calibrated wetting-film floor
pipe.setUseAdaptiveMinimumOnly(false);
pipe.setMinimumLiquidHoldup(0.01);  // 1% floor

Fixed-floor mode is opt-in and should be used only when a nonzero wetting film is supported by the fluid, wall-wetting, and flow-regime data. Even in this mode, an exactly absent phase remains exactly absent. setMinimumLiquidHoldup(0.0) therefore produces no absolute floor.

Minimum Holdup Correlations

The adaptive minimum uses Beggs-Brill type correlations:

Flow Regime Correlation Exponents
Stratified αL = 0.98 × λL^0.4846 / Fr^0.0868 Segregated flow
Slug/Churn αL = 0.845 × λL^0.5351 / Fr^0.0173 Intermittent flow
Annular Film model + 1.065 × λL^0.5824 / Fr^0.0609 Distributed flow

Where λL = no-slip liquid holdup, Fr = Froude number = v²/(g×D)

Phase Disappearance and Numerical Regularization

TwoFluidPipe keeps closure regularization separate from conserved state:

These values regularize local constitutive equations; they do not declare that a phase is present. The implementation is literature-inspired and does not claim numerical equivalence with OLGA, LedaFlow, or another commercial simulator.

Stratified Flow Holdup

The calculateStratifiedHoldupMomentumBalance() method calculates liquid holdup from momentum balance:

Holdup = f(τwG, τwL, τi, ∂P/∂x, geometry)

Implementation features:

Velocity-Dependent Slip Model

The model captures liquid accumulation at low velocities using Froude number correlation:

// Slip ratio as function of mixture Froude number
double baseSlip = 3.0;
double maxSlip = 25.0;
double exponent = 0.85;
double slip = baseSlip + (maxSlip - baseSlip) * Math.exp(-exponent * Frm);
Parameter Value Physical Meaning
baseSlip 3.0 Minimum slip at high velocity
maxSlip 25.0 Maximum slip at near-zero velocity
exponent 0.85 Velocity sensitivity factor

Terrain Tracking

Terrain Effects Model

The applyTerrainAccumulation() method implements terrain-induced multiphase flow effects:

1. Low Point Liquid Accumulation

Uses Froude number criterion (Fr < 0.5 indicates accumulation):

double Fr_liquid = vL / Math.sqrt(g * diameter * (rhoL - rhoG) / rhoL);
if (Fr_liquid < 0.5) {
    // Calculate accumulated volume based on velocity deficit
}

2. Flowline–Riser Severe-Slugging Stability

Severe slugging is a system instability, not a local pipe-section threshold. After solving a flowline–riser case, call the explicit diagnostic with the index of the first continuously rising section:

SevereSluggingSystemDiagnostic.Result stability =
    pipe.evaluateSevereSluggingSystem(riserBaseSection);

boolean severeSluggingPossible =
    stability.isApplicable() && stability.isSevereSluggingPossible();
double pressureMarginPa = stability.getPressureMarginPa();

The implementation uses the quasi-steady Taitel (1986) condition:

\[P_{top,crit} = \phi\,\rho_L\,g\left(\frac{V_G}{A_r\,\alpha'} - H\right)\]

where $P_{top}$ is absolute riser-outlet pressure [Pa], $\phi$ is average riser liquid holdup [-], $\rho_L$ is average riser liquid density [kg/m³], $V_G$ is upstream compressible gas volume [m³], $A_r$ is riser area [m²], $H$ is vertical riser height [m], and $\alpha’$ is the gas-cap void fraction [-]. The system is classified stable when $P_{top} + \Delta P_{choke} \ge P_{top,crit}$. A static choke pressure drop can be provided to represent one operating point; dynamic choke response is not modelled.

The diagnostic is applicable only to a two-phase, low-rate, stratified flowline followed by a continuously rising, constant-area riser. It assumes isothermal ideal-gas compression and neglects wall and interfacial shear during incipient gas penetration. It deliberately returns a not-applicable status for invalid topology, non-stratified feeders, single-phase states, and unvalidated oil–water–gas cases. It predicts a stability boundary, not slug frequency, slug length, or transient cycle amplitude.

The flowline and riser may have different diameters, and each flowline section contributes its own solved gas volume. Only the rising sections must have constant area. The public SevereSluggingSystemDiagnostic.fromSections(...) factory exposes the extracted descriptor for unit checking, audit, and reuse outside TwoFluidPipe.

The old getSevereSluggingNumberProfile() method is retained as a deprecated serialization/API alias. Its values are the local inclined-section gas-carryover screen now exposed accurately as getInclinedSectionGasCarryoverNumberProfile(); they must not be used as a flowline–riser stability criterion. The associated local flag is available from getInclinedSectionLiquidFallbackPotentialProfile(). getSevereSlugPotentialProfile() is reserved for the explicit system result and is cleared by the next transient step.

Reference: Taitel, Y. (1986), Stability of Severe Slugging, International Journal of Multiphase Flow 12(2), 203–217, doi:10.1016/0301-9322(86)90026-1.

Public severe-slugging benchmark

The diagnostic and transient solver are checked against Tengesdal’s public 2002 air–mineral-oil experiments in a 3-inch, -3-degree flowline and 14.94 m riser. The source data and the assumptions needed to reproduce them are recorded with the tests instead of being treated as an undocumented commercial-simulator comparison.

The diagnostic benchmark uses all 55 operating points in Figure 4-8 and the superficial velocities and uncertainties in Table A-3. Figure symbols were digitized as 26 severe-slug, 14 transition, and 15 stable observations. Transition points are reported separately and are not scored as either binary class. With homogeneous inlet holdup, the published effective upstream volume, 856 kg/m³ liquid density, atmospheric separator pressure, and a 0.89 gas-cap void fraction, the current Taitel screen gives:

Experimental class Predicted severe Predicted stable
Severe slug 22 4
Stable 8 7
Transition (not scored) 6 8

This is 70.7% binary accuracy, 84.6% severe-slug recall, and 46.7% stable recall. It is useful as a conservative screen but is not a high-specificity classifier and must not be described as quantitative dynamic validation.

The slow dynamic benchmark exercises large-facility Test 3 ($v_{SL}=0.50$ m/s and standard $v_{SG}=1.00$ m/s). It is currently a disabled qualification fixture, not a validated model. The coupled route now completes the 16-section 50 × 0.1 s progress probe and a 24-section 100 × 0.05 s refinement with phase and aggregate mass-balance residuals below $10^{-9}$ and without a rejected nonlinear substep. That is a numerical-progress result, not severe-slugging validation: the sticky pressure-correction limiter fires and the interval-average liquid outlet spans -18.55 to 6.88 kg/s, outside the stored 0.375 to 4.03 kg/s comparison range. The legacy route still activates isTransientOutletBackflowClamped(). Both trajectories remain disqualified before amplitude, period, or slug-length agreement is considered. The class remains in source so issues #2909, #2911, and #3298 can re-enable the unchanged public benchmark after the outlet/pressure coupling is resolved; unrelated CI must not treat a known-invalid trajectory as a passing or intentionally failing prediction.

This benchmark previously reported a riser-head-scaled swing of 0.40–0.54 heads and a tracked slug of about 2 m. Those results did not come from the momentum balance. The minimum-slip hold-up bound was written as $\alpha_L \ge \lambda_L \cdot S$, which is a slip statement only in the lean-gas limit: its exact slip ratio is $S\,v_{SG} / (v_{SG} + v_{SL}(1-S))$, which diverges at $v_{SL} = v_{SG}/(S-1)$, and as a hold-up it exceeds one for $\lambda_L > 1/S$. At this facility’s no-slip fraction of 0.33 with $S=2$ it evaluated to 0.67, fed back through the reduced gas area, and saturated at its 0.9 clamp in every section of both the flowline and the riser. The line was held liquid-full by a constant and the reported cycle was that constraint oscillating. The bound is now inverted properly, $\alpha_L \ge X/(1+X)$ with $X = S\,v_{SL}/v_{SG}$, which is below one at every liquid loading; the flowline then solves to a hold-up of 0.334 with the liquid running downhill at 1.51 m/s against a gas velocity of 0.60 m/s.

The riser is the part that is still wrong. The Taylor-bubble film in the slug closure is taken from an annular film balance that cannot close in a riser — the film weight exceeds the gas shear by more than two orders of magnitude — so the iteration stops at its 0.2 thickness clamp and returns a film hold-up of 0.64, which pins the riser slug unit at the 0.9 clamp. A riser that is always liquid-full cannot produce a riser-head pressure swing. Withdrawing the saturated film frees the riser and the swing rises to 1.59 heads at 16 sections, but 24 sections then gives 3.50 heads and reclassifies the riser as bubble flow, so that change is not mesh converged and is not applied.

Severe slugging in this configuration is additionally a deterministically chaotic limit cycle, so instantaneous extremes taken from a single trajectory are not asserted numerically.

Earlier closure-development runs evaluated a four-member ensemble — 16 sections at 0.1 s, the same case with a $10^{-12}$ inlet perturbation, 24 sections at 0.1 s, and 16 sections at 0.2 s — to separate trajectory-robust from trajectory-sensitive quantities. Those runs activated the outlet backflow clamp and therefore remain diagnostic development evidence only:

Quantity Observed across the ensemble How it is asserted
Phase-resolved and total mass closure below $10^{-15}$ below $10^{-10}$
Time-averaged riser-base pressure 176.5–178.3 kPa, spread below 1% mesh, outer-step and perturbation agreement within 8%
Outlet-liquid blowout and fallback present in every realization above 1.25 and below 0.75 of the liquid feed rate
Peak-to-peak riser-base pressure 7.5–10.4 kPa, or 0.060–0.083 riser heads inside 0.02–0.20 riser heads, recorded as a known limitation
Apparent cycle period 13.2–14.4 s each realization above the riser filling time, and the ensemble mean below the experimental 38 ± 2 s
Maximum tracked outlet slug 0 m in every realization asserted to be zero, so a non-zero length forces re-measurement

The green experimental SS trace in Tengesdal Figure 5-6 (printed page 91, physical PDF page 111) gives approximately 98 ± 5 kPa inlet-pressure amplitude and 38 ± 2 s cycle period by direct figure digitization. These are not printed tabular values, and the black Model trace is excluded. The values are independent benchmark targets only; they are not Troll C plant values. No clamped or non-progressing NeqSim trajectory may be compared with them as a valid prediction.

The dynamic reproduction uses the physical 19.81 m flowline plus riser, 0.0762 m diameter, atmospheric outlet, nitrogen as an air surrogate, and a single non-volatile TBP fraction fitted to the reported Crystex density. The source does not give a case-specific temperature or a full oil assay, so 25 °C and the TBP molecular weight are explicit modelling assumptions. The experimental upstream tank/plenum is not represented dynamically; this missing compressible volume is a likely contributor to the short period. Coarse grids are additionally sensitive to whether the flowline–riser boundary lands on a cell face, which is one reason the instantaneous amplitude is not mesh-converged even though the mean pressure is. The steady-state initialization runs with the wall-clock guard disabled and each realization asserts that the guard did not fire, so the reported results do not depend on the speed or load of the executing machine. The stochastic slug tracker uses a fixed benchmark seed; ordinary simulations retain its non-deterministic default. Only the explicit RK4 path is covered; no IMEX severe-slugging validation is claimed.

Run the public checks with:

./mvnw -Dtest=SevereSluggingBenchmarkHarnessTest test
./mvnw -DexcludedTestGroups= -Dtest=SevereSluggingExperimentalBenchmarkTest test

Source: S. Tengesdal, Investigation of Self-Lifting Concept for Severe Slugging Elimination in Deep-Water Pipeline/Riser Systems (2002), BSEE Technical Assessment Program report.

3. Uphill Liquid Fallback

Uses Turner droplet model for critical gas velocity:

double vG_critical = 3.0 * Math.pow(sigma * g * (rhoL - rhoG) / (rhoG * rhoG), 0.25);
if (vG < vG_critical) {
    // Liquid fallback occurs
}

4. Downhill Drainage

double drainageRate = Math.sqrt(2 * g * dz * holdup);

Multi-Layer Thermal Model

New Classes

  1. RadialThermalLayer - Represents a single thermal layer with material properties
  2. MultilayerThermalCalculator - Calculates U-value and transient heat transfer

Supported Layer Materials

Material k [W/(m·K)] ρ [kg/m³] Cp [J/(kg·K)]
Carbon Steel 50.0 7850 480
FBE Coating 0.3 1400 1000
PU Foam 0.035 80 1500
Syntactic Foam 0.15 650 1100
Aerogel 0.015 150 1000
Concrete 1.4 2400 880

Usage Example

TwoFluidPipe pipe = new TwoFluidPipe("subsea-export", inletStream);
pipe.setLength(20000.0); // 20 km
pipe.setDiameter(0.254); // 10 inch
pipe.setWallThickness(0.015);
pipe.setSurfaceTemperature(4.0, "C"); // Cold seabed

// Configure with 50mm PU foam + 40mm concrete
pipe.configureSubseaThermalModel(0.050, 0.040,
    RadialThermalLayer.MaterialType.PU_FOAM);

// Explicit shutdown assumption; the documented default is also 50 W/(m2 K)
pipe.setStagnantInnerHeatTransferCoefficient(50.0);

// Set hydrate formation temperature
pipe.setHydrateFormationTemperature(20.0, "C");

// Calculate cooldown time
double cooldownHours = pipe.calculateHydrateCooldownTime();
System.out.printf("Cooldown to hydrate: %.1f hours%n", cooldownHours);

// Run simulation
pipe.run();

// Get thermal summary
System.out.println(pipe.getThermalSummary());

Thermal Calculations

For validation, start from run(), close both boundaries, disable Joule–Thomson effects for an adiabatic invariant, and check that a uniform state remains uniform. For cooldown, report absolute pressure, composition and mixing rule, temperatures, heat-transfer coefficients, wall properties, mesh, time step, and units; verify that every cell approaches ambient monotonically without undershoot. For deterministic compact regressions, require the thermal report to meet an absolute residual of 1e-5 J or a relative residual of 1e-10, then demonstrate mesh and time-step refinement at nearby conditions. This implementation does not claim OLGA or LedaFlow equivalence.

Model Capabilities Summary

Category Feature Method/Correlation
Conservation Equations    
Gas mass Full continuity equation Phase-resolved flash transfer
Oil mass Full continuity equation Equilibrium mass split / donor inventory
Water mass Full continuity equation Equilibrium mass split / donor inventory
Gas momentum 1D momentum balance Wall and interfacial shear
Liquid momentum 1D momentum balance Wall and interfacial shear
Mixture energy Full energy balance Optional J-T effect
Closure Models    
Stratified holdup Momentum balance Taitel-Dukler geometry
Annular holdup Film model Ishii-Mishima entrainment
Slug holdup Empirical correlation Dukler correlation
Interfacial friction Flow-regime specific Multiple correlations
Oil-Water Models    
Oil-water flow regime OilWaterFlowRegimeDetector Trallero/Brauner/Angeli classification
Phase inversion Decarre-Fabre (1997) Viscosity/density-ratio model
Emulsion viscosity Brinkman correlation Continuous/dispersed mixture
Water wetting Per-section detection Corrosion risk indicator
Water dropout Velocity/holdup criterion Accumulation risk flag
Terrain Effects    
Low point accumulation Froude criterion Fr < 0.5 triggers accumulation
Riser-base liquid fallback Local gas-carryover screen Indicates possible local fallback only
Flowline-riser stability Taitel (1986) quasi-steady criterion Explicit topology-aware system diagnostic
Uphill fallback Turner model Critical gas velocity check
Thermal Model    
Multi-layer heat transfer Series resistance RadialThermalLayer class
Cooldown calculation Lumped capacitance MultilayerThermalCalculator
Hydrate/wax risk Temperature tracking Section-by-section monitoring
Numerical Methods    
Time stepping CFL-based RK4 (default), IMEX, adaptive dt
Spatial discretization Finite volume AUSM+ flux splitting, MUSCL reconstruction
Mesh Uniform or non-uniform generateRefinedMesh() or setSectionLengths()

Steady-State and Dynamic Simulation Workflow

TwoFluidPipe supports both steady-state initialization and transient simulation. In normal use, call run() first to build a physically consistent pressure, holdup, temperature, and flow-regime profile. Then change a boundary condition and advance time with runTransient(dt, id) until the profiles stop changing or match a new steady-state reference.

Steady-State Simulation

For a fixed inlet stream and a calculated outlet pressure, configure the pipe geometry and call run():

SystemInterface fluid = new neqsim.thermo.system.SystemSrkEos(293.15, 70.0);
fluid.addComponent("methane", 0.90);
fluid.addComponent("ethane", 0.06);
fluid.addComponent("propane", 0.04);
fluid.setMixingRule("classic");

Stream inlet = new Stream("inlet", fluid);
inlet.setFlowRate(4.0, "kg/sec");
inlet.setTemperature(20.0, "C");
inlet.setPressure(70.0, "bara");
inlet.run();

TwoFluidPipe pipe = new TwoFluidPipe("export line", inlet);
pipe.setLength(1000.0);
pipe.setDiameter(0.20);
pipe.setRoughness(1.0e-5);
pipe.setNumberOfSections(20);
pipe.run();

If the downstream pressure is known, set a constant outlet pressure before run(). The steady-state solver calculates the pressure-gradient shape from flow, friction, gravity, and holdup, then aligns the absolute profile to the specified pressure boundary:

pipe.setOutletPressure(55.0, "bara");
pipe.run();
double outletPressure = pipe.getPressureProfile()[pipe.getPressureProfile().length - 1] / 1.0e5;

Dynamic Simulation After a Boundary Change

Dynamic simulations should start from a steady state. After changing a boundary condition, run the transient solver long enough for the new stationary solution to be reached:

pipe.run();
pipe.setOutletPressure(52.0, "bara");

UUID id = UUID.randomUUID();
double elapsedTime = 0.0;
while (elapsedTime < 60.0) {
  pipe.runTransient(2.0, id);
  elapsedTime += 2.0;
}

For regression tests, compare the transient profile after the boundary change with a second TwoFluidPipe solved directly at the new boundary condition. A practical pressure-profile metric is root-mean-square pressure difference across all sections; for compact tests a limit of 1-2 bar is a useful sanity check. Treat acoustic pressure settling and material-inventory settling as separate checks. Pressure waves can settle quickly, while liquid, oil, and water holdup profiles may require one or more residence times to approach a stationary distribution. Do not impose the pressure-test horizon on holdup convergence or force agreement by reconstructing conservative phase masses from the stationary closure.

Transient holdup is reconstructed from the phase masses advanced by the finite-volume equations. There is no separate post-step projection toward the steady-state holdup correlation. Such a projection changes phase inventory without a boundary flux or mass-transfer source and makes the error scale with pipe length. The earlier unreferenced 4 s relaxation time has therefore been removed; steady-state closures remain part of initialization and the local closure/source terms. Auxiliary terrain and slug trackers may maintain primitive diagnostics, but they do not rebuild the finite-volume phase masses. Conservative source or flux coupling for those trackers remains future model-development work.

At the steady-to-transient handoff, run() converts the final pressure, phase-holdup, density, and velocity profiles into conservative phase mass, momentum, and energy exactly once. This conversion defines the initial condition and advances no simulation time. For three-phase flow, the oil and water momenta retain the independent phase velocities from the steady slip closure rather than being collapsed to the bulk-liquid velocity. After the transient solve starts, the conservative phase masses own cell inventory. A stream-connected inlet may update boundary composition and velocity for its inlet flux, but it must not replace the first finite-volume cell’s oil or water mass. This prevents an unchanged, near-zero-time handoff from producing an inventory or holdup jump that scales with pipe volume.

At an exact oil-only or water-only liquid endpoint, the active phase velocity is synchronized with the final bulk-liquid velocity before the conservative state is built. Inspect the phase-resolved steady fluxes when qualifying a new transient case:

double[] gasFlow = pipe.getGasMassFlowProfile();
double[] oilFlow = pipe.getOilMassFlowProfile();
double[] waterFlow = pipe.getWaterMassFlowProfile();
int outlet = gasFlow.length - 1;
double steadyOutletFlow = gasFlow[outlet] + oilFlow[outlet] + waterFlow[outlet];

For a no-transfer steady case, steadyOutletFlow should reproduce the inlet mass flow within the chosen numerical tolerance. This is a kinematic handoff check only. It does not prove that the phase-momentum sources and fixed-pressure outlet form a transient fixed point, and it does not qualify the current liquid-rich or severe-slugging trajectories.

For each phase $k$ and for the total domain, validate the discrete balance

\[M_k(t + \Delta t) - M_k(t) = \int_t^{t+\Delta t} \left(\dot m_{k,in} - \dot m_{k,out} + S_k\right)\,dt, \qquad M = \sum_i \left(m'_{g,i} + m'_{o,i} + m'_{w,i}\right)\Delta x_i .\]

Use getTotalMassInventory() to read $M$ in kg directly from the conservative gas, oil, and water masses. After each runTransient(...), getLastMassBalanceReport() provides initial and final inventory, integrated inlet and outlet fluxes, integrated sources, signed residual in kg, and relative residual for GAS, OIL, WATER, LIQUID, and TOTAL:

pipe.runTransient(0.1, UUID.randomUUID());
TwoFluidMassBalanceReport balance = pipe.getLastMassBalanceReport();
double totalResidualKg = balance.getResidualKg(TwoFluidMassBalanceReport.Phase.TOTAL);
double totalRelativeResidual = balance.getRelativeResidual(TwoFluidMassBalanceReport.Phase.TOTAL);
boolean closes = balance.isWithinTolerance(TwoFluidMassBalanceReport.Phase.TOTAL, 1.0e-7, 1.0e-10);

After runTransient(...), getOutletStream() publishes the accepted interval-average total outlet mass flux:

double outletMassFlow = balance.getOutletMassKg(TwoFluidMassBalanceReport.Phase.TOTAL)
    / balance.getElapsedTimeSeconds();

This preserves finite transport delay and inventory release when the pipe is coupled directly to downstream process equipment. Steady-state run() retains the inlet-balanced outlet-flow convention. Use the report itself for signed, phase-resolved flux integrals and conservation evidence.

The boundary and source integrals use the accepted internal substeps and the same Euler, Runge-Kutta, or IMEX stage weights as the conservative update. For deterministic regression cases, an absolute tolerance of $10^{-7}$ kg or a relative tolerance of $10^{-10}$ is appropriate; choose a larger engineering tolerance for long simulations after demonstrating time-step and mesh sensitivity. Positivity limiting or any future non-conservative correction appears explicitly as a non-zero residual rather than being hidden.

Flash-driven phase transfer is added to gas, oil, and water with $\Gamma_G+\Gamma_O+\Gamma_W=0$, so it cancels from the total source without losing liquid identity. Oil-water segregation is represented by separate phase momentum and face fluxes; it does not use a local oil-to-water mass relaxation. The IMEX pressure correction likewise changes phase momenta, not phase masses. Closed boundaries therefore have zero integrated boundary flux, while open boundaries close against the actual phase-resolved inlet and outlet fluxes.

A steady-state solution remains a useful long-time comparison, but agreement must result from the transient balances and closure forces rather than overwriting the conserved state.

Boundary Conditions

Use setInletBoundaryCondition(...) and setOutletBoundaryCondition(...) for explicit boundary types. Convenience methods such as closeOutlet() and openOutlet(...) update the type and value together.

Boundary condition Typical side Required value How to set it Use case
STREAM_CONNECTED Inlet Inlet stream flow, temperature, pressure, and composition Default inlet; pipe.openInlet() Pipe connected to upstream process equipment
CONSTANT_FLOW Inlet Mass flow pipe.setInletBoundaryCondition(BoundaryCondition.CONSTANT_FLOW); pipe.setInletMassFlow(4.0, "kg/sec"); Production-rate step or controlled inlet flow
CONSTANT_PRESSURE Inlet or outlet Pressure pipe.setInletPressure(70.0, "bara") or pipe.setOutletPressure(55.0, "bara") Known upstream pressure or downstream back-pressure
CLOSED Inlet or outlet None pipe.closeInlet() or pipe.closeOutlet() Shut-in, valve closure, blocked-in pipe
CHARACTERISTIC Inlet or outlet External pressure/flow state pipe.setOutletBoundaryCondition(BoundaryCondition.CHARACTERISTIC) Reduced wave reflection in fast transients

Example with explicit boundary settings:

pipe.setInletBoundaryCondition(TwoFluidPipe.BoundaryCondition.CONSTANT_FLOW);
pipe.setInletMassFlow(4.0, "kg/sec");
pipe.setOutletBoundaryCondition(TwoFluidPipe.BoundaryCondition.CONSTANT_PRESSURE);
pipe.setOutletPressure(55.0, "bara");
pipe.run();

pipe.openOutlet(52.0, "bara");
pipe.runTransient(2.0, UUID.randomUUID());

Avoid over-specifying the same side. For example, a STREAM_CONNECTED inlet already gets flow, temperature, pressure, and composition from the inlet stream; switch to CONSTANT_FLOW only when the flow should be independent of the stream flow rate.

Compressible upstream-volume boundary

A real well, flowline, or laboratory loop normally has compressible inventory upstream of the modeled pipe. A fixed stream rate removes that storage and can shorten a severe-slug cycle. UpstreamCompressibleVolume supplies a phase-resolved dynamic inlet pressure without adding pipeline closures to the boundary model.

Initialize the pipe to steady state, create the volume from its inlet section, and specify the external phase source rates:

pipe.run();
UpstreamCompressibleVolume sourceVolume =
    pipe.initializeUpstreamCompressibleVolume(25.0);
sourceVolume.setSourceMassFlowRates(gasSourceKgS, oilSourceKgS, waterSourceKgS);

pipe.runTransient(0.1, UUID.randomUUID());
double sourcePressurePa = sourceVolume.getPressurePa();
double volumeClosure = sourceVolume.getMaximumRelativeVolumeResidual();

After every accepted internal substep, the pipe passes its integration-weighted gas, oil, and water inlet masses to the volume. The volume applies

M_k(new) = M_k(old) + source_k dt - pipe_inlet_k
sum_k M_k / rho_k(p) = fixed volume,  with d rho_k / d p = 1 / c_k^2.

Signed pipe transfer is supported, so a phase returning from the pipe adds upstream inventory. Phase depletion and non-converged pressure closure throw rather than clamping mass. Connecting the volume selects a constant-pressure inlet whose value is updated from the volume. When a volume is attached, pipe.getInletPressure() reports that current boundary pressure after the accepted inventory transfer; the axial pressure profile remains the accepted internal pipe state. The ordinary stream remains the source of thermodynamic composition for the pipe boundary. Component mixing, heat transfer, and a momentum-resolved vessel/nozzle model are outside this lumped boundary’s current scope.

Multi-branch pressure-node network

TwoFluidPipeNetwork connects named TwoFluidPipe branches through fixed-pressure reservoirs and phase-resolved compressible storage nodes. It supports directed splits, merges, manifolds, and injection branches:

TwoFluidPipeNetwork network = new TwoFluidPipeNetwork("subsea network");
network.addFixedPressureNode("well", 120.0e5);
network.addCompressibleNode("manifold", manifoldVolume);
network.addFixedPressureNode("separator", 55.0e5);
network.addPipe("flowline A", "well", "manifold", flowlineA);
network.addPipe("riser", "manifold", "separator", riser);

network.runTransient(0.1, UUID.randomUUID());
double nodePressure = network.getNodePressurePa("manifold");
double closure = network.getLastBalanceReport()
    .getRelativeResidual(TwoFluidMassBalanceReport.Phase.TOTAL);

Every branch sees the node pressures at the beginning of the network step. After all branches advance, their accepted gas, oil, and water boundary transfers update every compressible node simultaneously. Branch iteration order therefore does not become an implicit pressure solve. The whole-network report includes pipe and node inventories, fixed-boundary transfers, phase sources, and gas/oil/water/liquid/total residuals.

This first network implementation requires finite compressible storage at internal junctions. Fixed-pressure nodes are external reservoirs or sinks. Algebraic zero-volume junctions, a global Newton pressure-flow solve, component and enthalpy mixing, reverse-flow boundary composition, and branch subcycling remain outside the validated scope.

Choosing Time Step, Sections, and Pipeline Length

runTransient(dt, id) takes a macro time step in seconds. Internally, the solver sub-steps to satisfy the CFL condition, so dt is the requested reporting/control interval, not necessarily the single numerical step.

Guidelines:

Choice Practical guidance
Number of sections Start with 10-20 sections for simple horizontal pipes, 30-100 for long pipelines, and refine around risers, low points, or sharp elevation changes.
Section length Keep dx = length / sections small enough to resolve terrain and holdup changes. A low point or riser should span several sections, not one cell.
Time step Start with 0.5-2 s for compact regression models and 5-60 s for long slow-flow pipelines when adaptive time stepping is enabled. Reduce it for valve closures, slug fronts, or fast pressure waves.
Adaptive stepping setEnableAdaptiveTimestepping(true) is recommended for difficult multiphase transients. It recomputes the stable internal step as velocities and holdups change.
CFL number setCflNumber(0.3-0.8) is typical. Lower values provide more stability margin and temporal resolution; higher values are faster.
Settling time Run until pressure and holdup profiles stop changing or match a new stationary reference. The required physical time scales with pipeline length, flow velocity, compressibility, and liquid inventory.

For a pressure-boundary step, a compact 300 m regression pipe may settle in a few seconds. A real long subsea line can require minutes to hours of simulated time, especially when liquid inventory, terrain accumulation, or thermal transients dominate. Always distinguish wall-clock runtime from physical simulated time.

Spatial Discretization

Uniform Mesh (default)

setNumberOfSections(N) creates N equal-length cells: $dx = L / N$.

Non-Uniform Mesh

Two approaches for variable cell sizes along the pipe:

Automatic refinement — generateRefinedMesh(baseSections, refinementFactor) analyses the elevation profile and creates shorter cells where the elevation gradient is steepest (risers, S-bends) and longer cells where the pipe is flat (flowlines):

\[\text{density}_i = 1 + (\text{factor} - 1) \cdot \frac{|\nabla z|_i}{\max |\nabla z|}\]

Section lengths are inversely proportional to density, then normalized to sum to $L$. The refinementFactor (clamped to 1.5–10) controls the coarsest/finest cell ratio.

Manual — setSectionLengths(double[]) sets explicit per-section lengths (must sum to total pipe length, minimum 2 sections).

All finite-volume calculations use per-section lengths:

Component Non-uniform treatment
AUSM+ flux assembly $-\frac{1}{dx_i}(F_{i+1/2} - F_{i-1/2})$
Pressure gradient Non-uniform central difference: $dx_c = \frac{1}{2} dx_{i-1} + dx_i + \frac{1}{2} dx_{i+1}$
CFL timestep $\Delta t = \min_i \left( \text{CFL} \cdot dx_i / c_i \right)$
Temperature updates Per-section exponential decay and advection
Pressure reconstruction Forward/backward march with per-section $dx$

Time Integration

Methods

Select the time integration method via setTimeIntegrationMethod(TimeIntegrator.Method):

Method CFL constraint Description
RK4 (default) Acoustic ($c + v$) Classical 4th-order Runge-Kutta. Stable for all geometries.
SSP_RK3 Acoustic Strong Stability Preserving RK3
RK2 Acoustic Heun’s method (2nd order)
EULER Acoustic Forward Euler (1st order)
IMEX_PRESSURE_CORRECTION Convective only Semi-implicit momentum pressure correction with conservative explicit phase-mass transport; ~10x larger dt. Not recommended for vertical risers. 
pipe.setTimeIntegrationMethod(TimeIntegrator.Method.RK4);       // default
pipe.setTimeIntegrationMethod(TimeIntegrator.Method.IMEX_PRESSURE_CORRECTION); // semi-implicit
TimeIntegrator.Method current = pipe.getTimeIntegrationMethod(); // query

Coupled pressure-momentum transient correction

The historical transient path advances conservative phase mass and momentum and then reconstructs pressure from the steady friction-and-gravity gradient. That sequential reconstruction cannot provide compressibility feedback to a momentum instability. It is retained as the default for compatibility.

The opt-in coupled correction solves a finite-volume pressure equation from the cell-volume constraint and phase compressibilities. The pressure correction changes phase mass fluxes conservatively and corrects gas, oil, and water momenta with the same face pressure gradients. It therefore advances pressure and momentum in one accepted substep and does not call the steady pressure reconstruction afterward.

pipe.setEnableInterfacialPressure(true);
pipe.setImplicitInterfacialPressureCoupling(true);
pipe.setEnableCoupledPressureMomentum(true);
pipe.setAllowOutletPhaseBackflow(true);
pipe.setCoupledPressureMomentumMaximumIterations(24); // default
pipe.setCoupledPressureMomentumRelativeVolumeTolerance(1.0e-7); // default

Use the four settings together for liquid-rich transients whose pressure outlet permits phase fallback. The interfacial-pressure term makes the phase-momentum system hyperbolic, its implicit treatment removes the small void-wave CFL limit, the coupled correction supplies compressibility feedback, and the signed outlet prevents a reversed liquid phase from being silently pinned at zero. Signed fallback extrapolates the interior phase state and remains opt-in because a one-way export boundary may instead require a check valve or a supplied external inflow composition. The mode reports isCoupledPressureMomentumConverged(), getCoupledPressureMomentumVolumeResidual(), and getCoupledPressureMomentumIterations(). The nonlinear gate is available through get/setCoupledPressureMomentumMaximumIterations() and get/setCoupledPressureMomentumRelativeVolumeTolerance(). The default iteration budget is 24; the former budget of 12 stopped the Tengesdal progress case at about $6\times10^{-7}$ relative cell-volume residual, above the $10^{-7}$ default gate.

After a transient window, also read the sticky diagnostics isTransientCoupledPressureMomentumFailureDetected(), isTransientCoupledPressureMomentumCorrectionLimited(), and getTransientCoupledPressureMomentumRejectedSubsteps(). They are cleared by the next steady run(). isCoupledPressureMomentumPressureCorrectionLimited() reports only the latest correction. Adaptive execution retries a non-converged correction at a smaller internal step. If the requested interval is still incomplete at the minimum factor, runTransient throws with accepted and requested elapsed time, residual, tolerance, iteration count/cap, and limiter state. It never returns a zero- or partial-progress coupled interval silently. A limiter event is diagnostic rather than an automatic rejection because the nonlinear residual can converge while a bounded correction is active; it must still be reported in qualification evidence.

The option remains off by default while long-horizon severe-slugging, mesh, timestep, and independent OLGA/public benchmark qualification is completed. Do not use short startup peaks as limit-cycle evidence; report period, the P10-P90 band, and completed cycle count over a settled window.

Standing benchmark acceptance metrics

Use TwoFluidBenchmarkMetrics to compare a rate sweep, a section profile, mesh realizations, and a settled transient with the same definitions in every benchmark:

double exponent = TwoFluidBenchmarkMetrics.fitRateExponent(rates, pressureDrops);
double terrainLocalization =
    TwoFluidBenchmarkMetrics.maximumToMedianRatio(liquidHoldupProfile);
double meshSpread =
    TwoFluidBenchmarkMetrics.relativeMeshSpread(coarseResult, refinedResult);
TwoFluidBenchmarkMetrics.LimitCycleMetrics cycle =
    TwoFluidBenchmarkMetrics.analyzeLimitCycle(times, pressure, settledWindowStart);

The limit-cycle result reports period, P10, median, P90, P10-P90 band, sample window, and completed median-upcrossing cycle count. A large startup peak followed by a flat trace reports zero completed cycles, so it cannot pass as a sustained oscillation.

The slow public Tengesdal benchmark now uses these settled-window definitions and requires at least two completed liquid-rate cycles in every mesh/timestep realization. It retains the published experimental source and its phase-mass, mean-pressure, mesh, amplitude, and deterministic-repeat checks.

A fresh OLGA 2025.1 execution for the same geometry reached normal stop and reported a 34.9234 kPa pressure amplitude and 21.7100 s liquid-trough period (21.7364 s from pressure), compared with the approximately 98 ± 5 kPa and 38 ± 2 s Figure 5-6 digitizations. That one-point run is 64.36% low in amplitude and 42.87% short in period, so it is a parseable comparator but not a qualified benchmark reproduction. These values are comparison evidence, not embedded commercial correlations. A NeqSim/OLGA comparison must record exported time series and evaluate both with the same settled window and metrics; it must not compare one startup peak or tune a closure to one trajectory.

Adaptive Timestepping

Adaptive timestepping provides robustness for challenging geometries. Enable via setEnableAdaptiveTimestepping(true).

Algorithm per macro-step:

  1. CFL recompute from current velocities (not fixed at initialization)
  2. Pre-check: reject if NaN or negative mass detected; rollback state, halve dtFactor
  3. Post-check: reject if pressure exceeds ceiling or velocities exceed 500 m/s
  4. Recovery: after each stable step, dtFactor grows by x1.02 back toward 1.0
  5. Floor: dtFactor cannot go below 0.001 to prevent stalling

Steady-State Solver Tuning

The initial steady-state solve iterates between the transient solver and thermodynamic flashes until convergence. Four parameters control this:

Parameter Setter Default Description
Under-relaxation setSteadyStateUnderRelaxation(double) 0.5 Update damping factor (0–1); lower = more damping, more stable
Flash interval setSteadyStateFlashInterval(int) 3 Re-flash thermodynamics every N iterations; higher = faster but less accurate
Max iterations setSteadyStateMaxIterations(int) 0 = mesh-scaled 0 uses max(100, 20 x sections); the sweep moves information about one section per iteration, so a fixed budget silently truncates long, finely-discretised lines
Max wall-clock time setSteadyStateMaxWallClockTime(double) 300 s Timeout for the SS solver; prevents runaway iterations 
pipe.setSteadyStateUnderRelaxation(0.3);   // More conservative damping
pipe.setSteadyStateFlashInterval(5);       // Flash every 5th iteration
pipe.setSteadyStateMaxWallClockTime(60.0); // Allow 60 seconds

Always check the steady-state outcome

run() does not throw when the steady state fails to settle, so the outcome has to be read back. Three independent flags describe it, and a profile is only trustworthy when the first is true:

Query Meaning when true
isSteadyStateConverged() The sweep met the tolerance and the profile is a solution
isSteadyStateWallClockLimited() The wall-clock guard stopped the sweep early
isSteadyStatePressureFloorLimited() One or more sections rest on the internal 1 bara pressure floor 
pipe.run();
if (!pipe.isSteadyStateConverged()) {
  if (pipe.isSteadyStatePressureFloorLimited()) {
    throw new IllegalStateException(
        "The line cannot deliver this rate at this inlet pressure");
  }
  throw new IllegalStateException("Steady state did not converge");
}

Why the pressure floor matters. The marching solver clamps every section at 1 bara so it stays numerically alive when a line has no deliverability. A profile resting on that clamp is a fixed point of the clamp, not of the momentum balance: the per-section change falls below tolerance and the sweep would otherwise report success on a line that cannot physically deliver the requested rate. isSteadyStatePressureFloorLimited() makes that case visible, and isSteadyStateConverged() is withheld. PipeBeggsAndBrills throws Outlet pressure is negative on the same condition, so both codes agree that such a case has no solution.

Always check the transient outcome too

runTransient(dt, id) exposes invalid outcomes that can still satisfy a discrete balance. Coupled non-convergence now throws if the requested interval cannot be completed; sticky diagnostics must still be read after successful calls.

Query Meaning when true
isTransientOutletBackflowClamped() A phase reversed at the outlet and its outflow is pinned at zero
isTransientCoupledPressureMomentumFailureDetected() At least one coupled nonlinear correction was rejected
isTransientCoupledPressureMomentumCorrectionLimited() At least one nonlinear pressure correction reached its limiter
getTransientCoupledPressureMomentumRejectedSubsteps() > 0 Coupled substeps were retried at a smaller adaptive factor 
pipe.runTransient(dt, null);
if (pipe.isTransientOutletBackflowClamped()) {
  throw new IllegalStateException("Transient liquid inventory is running away");
}
if (pipe.isTransientCoupledPressureMomentumFailureDetected()
    || pipe.isTransientCoupledPressureMomentumCorrectionLimited()) {
  throw new IllegalStateException("Coupled transient requires diagnostic review");
}

Any trajectory with one of these diagnostic flags is invalid for engineering comparison or model qualification until the event has been explained and accepted, even if its mass-balance residual, pressure amplitude, period, or other headline metric appears acceptable. The diagnostics are sticky for the transient run and must be checked after the full evaluated window.

Why outlet backflow matters. The outlet is transmissive: it can carry mass out but not in, so a reversed phase velocity is clamped to zero. That is correct as a boundary condition and is also a one-way trap. The phase momentum equations of the classical two-fluid system are ill-posed at high liquid fraction and can develop sustained backflow, after which the outflow of that phase pins at exactly 0 kg/s while the inlet keeps feeding the line and the inventory grows without bound. The finite-volume balance still closes to machine precision throughout, so the balance alone cannot detect the defect. Gas-dominated lines do not show it. Interfacial pressure alone is not a qualified remedy. The combined implicit-interfacial, coupled-pressure/momentum, and signed-outlet route now advances the short Tengesdal progress probe, but its sticky limiter and outlet-rate mismatch still disqualify the trajectory.

Direct electrical heating (DEH)

A uniform electrical heat input can be added to the energy equation, in steady state and in transient runs, and it works with wall heat transfer switched off:

pipe.setLength(73845.0);
pipe.setDirectElectricalHeatingPower(10.0e6);        // W, spread over the pipe length
// or, equivalently:
pipe.setDirectElectricalHeatingPowerPerMeter(135.4); // W/m

The power set here is what reaches the fluid, so cable and coating losses must already be deducted. The same convention is used by PipeBeggsAndBrills.setDirectElectricalHeatingPower(double), so the two models can be compared like for like.

In steady state the segment solution decays toward the wall-loss/DEH balance temperature

\[T_\infty = T_{surf} + \frac{q}{U \pi D}\]

rather than toward the surface temperature. This is exact for a uniform source, so the profile cannot overshoot the balance temperature — unlike explicit per-increment stepping, which does.

Validation Status

Evidence levels

Passing software tests establish numerical regressions, API behavior, and conservation. They do not by themselves establish agreement with experiment. Current external evidence is:

Commercial transient multiphase simulators are not used as a reference. Their licence terms generally prohibit publishing benchmark comparisons and prohibit using the software to develop the science, technology or product content of similar software, so no NeqSim closure is tuned to such a tool and no measured deviation against one is recorded here.

The public severe-slugging benchmark deliberately retains failed/limited metrics in its assertions and documentation. In particular, the riser-base pressure amplitude, the cycle period and the slug-length result all prevent a claim of fully quantitative severe-slugging validation.

Steady-state behaviour on a long gas-condensate export line

Measured on a 73.8 km subsea gas-condensate export line (ID 0.355 m, U = 3 W/m2K, seabed 4 C, 200 bara inlet) as an internal consistency check on the solver. The rate exponent in $\Delta P \sim \dot m^{\,n}$ rises from about 2.1 at low rate to about 3.1 at high rate, so the density feedback along the line is reproduced rather than merely the level at one rate. Adding 10 MW of direct electrical heating raises the arrival temperature 17.4 K and the pressure drop 15.0 per cent. PipeBeggsAndBrills sits far above TwoFluidPipe on the same cases because its two-phase friction multiplier is an extrapolation at this liquid loading. Results are grid-converged, at default settings.

Remaining limitations:

Implemented regression tests

Integration Tests (TwoFluidPipeIntegrationTest)

Validation Tests (TwoFluidPipeValidationTest)

Beggs-Brill Correlation Comparison:

Pipeline Scenario Validation:

Terrain-Induced Slugging Patterns:

The exact suite size changes as the model evolves. Use the Maven/JUnit result for the tested commit rather than a hard-coded historical count.

References

  1. Bendiksen, K.H., Maines, D., Moe, R., & Nuland, S. (1991). “The Dynamic Two-Fluid Model OLGA: Theory and Application.” SPE Production Engineering, 6(02), 171-180.

  2. Taitel, Y., & Dukler, A.E. (1976). “A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow.” AIChE Journal, 22(1), 47-55.

  3. Pots, B.F.M., Bromilow, I.G., & Konijn, M.J.W.F. (1987). “Severe Slug Flow in Offshore Flowline/Riser Systems.” SPE Production Engineering, 2(04), 319-324.

  4. Turner, R.G., Hubbard, M.G., & Dukler, A.E. (1969). “Analysis and Prediction of Minimum Flow Rate for the Continuous Removal of Liquids from Gas Wells.” Journal of Petroleum Technology, 21(11), 1475-1482.

  5. Bai, Y., & Bai, Q. (2010). “Subsea Pipelines and Risers.” Elsevier. Chapter on Thermal Design.

  6. Beggs, H.D. & Brill, J.P. (1973). “A Study of Two-Phase Flow in Inclined Pipes.” Journal of Petroleum Technology, SPE-4007-PA.