TwoFluidPipe Model: Detailed Review and External Comparison

Overview

The TwoFluidPipe class in NeqSim implements a transient two-fluid model for 1D multiphase pipeline flow. This document reviews the implementation and records traceable comparisons. API names containing OLGA are retained for compatibility; they identify literature-inspired NeqSim closure selections and do not claim numerical equivalence with that or any commercial simulator.

Table of Contents

  1. Mathematical Foundation
  2. Conservation Equations
  3. Closure Relations
  4. Slug Flow Modeling
  5. Lagrangian Slug Tracking
  6. Terrain-Induced Slugging
  7. Usage Examples
  8. Validation Results
  9. References

Mathematical Foundation

Two-Fluid Model Fundamentals

The two-fluid model treats gas and liquid as interpenetrating continua, each with their own velocity, density, and momentum. The model solves conservation equations for:

Governing Equations in Conservative Form

The 1D two-fluid equations in conservative form:

\[\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial x} = \mathbf{S}(\mathbf{U})\]

Where the state vector $\mathbf{U}$, flux vector $\mathbf{F}$, and source terms $\mathbf{S}$ are:

\[\mathbf{U} = \begin{pmatrix} \alpha_G \rho_G A \\ \alpha_L \rho_L A \\ \alpha_G \rho_G v_G A \\ \alpha_L \rho_L v_L A \\ E_{mix} A \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \alpha_G \rho_G v_G A \\ \alpha_L \rho_L v_L A \\ \alpha_G \rho_G v_G^2 A + \alpha_G P A \\ \alpha_L \rho_L v_L^2 A + \alpha_L P A \\ (E_{mix} + P) v_m A \end{pmatrix}\] \[\mathbf{S} = \begin{pmatrix} \Gamma_G \\ \Gamma_L \\ -\tau_{wG} S_G - \tau_i S_i - \alpha_G \rho_G g \sin\theta \cdot A \\ -\tau_{wL} S_L + \tau_i S_i - \alpha_L \rho_L g \sin\theta \cdot A \\ -q_{wall} \pi D + \dot{m} \Delta h \end{pmatrix}\]

Notation

Symbol Description Unit
$\alpha$ Phase holdup (volume fraction) -
$\rho$ Density kg/m³
$v$ Velocity m/s
$P$ Pressure Pa
$A$ Pipe cross-sectional area m²
$\tau_w$ Wall shear stress Pa
$\tau_i$ Interfacial shear stress Pa
$S$ Wetted/interfacial perimeter m
$g$ Gravitational acceleration m/s²
$\theta$ Pipe inclination angle rad
$\Gamma$ Mass transfer rate kg/(m·s)

Conservation Equations

Mass Conservation

Gas phase: \(\frac{\partial (\alpha_G \rho_G)}{\partial t} + \frac{\partial (\alpha_G \rho_G v_G)}{\partial x} = \Gamma_G\)

Liquid phase: \(\frac{\partial (\alpha_L \rho_L)}{\partial t} + \frac{\partial (\alpha_L \rho_L v_L)}{\partial x} = \Gamma_L\)

Where $\Gamma_G = -\Gamma_L$ (mass transfer between phases) and the constraint $\alpha_G + \alpha_L = 1$ must be satisfied.

Momentum Conservation

Gas phase: \(\frac{\partial (\alpha_G \rho_G v_G)}{\partial t} + \frac{\partial (\alpha_G \rho_G v_G^2)}{\partial x} = -\alpha_G \frac{\partial P}{\partial x} - \frac{\tau_{wG} S_G}{A} - \frac{\tau_i S_i}{A} - \alpha_G \rho_G g \sin\theta\)

Liquid phase: \(\frac{\partial (\alpha_L \rho_L v_L)}{\partial t} + \frac{\partial (\alpha_L \rho_L v_L^2)}{\partial x} = -\alpha_L \frac{\partial P}{\partial x} - \frac{\tau_{wL} S_L}{A} + \frac{\tau_i S_i}{A} - \alpha_L \rho_L g \sin\theta\)

Energy Conservation (Mixture)

\[\frac{\partial E_{mix}}{\partial t} + \frac{\partial}{\partial x}\left[(E_{mix} + P) v_m\right] = -\frac{q_{wall} \pi D}{A} + \dot{Q}_{source}\]

Where mixture energy: \(E_{mix} = \alpha_G \rho_G \left(e_G + \frac{v_G^2}{2}\right) + \alpha_L \rho_L \left(e_L + \frac{v_L^2}{2}\right)\)

Closure Relations

Wall Friction

Wall shear stress using Fanning friction factor:

\[\tau_{wk} = \frac{1}{2} f_k \rho_k v_k |v_k|\]

Friction factor (Haaland correlation): \(\frac{1}{\sqrt{f}} = -1.8 \log_{10}\left[\left(\frac{\epsilon/D}{3.7}\right)^{1.11} + \frac{6.9}{Re}\right]\)

Hydraulic diameter for stratified flow: \(D_{hG} = \frac{4 A_G}{S_G + S_i}, \quad D_{hL} = \frac{4 A_L}{S_L + S_i}\)

Interfacial Friction

Stratified smooth flow (Taitel-Dukler): \(\tau_i = \frac{1}{2} f_i \rho_G (v_G - v_L)|v_G - v_L|\)

Where $f_i = f_G$ (gas friction factor).

Stratified wavy flow (Andritsos-Hanratty): \(f_i = f_G \left[1 + 15 \sqrt{\frac{h_L}{D}} \left(\frac{v_G - v_{G,crit}}{v_{G,crit}}\right)^{0.5}\right]\)

Critical gas velocity: \(v_{G,crit} = 5 \sqrt{\frac{\rho_L - \rho_G}{\rho_G}} \sqrt{\frac{g h_L}{\cos\theta}}\)

Annular flow (Wallis): \(f_i = 0.005 \left[1 + 300 \frac{\delta}{D}\right]\)

Where $\delta$ is the liquid film thickness.

Stratified Flow Geometry

For a circular pipe with liquid height $h_L$:

Central angle: \(\phi = 2 \cos^{-1}\left(1 - \frac{2h_L}{D}\right)\)

Cross-sectional areas: \(A_L = \frac{D^2}{8}(\phi - \sin\phi), \quad A_G = A - A_L\)

Wetted perimeters: \(S_L = \frac{D\phi}{2}, \quad S_G = \frac{D(2\pi - \phi)}{2}, \quad S_i = D\sin\left(\frac{\phi}{2}\right)\)

Slug Flow Modeling

Slug Unit Model

A slug unit consists of:

  1. Slug body: High liquid holdup region ($H_{LS} \approx 0.7-1.0$)
  2. Taylor bubble: Elongated gas pocket
  3. Film region: Thin liquid film beneath the Taylor bubble
    ┌─────────────────────────────────────────────────────────┐
    │                                                         │
    │  ←── Taylor Bubble ──→  ←───── Slug Body ─────→        │
    │         (Gas)                  (Liquid)                 │
    │  ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓  ████████████████████████████    │
    │  ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓  ████████████████████████████    │
    │══════════════════════════════════════════════════════   │
    │    Film Region         │                                │
    │                        │                                │
    └────────────────────────┴────────────────────────────────┘
         ← L_bubble →          ←────── L_slug ──────→
         
    ◄──────────────── L_unit = L_slug + L_bubble ─────────────►

Taylor Bubble Velocity

Bendiksen (1984) correlation: \(v_{TB} = C_0 \cdot v_m + v_d\)

Distribution coefficient $C_0$: \(C_0 = \begin{cases} 1.2 & \text{if } Fr_m > 3.5 \\ 1.05 + 0.15\sin\theta & \text{if } Fr_m \leq 3.5 \end{cases}\)

Drift velocity $v_d$:

For horizontal flow (Zukoski 1966): \(v_{dH} = 0.54 \sqrt{\frac{g D (\rho_L - \rho_G)}{\rho_L}}\)

For vertical flow (Dumitrescu 1943): \(v_{dV} = 0.35 \sqrt{\frac{g D (\rho_L - \rho_G)}{\rho_L}}\)

Interpolation for inclined pipes: \(v_d = v_{dH} \cos\theta + v_{dV} \sin\theta\)

Slug Body Holdup

Gregory et al. (1978): \(H_{LS} = \frac{1}{1 + \left(\frac{v_m}{8.66}\right)^{1.39}}\)

This correlation accounts for gas entrainment in the slug body at high mixture velocities.

Slug Frequency

Zabaras (2000) correlation: \(f_s = \frac{0.0226 \cdot \lambda_L^{1.2} \cdot Fr_m^{2.0}}{D} \cdot (1 + \sin|\theta|)\)

Where:

Equilibrium Slug Length

Barnea-Taitel (1993): \(\frac{L_s}{D} = 25 + 10 \cdot \min(Fr_m, 2.0)\)

With inclination correction: \(L_s = D \cdot \frac{L_s}{D} \cdot (1 + 0.3\sin\theta) \quad \text{for } \theta > 0\)

Lagrangian Slug Tracking

Overview

The Lagrangian slug tracking model tracks individual slugs as discrete entities propagating through the pipeline. Each slug has:

Slug Dynamics Equations

Front velocity: \(v_{front} = C_0 \cdot v_m + v_d\)

Tail velocity (from mass balance): \(v_{tail} = v_{front} \cdot \phi_{shedding}\)

Where the shedding factor depends on slug length relative to equilibrium: \(\phi_{shedding} = \begin{cases} 0.95 & \text{if } L_s < 0.9 L_{eq} \text{ (growing)} \\ 0.98 & \text{if } 0.9 L_{eq} \leq L_s \leq 1.2 L_{eq} \text{ (stable)} \\ 1.0 + 0.1(L_s/L_{eq} - 1.2) & \text{if } L_s > 1.2 L_{eq} \text{ (decaying)} \end{cases}\)

Slug length evolution: \(\frac{dL_s}{dt} = v_{front} - v_{tail}\)

Position update: \(x_{front}^{n+1} = x_{front}^n + v_{front} \cdot \Delta t\) \(x_{tail}^{n+1} = x_{tail}^n + v_{tail} \cdot \Delta t\)

Mass Exchange

Pickup rate at front (liquid scooped from film): \(\dot{m}_{pickup} = \rho_L \cdot A \cdot H_{film} \cdot (v_{front} - v_{film})\)

Shedding rate at tail (liquid shed to film): \(\dot{m}_{shed} = \rho_L \cdot A \cdot (H_{LS} - H_{film}) \cdot (v_{tail} - v_{slug})\)

Net mass rate: \(\frac{dm_s}{dt} = \dot{m}_{pickup} - \dot{m}_{shed}\)

Wake Effects

Following slugs experience acceleration in the wake of preceding slugs:

\[v_{following} = v_{base} \cdot C_{wake}\]

Wake coefficient: \(C_{wake} = C_{max} - (C_{max} - 1) \cdot \frac{d}{L_{wake}}\)

Where:

Slug Merging

When the front of a following slug catches the tail of a preceding slug:

\[\text{if } x_{front,following} \geq x_{tail,preceding} - \epsilon_{merge}\]

The slugs merge:

Slug Tracking Modes

// Detailed NeqSim Lagrangian tracking (default)
pipe.setSlugTrackingMode(TwoFluidPipe.SlugTrackingMode.LAGRANGIAN);

// Simplified slug unit model
pipe.setSlugTrackingMode(TwoFluidPipe.SlugTrackingMode.SIMPLIFIED);

// Disable slug tracking
pipe.setSlugTrackingMode(TwoFluidPipe.SlugTrackingMode.DISABLED);

Terrain-Induced Slugging

Liquid Accumulation Model

At terrain low points, liquid accumulates when gas velocity is insufficient to sweep the liquid forward.

Gas Froude number: \(Fr_G = \frac{v_{SG}}{\sqrt{g D \frac{\rho_L - \rho_G}{\rho_G}}}\)

Critical Froude number: $Fr_{crit} \approx 1.5$

Below the critical Froude number, liquid accumulates: \(\text{Accumulation factor} = 1 + A \cdot \left(1 - \frac{Fr_G}{Fr_{crit}}\right)^{1.5}\)

Where $A \approx 10$ is an amplitude factor.

Slug Release Criterion

A terrain-induced slug is released when:

  1. Holdup exceeds threshold: $\alpha_L > \alpha_{crit}$ (typically 0.6)
  2. Sufficient volume accumulated: $V_{acc} > V_{min}$
  3. Gas velocity increases (pressure buildup behind liquid plug)

Explicit Flowline–Riser Stability Diagnostic

The former local “Bøe criterion” description did not match an implemented system model and has been removed. Severe slugging depends on the compressible upstream gas volume, the riser geometry, and the downstream pressure response. NeqSim now exposes the Taitel (1986) quasi-steady system screen:

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

Use evaluateSevereSluggingSystem(riserBaseSection) after solving the pipe. The result reports applicability, stable/unstable status, critical top pressure, pressure margin, stability ratio, and gas-expansion head. A flowline-to-riser diameter change is permitted, while the continuously rising riser itself must have constant area. The diagnostic rejects non-stratified, single-phase, three-phase, variable-area-riser, and invalid flowline–riser topologies. This screen does not predict dynamic slug cycles and is not a claim of equivalence to a commercial simulator.

Usage Examples

Basic Two-Phase Pipe Flow

import neqsim.process.equipment.pipeline.TwoFluidPipe;
import neqsim.process.equipment.pipeline.twophasepipe.SevereSluggingSystemDiagnostic;
import neqsim.process.equipment.stream.Stream;
import neqsim.thermo.system.SystemSrkEos;

// Create fluid system
SystemSrkEos fluid = new SystemSrkEos(298.15, 50.0);
fluid.addComponent("methane", 0.8);
fluid.addComponent("n-heptane", 0.2);
fluid.setMixingRule("classic");

// Create inlet stream
Stream inlet = new Stream("Inlet", fluid);
inlet.setFlowRate(100.0, "kg/hr");
inlet.setTemperature(25.0, "C");
inlet.setPressure(50.0, "bara");
inlet.run();

// Create TwoFluidPipe
TwoFluidPipe pipe = new TwoFluidPipe("Pipeline", inlet);
pipe.setLength(5000.0);           // 5 km
pipe.setDiameter(0.2);            // 200 mm (8 inch)
pipe.setNumberOfSections(50);     // 100 m per section
pipe.setInclination(0.0);         // Horizontal
pipe.setOutletPressure(45.0, "bara");

// Run steady-state
pipe.run();

// Print results
System.out.println("Pressure drop: " + 
    (pipe.getInletPressure() - pipe.getOutletPressure()) + " bar");
System.out.println("Outlet temperature: " + 
    pipe.getOutletStream().getTemperature("C") + " °C");

Transient Simulation with Slug Tracking

// Configure for transient simulation
pipe.setOLGAModelType(TwoFluidPipe.OLGAModelType.FULL);
pipe.setSlugTrackingMode(TwoFluidPipe.SlugTrackingMode.LAGRANGIAN);

// Configure Lagrangian tracking
pipe.configureLagrangianSlugTracking(
    true,   // enableInletGeneration
    true,   // enableTerrainGeneration
    true    // enableWakeEffects
);

// Advanced tracker configuration
LagrangianSlugTracker tracker = pipe.getLagrangianSlugTracker();
tracker.setMinSlugLengthDiameters(12.0);    // Minimum stable slug
tracker.setMaxSlugLengthDiameters(300.0);   // Maximum slug length
tracker.setInitialSlugLengthDiameters(20.0); // Initial slug length
tracker.setWakeLengthDiameters(30.0);        // Wake region
tracker.setMaxWakeAcceleration(1.3);         // Wake acceleration

// Run transient for 1 hour
pipe.runTransient(3600.0);

// Get slug statistics
System.out.println(pipe.getSlugStatisticsSummary());

// Access detailed slug data
System.out.println("\nActive slugs:");
for (LagrangianSlugTracker.SlugBubbleUnit slug : tracker.getSlugs()) {
    System.out.printf("  Slug #%d: pos=%.1fm, L=%.1fm, v=%.2fm/s, H=%.2f%n",
        slug.id, slug.frontPosition, slug.slugLength, 
        slug.frontVelocity, slug.slugHoldup);
}

// Outlet statistics
System.out.println("\nOutlet slug statistics:");
System.out.printf("  Slugs exited: %d%n", tracker.getTotalSlugsExited());
System.out.printf("  Max volume: %.4f m³%n", tracker.getMaxSlugVolumeAtOutlet());
System.out.printf("  Outlet frequency: %.4f Hz%n", tracker.getOutletSlugFrequency());

The tracker-level liquid-mass audit is borrowed - returned to Eulerian - exited at outlet - active slug mass. Use getMassConservationError() for the residual and getTotalMassExitedAtOutlet() for the cumulative interval-accounted outlet mass. This diagnostic audit is distinct from the conservative phase-resolved TwoFluidMassBalanceReport produced by TwoFluidPipe.

Terrain Profile with Slugging

// Create terrain profile (undulating pipeline)
double[] distances = {0, 1000, 2000, 3000, 4000, 5000};  // m
double[] elevations = {0, -50, -100, -50, -150, 0};      // m (relative)

// Set terrain profile
pipe.setTerrainProfile(distances, elevations);
pipe.setEnableTerrainTracking(true);
pipe.setEnableTerrainSlugClosures(true);

// Configure terrain parameters
pipe.setTerrainSlugCriticalHoldup(0.6);
pipe.setLiquidFallbackCoefficient(0.3);

// Run transient simulation
pipe.runTransient(7200.0);  // 2 hours

// Evaluate the explicit system after selecting the first rising section
int riserBaseSection = 40;
SevereSluggingSystemDiagnostic.Result stability =
    pipe.evaluateSevereSluggingSystem(riserBaseSection);
boolean severeSluggingPossible = stability.isSevereSluggingPossible();
double pressureMarginPa = stability.getPressureMarginPa();
// Pass the status and margin to the application's logger or monitoring system.

// Get holdup profile
double[] positions = pipe.getPositionProfile();
double[] holdups = pipe.getLiquidHoldupProfile();

System.out.println("\nHoldup along pipe:");
for (int i = 0; i < positions.length; i++) {
    System.out.printf("  x=%.0fm: αL=%.3f%n", positions[i], holdups[i]);
}

Heat Transfer with Insulation

// Enable heat transfer
pipe.enableHeatTransfer(true);

// Configure multi-layer insulation
pipe.setInsulationType(TwoFluidPipe.InsulationType.SUBSEA_INSULATED);

// Or manual layer configuration
MultilayerThermalCalculator thermal = pipe.getThermalCalculator();
thermal.clearLayers();
thermal.addLayer(MultilayerThermalCalculator.LayerMaterial.CARBON_STEEL, 0.020);  // 20mm wall
thermal.addLayer(MultilayerThermalCalculator.LayerMaterial.FBE_COATING, 0.0004);  // 0.4mm FBE
thermal.addLayer(MultilayerThermalCalculator.LayerMaterial.PU_FOAM, 0.060);       // 60mm PU foam
thermal.addLayer(MultilayerThermalCalculator.LayerMaterial.CONCRETE, 0.040);      // 40mm concrete

// Set ambient conditions
pipe.setSurfaceTemperature(4.0, "C");  // Seabed temperature
pipe.setHeatTransferCoefficient(50.0); // W/(m²·K) outer HTC

// Run with heat transfer
pipe.run();

// Get temperature profile
double[] temps = pipe.getTemperatureProfile();
System.out.printf("Temperature: %.1f°C (inlet) → %.1f°C (outlet)%n",
    temps[0] - 273.15, temps[temps.length-1] - 273.15);

// Check hydrate risk
System.out.printf("Hydrate formation temperature: %.1f°C%n", 
    thermal.getHydrateFormationTemperature() - 273.15);
System.out.printf("Cooldown time to hydrate: %.1f hours%n", 
    thermal.getCooldownTimeToHydrate());

Python/Jupyter Integration

# Using neqsim-python with direct Java access
from jpype import JClass

# Import NeqSim classes
SystemSrkEos = JClass('neqsim.thermo.system.SystemSrkEos')
Stream = JClass('neqsim.process.equipment.stream.Stream')
TwoFluidPipe = JClass('neqsim.process.equipment.pipeline.TwoFluidPipe')

# Create fluid
fluid = SystemSrkEos(298.15, 50.0)
fluid.addComponent("methane", 0.85)
fluid.addComponent("ethane", 0.10)
fluid.addComponent("propane", 0.05)
fluid.setMixingRule("classic")

# Create stream and pipe
inlet = Stream("Inlet", fluid)
inlet.setFlowRate(5000.0, "kg/hr")
inlet.setTemperature(40.0, "C")
inlet.setPressure(80.0, "bara")
inlet.run()

pipe = TwoFluidPipe("Subsea Pipeline", inlet)
pipe.setLength(20000.0)  # 20 km
pipe.setDiameter(0.254)  # 10 inch
pipe.setNumberOfSections(100)
pipe.setOutletPressure(50.0, "bara")

# Enable Lagrangian slug tracking
SlugTrackingMode = JClass('neqsim.process.equipment.pipeline.TwoFluidPipe$SlugTrackingMode')
pipe.setSlugTrackingMode(SlugTrackingMode.LAGRANGIAN)

# Run transient
pipe.runTransient(3600.0)

# Get results for plotting
import numpy as np
positions = np.array(pipe.getPositionProfile())
pressures = np.array(pipe.getPressureProfile()) / 1e5  # Convert to bar
holdups = np.array(pipe.getLiquidHoldupProfile())

# Plot results
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

ax1.plot(positions/1000, pressures)
ax1.set_xlabel('Distance (km)')
ax1.set_ylabel('Pressure (bar)')
ax1.set_title('Pressure Profile')
ax1.grid(True)

ax2.plot(positions/1000, holdups)
ax2.set_xlabel('Distance (km)')
ax2.set_ylabel('Liquid Holdup (-)')
ax2.set_title('Liquid Holdup Profile')
ax2.grid(True)

plt.tight_layout()
plt.show()

# Print slug statistics
print(pipe.getSlugStatisticsSummary())

Validation Results

From the comparison tests (TwoFluidVsBeggsBrillComparisonTest):

Test Case Beggs-Brill ΔP TwoFluid ΔP Difference
Horizontal gas-dominant 0.145 bar 0.138 bar 5.3%
150mm diameter 0.263 bar 0.249 bar 5.5%
200mm diameter 0.060 bar 0.057 bar 5.2%
300mm diameter 0.008 bar 0.007 bar 5.1%

Terrain slug detection successfully identifies:

These software scenarios are regression checks, not experimental validation.

Public severe-slugging benchmark

Tengesdal’s 2002 public air–mineral-oil experiments provide the current external check. On the 55-point -3-degree velocity map, the Taitel system diagnostic scores 22/26 severe observations and 7/15 stable observations correctly; the 14 transition observations remain a separate category. Dynamic large-facility Test 3 is a deterministically chaotic limit cycle, so it is evaluated as a four-member ensemble. The time-averaged riser-base pressure (171–176 kPa) and the blowout/fallback regime signature are reproducible, while the peak-to-peak pressure spans 42–300 kPa against 98 ± 5 kPa digitized and the apparent period spans 14–35 s against 38 ± 2 s. The current slug tracker also underpredicts the experimental severe-slug length. These limits preclude a claim of quantitative severe-slugging validation.

Source: Tengesdal (2002), BSEE Technical Assessment Program.

Lagrangian slug-tracker regression evidence

LagrangianSlugTrackerTest checks the implementation against its selected correlations and internal invariants:

  1. Correlation consistency: slug frequency, velocity, and holdup calculations
  2. Bookkeeping: tracked, returned, and exited liquid mass
  3. Merging behavior: coalescence when a following slug catches a preceding slug
  4. Wake behavior: the configured acceleration response

These are unit/regression checks, not independent experimental validation of predicted slug frequency or length. The public Tengesdal comparison above is the relevant external evidence.

Numerical Methods

AUSM+ Scheme

The TwoFluidPipe uses the AUSM+ (Advection Upstream Splitting Method Plus) numerical scheme for flux computation:

\[\mathbf{F}_{i+1/2} = \dot{m}_{i+1/2}^+ \boldsymbol{\phi}_L + \dot{m}_{i+1/2}^- \boldsymbol{\phi}_R + P_{i+1/2} \mathbf{n}\]

Where:

Mach number splitting (van Leer): \(\mathcal{M}^{\pm} = \pm\frac{1}{4}(M \pm 1)^2 \quad \text{for } |M| \leq 1\)

MUSCL Reconstruction

For second-order spatial accuracy, Monotonic Upstream-centered Scheme for Conservation Laws:

\[\mathbf{U}_{i+1/2}^L = \mathbf{U}_i + \frac{1}{4}\left[(1-\kappa)\tilde{\Delta}_{i-1/2} + (1+\kappa)\tilde{\Delta}_{i+1/2}\right]\] \[\mathbf{U}_{i+1/2}^R = \mathbf{U}_{i+1} - \frac{1}{4}\left[(1+\kappa)\tilde{\Delta}_{i+1/2} + (1-\kappa)\tilde{\Delta}_{i+3/2}\right]\]

Where $\kappa = 1/3$ gives third-order upwind bias, and $\tilde{\Delta}$ are slope-limited differences.

van Leer slope limiter: \(\psi(r) = \frac{r + |r|}{1 + |r|}\)

Time Integration

For transient simulations, explicit time-stepping with CFL-based time step:

\[\Delta t = \text{CFL} \cdot \min_i \left(\frac{\Delta x_i}{|v_i| + c_i}\right)\]

Where $c$ is the mixture sound speed. Typical CFL = 0.5-0.8 for stability.

Recommendations for Reproducible NeqSim Results

  1. Use FULL model type for best accuracy: 
    pipe.setOLGAModelType(TwoFluidPipe.OLGAModelType.FULL);
  2. Enable Lagrangian slug tracking for detailed slug analysis: 
    pipe.setSlugTrackingMode(TwoFluidPipe.SlugTrackingMode.LAGRANGIAN);
pipe.configureLagrangianSlugTracking(true, true, true);
  3. Enable terrain tracking for undulating pipelines: 
    pipe.setEnableTerrainTracking(true);
pipe.setEnableTerrainSlugClosures(true);
  4. Configure minimum holdup based on system: 
    // For lean gas systems
pipe.setMinimumLiquidHoldup(0.01); // 1%
// For rich gas/condensate
pipe.setMinimumLiquidHoldup(0.02); // 2%
  5. Enable heat transfer for long pipelines: 
    pipe.enableHeatTransfer(true);
pipe.setHeatTransferCoefficient(10.0); // W/(m²·K)
pipe.setSurfaceTemperature(4.0, "C"); // Seabed temperature

References

Two-Fluid Model Theory

  1. Bendiksen, K.H. et al. (1991) “The Dynamic Two-Fluid Model OLGA: Theory and Application” SPE Production Engineering - Foundational OLGA paper

  2. Taitel, Y. and 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. Barnea, D. (1987) “A Unified Model for Predicting Flow-Pattern Transitions for the Whole Range of Pipe Inclinations” Int. J. Multiphase Flow 13(1):1-12

  4. Ishii, M. and Hibiki, T. (2011) “Thermo-Fluid Dynamics of Two-Phase Flow” Springer - Comprehensive two-fluid model reference

Slug Flow Correlations

  1. Bendiksen, K.H. (1984) “An Experimental Investigation of the Motion of Long Bubbles in Inclined Tubes” Int. J. Multiphase Flow 10(4):467-483 - Taylor bubble velocity

  2. Gregory, G.A., Nicholson, M.K., and Aziz, K. (1978) “Correlation of the Liquid Volume Fraction in the Slug for Horizontal Gas-Liquid Slug Flow” Int. J. Multiphase Flow 4(1):33-39 - Slug holdup

  3. Zabaras, G.J. (2000) “Prediction of Slug Frequency for Gas/Liquid Flows” SPE Journal 5(3):252-258 - Slug frequency correlation

  4. Barnea, D. and Taitel, Y. (1993) “A Model for Slug Length Distribution in Gas-Liquid Slug Flow” Int. J. Multiphase Flow 19(5):829-838 - Equilibrium slug length

Friction and Closure

  1. Andritsos, N. and Hanratty, T.J. (1987) “Influence of interfacial waves in stratified gas-liquid flows” AIChE Journal 33(3):444-454 - Wavy flow interfacial friction

  2. Wallis, G.B. (1969) “One-Dimensional Two-Phase Flow” McGraw-Hill - Classic two-phase flow reference

  3. Beggs, H.D. and Brill, J.P. (1973) “A Study of Two-Phase Flow in Inclined Pipes” Journal of Petroleum Technology 25(5):607-617

Terrain and Severe Slugging

  1. Bøe, A. (1981) “Severe Slugging Characteristics” Selected Topics in Two-Phase Flow, NTH, Trondheim

  2. Taitel, Y. (1986) “Stability of Severe Slugging” Int. J. Multiphase Flow 12(2):203-217

Numerical Methods

  1. Liou, M.S. (1996) “A Sequel to AUSM: AUSM+” Journal of Computational Physics 129(2):364-382 - AUSM+ scheme

  2. van Leer, B. (1979) “Towards the Ultimate Conservative Difference Scheme V: A Second-Order Sequel to Godunov’s Method” Journal of Computational Physics 32(1):101-136 - MUSCL reconstruction

Appendix A: Nomenclature

Symbol Description Unit
$A$ Pipe cross-sectional area m²
$\alpha$ Phase volume fraction (holdup) -
$c$ Sound speed m/s
$C_0$ Distribution coefficient -
$D$ Pipe diameter m
$E$ Total energy per unit volume J/m³
$f$ Friction factor -
$Fr$ Froude number $= v/\sqrt{gD}$ -
$g$ Gravitational acceleration m/s²
$h$ Enthalpy J/kg
$H$ Liquid height in stratified flow m
$L$ Length m
$\dot{m}$ Mass flow rate kg/s
$M$ Mach number -
$P$ Pressure Pa
$q$ Heat flux W/m²
$Re$ Reynolds number -
$\rho$ Density kg/m³
$S$ Wetted/interfacial perimeter m
$t$ Time s
$\tau$ Shear stress Pa
$\theta$ Pipe inclination rad
$v$ Velocity m/s
$v_d$ Drift velocity m/s
$v_m$ Mixture velocity m/s
$v_s$ Superficial velocity m/s
$v_{TB}$ Taylor bubble velocity m/s
$x$ Axial position m
$\Gamma$ Mass transfer rate kg/(m·s)
$\lambda$ No-slip holdup (input fraction) -
$\mu$ Dynamic viscosity Pa·s
$\phi$ Central angle (stratified geometry) rad

Subscripts

Subscript Meaning
G, g Gas phase
L, l Liquid phase
O, o Oil phase
W, w Water phase
i Interface
m Mixture
S Superficial
TB Taylor bubble
w Wall

Appendix B: Flow Regime Map

                    Superficial Gas Velocity (m/s)
                0.1     1       10      100
    0.01  ┌─────────────────────────────────┐
          │         STRATIFIED              │
          │    SMOOTH    │    WAVY          │
    0.1   ├─────────────────────────────────┤
          │              │                  │
Super-    │   SLUG       │     ANNULAR      │
ficial    │              │                  │
Liquid    ├─────────────────────────────────┤
Velocity  │              │                  │
(m/s)     │              │     MIST         │
    1.0   │   ELONGATED  │                  │
          │   BUBBLE     │                  │
          ├─────────────────────────────────┤
    10    │         DISPERSED BUBBLE        │
          └─────────────────────────────────┘

Transition criteria implemented:

Document generated for NeqSim TwoFluidPipe model. Last updated with comprehensive mathematical documentation and Lagrangian slug tracking implementation.