TwoFluidPipe Model Documentation

Overview

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

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

Conservation Equations

Mass Conservation

Separate mass conservation equations for gas and liquid phases:

Equation Mathematical Form Description
Gas mass ∂(αG ρG)/∂t + ∂(αG ρG vG)/∂x = ΓG Gas phase continuity with mass transfer
Liquid mass ∂(αL ρL)/∂t + ∂(αL ρL vL)/∂x = -ΓG Liquid phase continuity with mass transfer
Mass transfer ΓG Flash-based calculation Evaporation/condensation with optional kinetic limits

Where:

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

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 model enforces a minimum liquid holdup to prevent unrealistically low values in gas-dominant systems. By default, an adaptive minimum is used that scales with the no-slip holdup, making it suitable for both lean gas and rich condensate systems.

Configuration Methods

Method Default Description
setUseAdaptiveMinimumOnly(boolean) true Use correlation-based minimum only
setMinimumLiquidHoldup(double) 0.001 Absolute floor (when adaptive-only = false)
setMinimumSlipFactor(double) 2.0 Multiplier for no-slip holdup
setEnforceMinimumSlip(boolean) true Enable/disable minimum constraint

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: Fixed floor (OLGA-style)
pipe.setUseAdaptiveMinimumOnly(false);
pipe.setMinimumLiquidHoldup(0.01);  // 1% 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)

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. Riser Base Severe Slugging

Detects severe slugging potential using Pots criterion:

double pi_ss = (inletPressure - outletPressure) / (rhoL * g * riserHeight);
if (pi_ss > 1.0) {
    // Severe slugging potential flagged
}

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);

// 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

Model Capabilities Summary

Category Feature Method/Correlation
Conservation Equations    
Gas mass Full continuity equation Flash-based mass transfer
Liquid mass Full continuity equation Flash-based mass transfer
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 slugging Pots criterion πSS > 1.0 indicates severe slugging
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. For multiphase cases, also compare liquid, oil, and water holdup profiles so the dynamic calculation has settled hydraulically, not only acoustically.

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.

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 are more stable; higher values are faster but less conservative.
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; ~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

Adaptive Timestepping

OLGA-style 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. Three 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 wall-clock time setSteadyStateMaxWallClockTime(double) 30 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

Validation Status

Implemented Tests

Integration Tests (TwoFluidPipeIntegrationTest)

Validation Tests (TwoFluidPipeValidationTest)

Beggs-Brill Correlation Comparison:

Pipeline Scenario Validation:

Terrain-Induced Slugging Patterns:

Test Coverage Summary

Test Category Tests Status
Integration Tests 24 ✅ All passing
Validation Tests 13 ✅ All passing
Total 37 ✅ All passing

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.