neqsim-flow-assurance

Flow Assurance Analysis with NeqSim

Consolidated guide for all flow assurance threats — hydrate, wax, asphaltene, corrosion, hydraulics, water/liquid hammer screening, slugging, and thermal management — with NeqSim code patterns.

When to Use This Skill

• Hydrate formation temperature/pressure prediction
• Hydrate inhibitor dosing (MEG, methanol, ethanol)
• Wax appearance temperature (WAT) and wax deposition risk
• Asphaltene stability screening (de Boer, CII)
• CO2 and H2S corrosion rate estimation
• Elemental sulfur (S8) deposition risk at pressure/temperature letdown (compressor inlets, control/letdown valves, dry-gas seals, filters)
• Pipeline pressure drop and temperature profile
• DNV-RP-F109 vertical and lateral on-bottom stability screening
• Water hammer/liquid hammer screening for fast valve closure, pump trip, or check-valve slam
• Multiphase flow pattern prediction (slug, annular, stratified)
• Thermal insulation sizing for subsea pipelines
• Arrival temperature and cooldown calculations

Applicable Standards

For fast acoustic transients, also load neqsim-water-hammer. Use WaterHammerStudy or MCP runWaterHammer with STID route geometry, tagreader event windows, and valve/pump event schedules; use this flow-assurance skill for the broader operating-envelope and mitigation context.

For on-bottom stability, load neqsim-subsea-and-wells and use the typed DnvRpF109OnBottomStabilityKernel. It provides a transparent absolute-static screen and checks externally calculated response displacements. It does not contain generalized design tables, produce dynamic response, qualify environmental or soil models, or claim DNV conformity. A pass still requires independent review.

1. Hydrate Analysis

EOS Selection for Hydrate Calculations

CPA is required for water–hydrocarbon hydrate modeling. When dissolved salts / electrolytes are present, use SystemElectrolyteCPAstatoil so the salt thermodynamic (hydrate-suppression) effect is captured — plain SystemSrkCPAstatoil ignores ion activity and underestimates subcooling margin.

Hydrate Formation Temperature

// CPA EOS required for accurate water-hydrocarbon modeling
SystemInterface fluid = new SystemSrkCPAstatoil(273.15 + 10, 100.0);
fluid.addComponent("methane", 0.80);
fluid.addComponent("ethane", 0.05);
fluid.addComponent("propane", 0.03);
fluid.addComponent("CO2", 0.02);
fluid.addComponent("water", 0.10);
fluid.setMixingRule(10);  // CPA mixing rule
fluid.setMultiPhaseCheck(true);
fluid.setHydrateCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.hydrateFormationTemperature();
double hydrateT_C = fluid.getTemperature() - 273.15;

Hydrate Equilibrium Curve (Multiple Pressures)

// Calculate hydrate T at several pressures for the full curve
double[] pressures = {20, 40, 60, 80, 100, 120, 150, 200};
double[] hydrateTemps = new double[pressures.length];

for (int i = 0; i < pressures.length; i++) {
    SystemInterface testFluid = fluid.clone();
    testFluid.setPressure(pressures[i]);
    ThermodynamicOperations testOps = new ThermodynamicOperations(testFluid);
    testOps.hydrateFormationTemperature();
    hydrateTemps[i] = testFluid.getTemperature() - 273.15;
}

Hydrate Inhibitor Dosing (MEG)

// Add MEG to suppress hydrate formation temperature
SystemInterface inhibitedFluid = new SystemSrkCPAstatoil(273.15 + 4, 100.0);
inhibitedFluid.addComponent("methane", 0.80);
inhibitedFluid.addComponent("water", 0.15);
inhibitedFluid.addComponent("MEG", 0.05);  // 25 wt% MEG in water phase
inhibitedFluid.setMixingRule(10);
inhibitedFluid.setMultiPhaseCheck(true);
inhibitedFluid.setHydrateCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(inhibitedFluid);
ops.hydrateFormationTemperature();
double inhibitedHydrateT = inhibitedFluid.getTemperature() - 273.15;
// Compare with uninhibited to get subcooling margin

Salt-Inhibited Hydrate (Formation Water / Brine)

// Dissolved salts depress the hydrate temperature — use the electrolyte CPA model
SystemInterface brineGas = new SystemElectrolyteCPAstatoil(273.15 + 4, 100.0);
brineGas.addComponent("methane", 0.80);
brineGas.addComponent("water", 0.18);
brineGas.addComponent("Na+", 0.01);   // dissociated NaCl
brineGas.addComponent("Cl-", 0.01);
brineGas.setMixingRule(10);
brineGas.setMultiPhaseCheck(true);
brineGas.setHydrateCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(brineGas);
ops.hydrateFormationTemperature();
double brineHydrateT = brineGas.getTemperature() - 273.15;
// SystemSrkCPAstatoil would miss the salt suppression — always use electrolyte CPA with ions

MEG Concentration Sweep

// Find required MEG concentration for target subcooling
double[] megWtPct = {0, 10, 20, 30, 40, 50};
for (double wt : megWtPct) {
    // Create fluid with appropriate MEG/water ratio
    double waterFrac = 0.20 * (1.0 - wt / 100.0);
    double megFrac = 0.20 * (wt / 100.0);
    SystemInterface testFluid = new SystemSrkCPAstatoil(273.15, 100.0);
    testFluid.addComponent("methane", 0.80);
    testFluid.addComponent("water", waterFrac);
    testFluid.addComponent("MEG", megFrac);
    testFluid.setMixingRule(10);
    testFluid.setHydrateCheck(true);
    ThermodynamicOperations testOps = new ThermodynamicOperations(testFluid);
    testOps.hydrateFormationTemperature();
    // Record hydrate T vs MEG concentration
}

2. Wax Analysis

Wax Appearance Temperature (WAT)

// Oil system with C7+ fractions for wax prediction
SystemInterface oil = new SystemSrkEos(273.15 + 60, 50.0);
oil.addComponent("methane", 0.30);
oil.addComponent("ethane", 0.10);
oil.addTBPfraction("C7", 0.10, 92.0 / 1000, 0.727);
oil.addTBPfraction("C10", 0.15, 134.0 / 1000, 0.78);
oil.addTBPfraction("C15", 0.15, 206.0 / 1000, 0.83);
oil.addPlusFraction("C20", 0.20, 350.0 / 1000, 0.88);
oil.getCharacterization().setWaxModel(true);
oil.getCharacterization().characterisePlusFraction();
oil.setMixingRule("classic");

ThermodynamicOperations ops = new ThermodynamicOperations(oil);
ops.calcWAT();
double wat_C = oil.getTemperature() - 273.15;

Wax Fraction vs Temperature (PVT Simulation)

import neqsim.pvtsimulation.simulation.WaxFractionSim;

WaxFractionSim waxSim = new WaxFractionSim(oil);
waxSim.setTemperatures(new double[]{333.15, 313.15, 293.15, 273.15});
waxSim.run();
double[] waxFractions = waxSim.getWaxFraction();

3. Asphaltene Stability

de Boer Screening

// Assess asphaltene precipitation risk
// Key parameters: reservoir pressure, bubble point, density difference
// Risk increases when operating pressure approaches bubble point
// High-risk zone: ΔP > 200 bar above bubble point for light oils

// Use CPA for asphaltene modeling
SystemInterface aspFluid = new SystemSrkCPAstatoil(273.15 + 90, 300.0);
// Add components including heavy asphaltenic fractions
aspFluid.setMixingRule(10);

ThermodynamicOperations ops = new ThermodynamicOperations(aspFluid);
ops.TPflash();
aspFluid.initProperties();

// Check if asphaltene phase is stable
// Compare upper/lower asphaltene onset pressures vs operating P

4. Pipeline Hydraulics

Simple Adiabatic Pipe

AdiabaticPipe pipe = new AdiabaticPipe("Export Pipeline", feedStream);
pipe.setLength(50000.0);       // 50 km in meters
pipe.setDiameter(0.508);       // 20 inch in meters
pipe.setInletElevation(0.0);
pipe.setOutletElevation(-350.0);  // negative = downhill (subsea)
pipe.run();

double outletP = pipe.getOutletStream().getPressure();  // bara
double outletT = pipe.getOutletStream().getTemperature() - 273.15;  // C
double dP = feedStream.getPressure() - outletP;  // pressure drop

DNV-RP-F105 free-span routing

When a hydraulics or environment study feeds an explicit current DNV-RP-F105 2025-12 free-span screen, route verified structural and environmental inputs through DnvRpF105FreeSpanScreeningKernel. Keep steel and hydrodynamic diameters distinct and use velocities normal to the span. Effective modal mass, axial force, span geometry, and response-trigger basis are external structural inputs, not quantities inferred silently from a hydraulic pipe object.

The kernel is a simply supported first-mode/dimensionless escalation screen. Its Strouhal number, frequency-ratio band, and reduced-velocity triggers are project-controlled and cannot be called DNV limits. Keep soil/shoulder and multi-span response, VIV amplitudes, direct wave loading, ULS/FLS, fatigue, monitoring, and intervention external. Never relabel PipeMechanicalDesignCalculator.calculateAllowableSpanLength(...) as F105 evidence.

DNV-RP-F101 inspected metal-loss routing

When inspection data feeds an explicit current DNV-RP-F101 2019-09+AMD:2025-09 screen, route one verified isolated longitudinal metal-loss defect under internal pressure through DnvRpF101CorrodedPipelineScreeningKernel. Require the measured depth and axial length, assessment wall-thickness definition, inspection/growth allowance, characteristic ultimate tensile strength, internal/external absolute pressures, and project-controlled pressure factor.

Do not infer defect geometry from hydraulic corrosion-rate calculations or projected uniform wall loss. The typed kernel does not handle defect interaction or complex profiles, longitudinal compression, probabilistic assessment, crack-like damage, repair, or fitness-for-service approval. It also does not replace DNV-ST-F101 original-design checks.

DNV-RP-F104 CO2 pipeline routing

When an actual-composition phase-envelope and hydraulic/thermal study feeds a current DNV-RP-F104 2021-02+AMD:2021-09 screen, route the bounded project composition, CO2/water limits, ordered profile, absolute MAOP, design temperatures, and a separately verified minimum single-phase pressure boundary at each point through DnvRpF104Co2PipelineEnvelopeScreeningKernel.

The external thermodynamic basis must establish that pressure above each boundary represents the intended single-phase region for the specific composition, temperature, path, EOS, and uncertainty. Do not substitute pure-CO2 critical conditions or CO2FlowCorrections.isDensePhase(...). Treat composition, phase-boundary, MAOP, and temperature margins as screening findings. Keep transient cases, F104 decompression/fracture and crack arrest, materials/corrosion, release consequences, construction, operation, requalification, and all DNV-ST-F101 structural checks external.

DNV-RP-F114 pipe-soil interaction routing

When route, hydraulic/thermal, or installation work feeds a current DNV-RP-F114 2021-05 screen, route named design situations with externally verified vertical, axial, and lateral action and resistance magnitudes through DnvRpF114PipeSoilInteractionScreeningKernel. Treat margin and utilization outputs as caller-controlled screening findings.

Do not convert burial depth, soil thermal resistance, or a generic friction factor into geotechnical resistance. Keep site investigation, soil interpretation, penetration/burial and load-displacement response, time/cyclic effects, characteristic values, uncertainty, structural actions, and F109/F110/F105/ST-F101 acceptance external.

DNV-RP-F110 global-buckling response routing

When hydraulic/thermal, route, or installation work feeds a current DNV-RP-F110 2019-09+AMD:2021-09 screen, route named external structural-analysis cases through DnvRpF110GlobalBucklingResponseScreeningKernel. Supply effective force, peak longitudinal strain, peak global displacement, and required feed-in length with caller-controlled allowable or available values. Treat margins and utilizations as screening findings.

Require external evidence for the effective-force derivation, pipe/as-laid geometry, pipe-soil response, imperfections/triggers/strategy, global structural model, load combinations, local capacity and strain criteria, uncertainty/sensitivity/buckle sharing, and lifecycle actions. Never derive critical buckling, initiation/prevention criteria, structural response, or soil springs from NeqSim hydraulic or thermal output. Keep F109/F114/F105 and all ST-F101 acceptance external.

Beggs and Brill Multiphase Correlation

PipeBeggsAndBrills pipe = new PipeBeggsAndBrills("Subsea Flowline", feedStream);
pipe.setPipeWallRoughness(5e-5);  // meters
pipe.setLength(50000.0);          // meters
pipe.setAngle(0.0);               // horizontal
pipe.setDiameter(0.254);          // 10 inch

// For subsea with heat loss
pipe.setOuterTemperature(277.15);  // 4°C seawater
pipe.run();

// Get flow regime, liquid holdup, pressure profile
double outP = pipe.getOutletStream().getPressure();
double outT = pipe.getOutletStream().getTemperature() - 273.15;

Liquid Holdup, Flow Regime & Liquid-Loading (gravity-dominated screening)

Verified reader methods on PipeBeggsAndBrills (after run()):

String regime = pipe.getFlowRegime();          // SEGREGATED / TRANSITION / INTERMITTENT / DISTRIBUTED
double dP     = pipe.getPressureDrop();         // bar (inlet - outlet)
double vmix   = pipe.getMixtureVelocity();      // m/s
double[] holdupProfile = pipe.getLiquidHoldupProfile();   // fraction per segment (0-1)
// per-segment access (0 .. numberOfIncrements-1):
Double hSeg   = pipe.getSegmentLiquidHoldup(i);
Double vsgSeg = pipe.getSegmentGasSuperficialVelocity(i);
Double vslSeg = pipe.getSegmentLiquidSuperficialVelocity(i);
Double elevSeg = pipe.getSegmentElevation(i);
// average holdup = mean(holdupProfile); "liquid content" for the line
// liquid inventory (m3) = sum_i holdup_i * (pi/4*D^2) * segmentLength

Liquid-loading / gravity-dominated screening (PEPR-style "is the line filling with liquid?"): sweep gas rate (and water cut) and read average holdup + regime. As gas rate falls, holdup rises and the regime moves INTERMITTENT -> TRANSITION -> SEGREGATED (stratified) = gravity-dominated / liquid loading. Higher water cut lifts holdup at every rate and pushes toward INTERMITTENT (slugging). Define liquid-loading onset as the gas rate where average holdup crosses a threshold (e.g. 25%).

for (double qgMSm3d : gasRates) {
  SystemInterface feed = fluidTemplate.clone();      // gas-condensate + water (CPA, rule 10)
  Stream s = new Stream("feed", feed);
  s.setFlowRate(qgMSm3d, "MSm3/day");                // wellstream gas standard volume
  s.setTemperature(50.0, "C"); s.setPressure(150.0, "bara");
  s.run();
  PipeBeggsAndBrills p = new PipeBeggsAndBrills("line", s);
  p.setLength(21000.0); p.setDiameter(0.254); p.setAngle(0.0);
  p.setPipeWallRoughness(5e-5); p.setNumberOfIncrements(40);
  try {
    p.run();
    double[] h = p.getLiquidHoldupProfile();
    // record mean(h), p.getFlowRegime(), p.getMixtureVelocity()
  } catch (RuntimeException e) {
    // "Outlet pressure is negative" = DELIVERABILITY LIMIT, not a bug (see gotcha)
  }
}

GOTCHA — deliverability limit vs bug. PipeBeggsAndBrills uses a fixed inlet pressure. If frictional ΔP over a long/small line exceeds the inlet pressure, run() throws InvalidOutputException: ... Outlet pressure is negative. That is a genuine deliverability limit (the line cannot pass that rate at that inlet P), not a solver failure — catch it and report the max deliverable rate. To model to a fixed arrival (outlet) pressure instead, raise the inlet pressure until the delivered rate matches, or iterate inlet P per rate.

GOTCHA — getFlowRegime() naming. Beggs & Brill regimes are returned as SEGREGATED (stratified/annular — gravity-dominated), TRANSITION, INTERMITTENT (plug/slug), DISTRIBUTED (bubble/mist). "Gravity-dominated / liquid loading" = SEGREGATED (+ low-velocity TRANSITION).

GOTCHA — profile index vs getFlowRegime(). run() evaluates the correlation once per increment and then once more at the outlet state, so every correlation profile (getLiquidHoldupProfile(), getFlowRegimeProfileList(), …) has numberOfIncrements + 1 entries and index i belongs to getPressureProfile()[i]. getFlowRegime() returns the outlet state. Reading getFlowRegime() next to getSegmentLiquidHoldup(0) compares two different states and can look like a discontinuity. Always pair getSegmentFlowRegime(i) with getSegmentLiquidHoldup(i).

VALIDITY — low liquid loading. On a long wet-gas export line (ID 0.355 m, no-slip liquid fraction λ_L ≈ 0.006), Beggs & Brill over-predicts ΔP well above a single-phase Darcy-Weisbach integration of the same line, and the gap is entirely the two-phase friction multiplier exp(S): removing it brings ΔP within 4% of that analytic integration. S is monotonically increasing in y = λ_L / H_L², so less liquid gives a larger multiplier — the opposite of the physical trend — up to a bounded maximum of exp(S) = 3.19 at y = 52.1. B&B is calibrated for λ_L down to roughly 0.01–0.02; below that, use TwoFluidPipe or a mechanistic simulator and treat B&B as a conservative upper bound. The published map also has a genuine step where the segregated and distributed correlations meet at L1 for λ_L < 0.01 (no transition band exists there): hold-up ×0.70 and ΔP ×1.12 across a 1 bara change. That step is in the correlation — do not smooth it.

Fixed defects worth knowing (all affected the inclination correction). Older NeqSim builds silently returned the horizontal hold-up on an inclined leg when (a) the pipe angle was converted degrees→radians twice, (b) the Baker-Swerdloff surface tension went negative above 274 bara or for a very light liquid, making the liquid velocity number NaN, or (c) the flow regime was TRANSITION, which had no inclination branch at all. All three are fixed and locked by PipeBeggsAndBrillsCorrelationTest. If you are on an older build, sanity-check that uphill hold-up clearly exceeds horizontal hold-up before trusting an inclined result.

Phase-count code paths. PipeBeggsAndBrills assumes phase 0 is the gas whenever the stream has more than one phase. A gas-free stream that splits into oil and water (dead oil with free water at high pressure) used to be modelled as gas–liquid, with the oil phase acting as the gas: fictitious flow regime, hold-up 0.21 instead of 1, and ΔP 44% above the homogeneous liquid value. It is now carried as a homogeneous liquid with a volume-weighted density and viscosity. The three-phase liquid density is also now combined on volume fractions (total mass over total volume) rather than mass fractions. Locked by PipeBeggsAndBrillsPhasePathsTest.

Where the error sits. On a single-phase gas line Beggs & Brill matches a Darcy-Weisbach integration to a few tenths of a per cent on pressure drop and to ~0.15 K on arrival temperature — the friction and energy paths are sound. The over-prediction seen on the two-phase version of the same line is therefore entirely the two-phase friction multiplier, not the solver.

Comparisons against commercial transient multiphase simulators are not published in this repository. Their licence terms generally prohibit publishing benchmark results and prohibit using the software to develop competing software, so no NeqSim closure is tuned to such a tool and no measured deviation against one is recorded here. Validate against experimental, laboratory or field data, or against an analytic/first-principles check.

TwoFluidPipe steady-state usage. The refinement loop needs a few hundred sweeps to settle the pressure profile against the updated section densities (74 km at 160 sections ≈ 25 s, at 320 sections ≈ 60 s). Always assert pipe.isSteadyStateConverged() and check pipe.getSteadyStateIterationsUsed() > 1 — a result reported after one sweep still carries the densities the sections were initialised with and understates the pressure drop of a gas line by roughly ten per cent. Around 160 sections (≈450 m) is the sweet spot on a long transmission line: it converges reliably and is grid-converged to 0.4% against 320 sections. Raise setSteadyStateMaxWallClockTime(...) (default 300 s) before blaming the model if isSteadyStateWallClockLimited() is true.

MCP runPipeline solver and response handoff. Omit solver (or use beggsBrill) for the established correlation path. Use solver: "twoFluid" when the task needs finite-volume pressure, temperature, holdup, phase-velocity, flow-regime, inventory, erosion-margin, hydrate/wax, or slug profiles. Pass sectionLengths_m, elevationProfile_m, heatTransferProfile_W_m2K, and surfaceTemperatureProfile_C or _K as equal-length per-section arrays; the lengths must sum to length_m. Bound expensive solves with steadyStateMaxWallClockTime_s. Request detailLevel: FULL for spatial profiles, SUMMARY for engineering KPIs without profiles, or MINIMUM for the compact core result. The authoritative object is TwoFluidPipeResponse (a BaseResponse): preserve it through MCP/report handoffs instead of independently reconstructing fields from the equipment. Always propagate its convergence/pressure-floor/wall-clock validation findings and apply the limitations below before quoting a result.

Three defects found by auditing the solve against its own governing equations — all fixed, but the symptoms are generic and worth recognising elsewhere. (1) A time integrator inside a fixed-point sweep. The steady solve integrated LiquidAccumulationTracker once per iteration with a nominal dt; that tracker only ever adds liquid, ratchets its volume up to what the sections already hold, then adds it back on top of that holdup, so it has no fixed point. Valley sections climbed to the 0.85/0.95 cap and the profile never settled. (2) A correlation overriding a solved momentum balance. The minimum-slip constraint applied the Beggs and Brill horizontal holdup correlation as a lower bound in every regime. It is fitted to 1–1.5 inch air-water loops at near-atmospheric pressure with λ_L ≥ 0.01; the line runs λ_L ≈ 0.008 in a 14-inch pipe at 200 bara, and the floor was binding in EVERY section — so the reported holdup was the correlation, not the solved momentum balance, worth ≈+20% on ΔP through the mixture density. Only the scale-free lambdaL * minimumSlipFactor bound remains. (3) Convergence declared without re-evaluating thermodynamics. The flash runs every ssFlashInterval sweeps but the "flash moved nothing" flag started false, so a non-flash sweep read as settled. The solve could exit after ONE sweep on densities it had never revisited. If you see a steady profile with a section sitting exactly on a cap, or getSteadyStateIterationsUsed() returning 1 on a long line, suspect these.

Behaviour after those fixes, on that line at default settings: the rate exponent in ΔP ~ rate^n rises from about 2.1 to about 3.1 across a 3x rate range, i.e. the density feedback along the line is reproduced, not just the level at one rate; 10 MW of DEH raises arrival T by ~17 K and ΔP by ~15%, so the energy equation feeds the momentum balance. Grid-converged (160 vs 320 sections within 0.4%). Terrain response is solved, not tuned: the annular film balance carries tau_i = tau_wL + rhoL*g*sin(theta)*delta, so holdup responds to inclination and scales with sin(theta). Still open: the three-phase free-water case does not converge — 15 m3/hr of free water on the same line is wall-clock limited after 4078 iterations at a 1200 s budget. ΔP is identical between a 300 s and a 1200 s budget, so the profile is stationary and the criterion is stalling on the three-phase liquid split — but isSteadyStateConverged() is false, so the number must not be quoted. ALWAYS check it on a water-bearing line.

TwoFluidPipe also fails silently when a line has no deliverability. The marching solver clamps section pressure at a 1 bara floor; that clamp is a fixed point of itself, so the per-section change falls below tolerance and the sweep would report success on a case with no physical solution. Check pipe.isSteadyStatePressureFloorLimited() — when it is true, isSteadyStateConverged() is withheld and the profile must be discarded, not reported. PipeBeggsAndBrills throws Outlet pressure is negative on the same condition.

TwoFluidPipe holdup at low rate is dominated by terrain trap sections. The MAXIMUM holdup at 4 and 7 MSm3/d sits well above the line mean (a single valley section), so do not quote local TwoFluidPipe holdup or valley inventory as a design number. The three-phase bookkeeping is sound — gas/oil/water fractions sum to one and stay in range at every node.

The minimum-slip bound applies only level and uphill. alphaL >= lambdaL * minimumSlipFactor (default 2.0) states that the gas outruns the liquid, which is a property of gas-driven transport. On a downhill section gravity moves the liquid and the slip ratio legitimately falls, so the bound is not applied there. It used to be applied everywhere and was binding on 39 of 42 downhill sections of an undulating fixture, replacing the momentum balance with a constant.

The horizontal annular criterion is the Taitel-Dukler equilibrium level. pipe.setUseEquilibriumLevelAnnularTransition(false) restores the earlier path, the vertical droplet-entrainment threshold U_SG > 3.1*(sigma*g*drho/rhoG^2)^0.25, which is about 0.75 m/s on a 14-inch export line and so classified essentially any horizontal gas pipeline as annular, solving a shallow stratified layer with a thin-film closure. The two agree wherever the Kelvin-Helmholtz margin exceeds one — identical at 10 MSm3/d on the export line — and differ at 4 MSm3/d, where the equilibrium branch reclassifies 272 of 320 sections as stratified-wavy.

Friction is per-phase wall shear in stratified flow. setSeparatedFrictionModel(false) restores the mixture correlation, which charges the whole perimeter with a holdup-weighted density; on a stratified line at 41% holdup that over-predicts ΔP by ~2.3x, and because it scales as G^2/rho_mix it makes extra liquid REDUCE the gradient, inverting the terrain response. The separated form is scoped to stratified flow because its perimeters come from a circular-segment layer; annular flow, whose film wets the whole perimeter, is not that geometry.

Measured accuracy on the 73.8 km reference line, across a threefold rate range: ΔP +1.4 / +1.6 / +0.1 / −2.7 % and maximum holdup −2.4 / −7.1 / −6.6 / +3.5 % at 4 / 7 / 10 / 12 MSm3/d, all converged and grid-converged. The earlier defaults gave ΔP +5.7 / +5.6 / +1.4 / −0.0 % and holdup −25.5 / −18.6 / −6.2 / +3.3 %.

Direct electrical heating (DEH) is available on both models with the same convention — the power set is what reaches the fluid, so cable and coating losses must already be deducted: pipe.setDirectElectricalHeatingPower(watts) or setDirectElectricalHeatingPowerPerMeter(wattsPerMetre). In TwoFluidPipe steady state each segment decays toward the balance temperature T_surface + q/(U·π·D) (exact for a uniform source, cannot overshoot); it also works in transient and with wall heat transfer switched off. A tool with no distributed-heating input can represent the same source through the identity −UπD(T−T_surf) + q ≡ −UπD(T − [T_surf + q/(UπD)]), i.e. by raising the ambient temperature by q/(UπD).

TwoFluidPipe liquid-rich and severe-slugging transients remain unqualified. The legacy route can still develop phase backflow and clamp the outlet flux at zero while its finite-volume balance closes exactly. Gas-dominated null cases remain usable, but do not promote a liquid-rich or severe-slugging trajectory from either a small residual or a completed time loop alone. Use the analytical SevereSluggingBenchmarkHarnessTest screen and the public Tengesdal evidence until the dynamic qualification gates below pass.

The coupled route is an opt-in four-part configuration. For a pressure outlet that physically permits phase fallback, use setEnableInterfacialPressure(true), setImplicitInterfacialPressureCoupling(true), setEnableCoupledPressureMomentum(true), and setAllowOutletPhaseBackflow(true) together. The nonlinear controls are public: setCoupledPressureMomentumMaximumIterations(int) and setCoupledPressureMomentumRelativeVolumeTolerance(double), with defaults 24 and 1e-7.

WS3 restores progress, not physical parity. The 16-section Tengesdal Test 3 probe now completes 50/50 calls of 0.1 s; a 24-section refinement completes 100/100 calls of 0.05 s. Neither rejects a nonlinear substep, and gas, oil, water, liquid, and total mass residuals are below 1e-9. The former 12-iteration cap stopped around 6e-7, above the 1e-7 gate. The sticky pressure limiter still fires, however, and the liquid outlet spans -18.55 to 6.88 kg/s versus the stored 0.375 to 4.03 kg/s comparison. This is a disclosed boundary/ pressure-coupling gap, not a reason to tune a public closure to a commercial trace. Subsequent validation must use the public Tengesdal experiment, nearby operating points, conservation, and mesh/time-step refinement.

What did help: fixing the regime. The same case classified SLUG rather than ANNULAR (PR #3086, Barnea bridging limit) cuts the 30-minute inventory runaway from +56.3% to +20.1% at default settings. Still unusable, still gated by the flag, but the regime branch is a bigger lever on the runaway than the interfacial-pressure term is.

The transient tells you when it has failed — always check every diagnostic. Read isTransientOutletBackflowClamped(), isTransientCoupledPressureMomentumFailureDetected(), isTransientCoupledPressureMomentumCorrectionLimited(), and getTransientCoupledPressureMomentumRejectedSubsteps() after the full window. isCoupledPressureMomentumPressureCorrectionLimited() is the non-sticky view of only the latest correction. The flags/counter are sticky and reset on the next steady run(). A coupled call that cannot complete its requested interval now throws with accepted time, requested time, residual/tolerance, iterations/cap, and limiter state; never catch it and advance the engineering timeline. The mass-balance report alone cannot qualify the result because a clamped or limited route can still conserve its own discrete fluxes exactly.

Three-phase (gas/oil/water) steady state is fixed. The oil/water slip ratio uses S = 1 + 1.75·max(0, 1 − (Fr/3)²), a stratified plateau that rolls off to no slip once the liquid disperses above a liquid Froude number of about 3; the previous form cut off at Fr = 2 and under-predicted water holdup badly. One gap remains open: the pressure drop is over-predicted in this liquid-rich regime, far more than on a gas-dominated line.

Exporting NeqSim to OLGA — the fluid basis is two files, not one. The .tab PVT table fixes the phase behaviour; the hydrate boundary is separate and OLGA does not compute it. Without a HYDRATECURVE in the case, OLGA falls back to the Hammerschmidt correlation, so a study whose NeqSim half uses CPA hydrate equilibrium and whose OLGA half uses Hammerschmidt disagrees about where hydrates form, invisibly. Export both from the same fluid: OLGAhydrateCurveGenerator writes the HYDRATECURVE LABEL=..., PRESSURE=(...) bara, TEMPERATURE=(...) C block and returns the matching HYDRATECHECK HYDRATECURVE="..." line for the flowpath. OLGA interpolates that curve linearly, so span the pressures the case actually visits and use ≥20 points when the range reaches below ~50 bara (4 points over 10–200 bara costs 4.1 K of hydrate temperature; 20 points costs 0.48 K). The OLGA output variable is DTHYD in °C, and it is T_hydrate − T_fluid — positive means inside the hydrate region, negative is the safe margin, which is the opposite of the intuitive reading. Full OLGA-side workflow in the community neqsim-olga-multiphase-simulator skill.

Never build a volumetric phase fraction from phase.getVolume(). With a Peneloux volume shift active it disagrees with getDensity() by the shift — +16.6% for oil and +31.7% for water on a typical SRK three-phase system, while gas matches to 0.25%. Use getFlowRate("kg/sec")/getDensity("kg/m3").

advection-relaxation transport lag, not a conservation-law solver, and the Beggs & Brill correlation is not even used on that path (friction reverts to single-phase Darcy-Weisbach, viscosity is frozen at the inlet). It does not store mass: on a rate step the outlet mass flow equals the inlet at every timestep, so ∫(ṁ_in − ṁ_out)dt = 0 while the inventory implied by its own profile moves 171 t — about 192 t of gas appears from nowhere on a 74 km line. This is not repairable in the class as written: it takes only an inlet boundary condition, so there is nothing to pin the arrival pressure and drive line pack. It also does not exactly preserve its own steady state — with the boundary conditions held constant it drifts −6.3 bar on that line (was +30 bar before the cell-density fix), because the transient friction closure differs from the steady one. Note also that calculateSteadyState defaults to true, so without setCalculateSteadyState(false) runTransient is only a steady-state solve with the clock advanced; and the time step must be shorter than the segment transit time L/numberOfIncrements/v, otherwise the relaxation factor saturates at 1 and the whole line responds in a single step. Use it for transport delay in a flowsheet; use TwoFluidPipe for line pack, shut-in, ramps and blowdown (0.00 bar drift on the same null test, mass balance closing to the digit), and WaterHammerPipe for surge.

Gray (1974) Correlation — Gas / Gas-Condensate Vertical Wells

PipeGray implements the Gray (1974) correlation, the industry standard for gas-dominated vertical wells producing condensate and/or water (API 14B program). Prefer it over Beggs & Brill for vertical/near-vertical gas-condensate tubing where the superficial gas velocity is high (> ~4.6 m/s), the tubing is small (< ~3.5 in), and condensate loading is low (< ~50 bbl/MMscf). It predicts in-situ liquid holdup and a condensate-film effective roughness.

PipeGray well = new PipeGray("Gray well", inletStream);   // gas-condensate wellstream
well.setDiameter(0.0889);        // 3.5 inch tubing
well.setLength(3000.0);
well.setElevation(3000.0);       // vertical well (upward flow)
well.setNumberOfIncrements(10);
// Optional: swap the holdup closure to Woldesemayat-Ghajar (2007)
well.setHoldupMethod(PipeGray.HoldupMethod.WOLDESEMAYAT_GHAJAR);
well.run();

double dP       = well.getTotalPressureDrop();          // bar
double holdup   = well.getLiquidHoldup();               // fraction (0-1)
double vsg      = well.getSuperficialGasVelocity();     // m/s
double ke       = well.getEffectiveRoughness();         // m (Gray condensate-film roughness)

Single-phase gas and single-phase liquid segments fall back to a Haaland friction-factor Darcy-Weisbach drop, so the same model spans wet-gas wells that drop out condensate along the tubing. The void-fraction closures are exposed directly via VoidFractionCorrelations.woldesemayatGhajar(...).

Pipeline with Formation Temperature Gradient (Wells / Risers)

PipeBeggsAndBrills wellbore = new PipeBeggsAndBrills("Production Well", feedStream);
wellbore.setLength(3000.0);
wellbore.setElevation(-3000.0);  // vertical well
wellbore.setDiameter(0.1571);    // 6-5/8 inch
wellbore.setPipeWallRoughness(5e-5);

// Formation temperature: 90°C at bottom, gradient of -0.03°C/m going up
wellbore.setFormationTemperatureGradient(4.0, -0.03, "C");
wellbore.run();

Pipeline Sizing (Iterative)

// Iterate over diameters to find optimal size
double[] diameters = {0.1524, 0.2032, 0.254, 0.3048, 0.3556, 0.4064, 0.508};
// 6", 8", 10", 12", 14", 16", 20"

for (double d : diameters) {
    Stream testFeed = new Stream("feed", feedFluid.clone());
    testFeed.setFlowRate(flowRate, "kg/hr");
    testFeed.run();

    PipeBeggsAndBrills testPipe = new PipeBeggsAndBrills("test", testFeed);
    testPipe.setLength(pipeLength);
    testPipe.setDiameter(d);
    testPipe.setPipeWallRoughness(5e-5);
    testPipe.run();

    double dP = testFeed.getPressure() - testPipe.getOutletStream().getPressure();
    double velocity = testPipe.getSuperficialVelocity();
    // Check: erosional velocity < API RP 14E limit, dP within allowable
}

5. CO2 / H2S Corrosion Assessment

CO2 Partial Pressure Based Screening

// After flash calculation
fluid.initProperties();
double pCO2 = fluid.getPhase("gas").getComponent("CO2").getx()
              * fluid.getPressure();  // CO2 partial pressure in bara

// NORSOK M-001 / DNV-RP-F112 screening:
// pCO2 < 0.02 bar  → low risk (carbon steel OK)
// 0.02 < pCO2 < 0.2 → moderate (corrosion allowance or inhibitor)
// pCO2 > 0.2        → high risk (CRA or heavy inhibition)

Temperature and pH Effects

// Corrosion rate increases with temperature up to ~80°C
// then decreases due to protective FeCO3 film
// Lower pH (more acidic) → higher corrosion rate
// Water cut affects wetting: >30% water cut → higher risk

Network-Level Corrosion (LoopedPipeNetwork)

For production gathering networks, LoopedPipeNetwork has inline corrosion models that compute rates per element during network solution:

// de Waard-Milliams (default) or NORSOK M-506
net.setCorrosiveGas("trunk", 0.035, 0.002);  // CO2 mol%, H2S mol%
net.setCorrosionModel("trunk", "NORSOK");     // "DEWAARD" or "NORSOK"
net.setMinAllowableWallLife(20.0);            // years

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

Models: de Waard-Milliams (log10(Vcorr) = 5.8 - 1710/T + 0.67*log10(pCO2)) and NORSOK M-506 (Vcorr = Kt * fCO2^0.62 * (S/19)^0.146).

Rigorous NORSOK M-506 from a brine (electrolyte pH + FeCO3 film)

NorsokM506CorrosionRate (neqsim.process.corrosion) is the standalone NORSOK M-506 model (fugacity, in-situ pH, FeCO3 scaling temperature, wall shear, glycol, inhibitor). By default it estimates pH from a CO2-water correlation. To drive it from rigorous electrolyte thermodynamics instead, use NorsokM506ElectrolyteBridge — the M-506 analogue of CO2CorrosionAnalyzer:

// fluid = SystemElectrolyteCPAstatoil with CO2, water, ions (+ optional Fe++)
NorsokM506ElectrolyteBridge bridge = new NorsokM506ElectrolyteBridge(fluid);
bridge.setFlowVelocityMs(3.0);
bridge.setPipeDiameterM(0.254);
bridge.run();                       // clones + flashes; does not mutate `fluid`
double rate = bridge.getModel().getCorrectedCorrosionRate();  // mm/yr
double pH   = bridge.getInSituPH();          // rigorous aqueous pH (getpH())
double sr   = bridge.getFeCO3SaturationRatio();  // IAP/Ksp, or -1 if no Fe++/CO3--

FeCO3 film feedback (closes the corrosion↔scaling loop): when the aqueous phase is supersaturated in siderite (SR>1), a protective film suppresses corrosion even below the scaling temperature. Supply the ratio directly on the model, or let the bridge compute it from aqueous Fe++/CO3-- molalities (Sun & Nesic 2009 Ksp):

model.setFeCO3SaturationRatio(sr);   // >1 = protective; -1 = disabled (default)
double film = model.calculateFeCO3FilmFactor();  // 1.0 = no credit, <1 = protective

For an auditable standards calculation, pass the resulting in-situ pH and optional saturation ratio to NorsokM506CorrosionDesignKernel.Input. The kernel is the preferred public path when the task names NORSOK M-506: it enforces the unamended 2017 edition and model envelope, preserves raw inputs without setter clamping, and returns CALCULATED_REVIEW_REQUIRED. Continue to use the bridge to derive chemistry and the legacy model for sweeps, but do not present either route as a conformity assessment. The optional saturation-ratio film factor and projected wall loss are NeqSim screening extensions, not code acceptance criteria.

Gotchas (verified): SystemInterface.clone() drops the chemical-reaction setup — re-run chemicalReactionInit() on the clone before flashing or the CO2-brine pH comes out unphysically basic (~10). setActualPH() is read back via getEffectivePH(), NOT getCalculatedPH(). Proven converging brine: CO2 0.10 / water 0.88 / Na+ 0.01 / Cl- 0.01, setMixingRule(10).

Robust in-situ pH for investigations

getpH() on the aqueous phase (or SystemInterface.getpH()) returns the in-situ pH. It now has a built-in acid-gas dissociation fallback: when the aqueous phase has water + dissolved CO2/H2S but no explicit H3O+ species (i.e. chemicalReactionInit() was not run) or the electrolyte solver is unstable and returns NaN at low pressure, getpH() estimates pH from the carbonic/sulfide acid first-dissociation equilibria — [H+]=sqrt(K1_CO2·C_CO2 + K1_H2S·C_H2S + Kw) — instead of the old silent, unphysical 7.0. CO2-saturated water → pH ≈ 3.9. Force it explicitly with getpH("acidgas"). This is a screening estimate (ignores bicarbonate buffering and salt-ion alkalinity); for a rigorous speciated pH in a buffered brine still run chemicalReactionInit().

For corrosion/scale investigations that must always return a finite, bounded, source-tagged value with an explicit pCO2 basis, use RobustAqueousPH.estimate(fluid, pCO2Bar): it takes the rigorous getpH() when valid, else a CO2-water correlation, and records which source was used (getSource(), isFellBack()).

Regression test for the fallback: neqsim.chemicalreactions.AcidGasPHTest (CO2-saturated water pH ≈ 3.9, neutral water ≈ 7, CO2+H2S acidic — all without chemicalReactionInit()).

Per-segment corrosion + scale along a line (PipeSegmentIntegrity)

PipeSegmentIntegrity (neqsim.process.corrosion) walks a temperature/pressure/ velocity profile and reports, per segment, the NORSOK M-506 corrosion rate AND the CaCO3 scale saturation, then ranks the worst corrosion and worst scale segments — so mitigation can be located along the line. Feed it a run PipeBeggsAndBrills:

PipeSegmentIntegrity integrity = new PipeSegmentIntegrity();
integrity.fromPipe(pipe)                       // extracts T, P, mixture velocity, diameter
    .setPipeAndGas(0.2, 0.05)                  // diameter [m], CO2 mole fraction
    .setBrineChemistry(1500, 400, 0, 0, 12000, 35000); // Ca,HCO3,SO4,Ba,Na,Cl [mg/L]
integrity.evaluate();
int worstCorr = integrity.getWorstCorrosionIndex();
int worstScale = integrity.getWorstScaleIndex();
// or supply an explicit profile: integrity.setProfile(tC[], pBara[], vMs[])

Network-Level Sand Erosion (DNV RP O501)

net.setSandRate("W1", 3.0);             // kg/hr
net.setMaxAllowableSandRate(10.0);
net.setMaxAllowableErosionRate(5.0);    // mm/yr

net.run();
Map<String, double[]> sand = net.calculateSandTransport();
// Per element: [0] = rate, [1] = concentration, [2] = erosion, [3] = deposition
List<String> violations = net.getSandViolations();

Erosion per DNV RP O501: E = K * Csand * v^2.6 * dp^0.2. Deposition flagged when v < 1 m/s.

Mineral Scale (thermodynamics + kinetics + brine mixing)

NeqSim offers two routes to mineral scale. Pick by how much you know about the brine and whether you need a rigorous precipitated amount.

(A) Screening SI from an ion analysis (fast, no flash). Activity-corrected saturation index (Davies + Ksp(T)) for the common scales directly from a produced-water ion table in mg/L:

ElectrolyteScaleCalculator scale = new ElectrolyteScaleCalculator();
scale.setTemperatureCelsius(60).setPressureBara(100).setPH(6.0);
scale.setCations(caMgL, baMgL, srMgL, mgMgL, naMgL, kMgL, feMgL);
scale.setAnions(clMgL, so4MgL, hco3MgL, co3MgL);
scale.calculate();
double siCaCO3 = scale.getCaCO3SaturationIndex();  // also BaSO4, CaSO4, SrSO4
// SI = log10(IAP/Ksp): >0 supersaturated (scale risk), 0 at equilibrium, <0 undersaturated

(B) Rigorous scale potential from a speciated electrolyte fluid — REQUIRES chemical reactions + speciation. The rigorous saturation ratio needs the in situ ion molalities (Ca²⁺, HCO₃⁻, CO₃²⁻, Ba²⁺, SO₄²⁻, …), which only exist after the aqueous speciation equilibria have been solved. So you MUST call chemicalReactionInit() before createDatabase(true) / setMixingRule(10), and you should get the in-situ pH from the same speciated flash (see the in-situ pH section above). Then checkScalePotential(phaseNumber) returns, per salt, the saturation ratio SR = IAP/Ksp ("relative solubility"): SR > 1 means supersaturated and at scale risk.

SystemInterface brine = new SystemElectrolyteCPAstatoil(273.15 + 60.0, 100.0);
brine.addComponent("CO2", 0.02);
brine.addComponent("water", 55.5);
brine.addComponent("Na+", 1.0);
brine.addComponent("Cl-", 1.0);
brine.addComponent("Ca++", 0.02);
brine.addComponent("HCO3-", 0.04);
brine.chemicalReactionInit();     // MANDATORY: builds carbonate/bicarbonate speciation
brine.createDatabase(true);
brine.setMixingRule(10);          // numeric CPA rule (required)
brine.setMultiPhaseCheck(true);   // allow gas + aqueous (+ solid) phases

ThermodynamicOperations ops = new ThermodynamicOperations(brine);
ops.TPflash();                    // solves speciation in the aqueous phase
brine.initProperties();

int aq = brine.getPhaseNumberOfPhase("aqueous");
ops.checkScalePotential(aq);      // uses speciated ion molalities in phase `aq`
String[][] sr = ops.getResultTable();   // rows: {saltName, SR (=IAP/Ksp), ""}
// Read SR per salt; SR > 1.0 => supersaturated => scale precipitation likely.
double pH = brine.getpH();        // rigorous speciated in-situ pH from the same flash

Convenience equipment wrapper: ScalePotentialCheckStream runs the same SR check on a process Stream.

Precipitation amount (how much solid forms). Two options:

• Rigorous (solid phase flash): with setMultiPhaseCheck(true), a supersaturated brine drops a solid phase during TPflash(). Read the precipitated mineral directly:

  if (brine.hasPhaseType("solid")) {
    double solidMoles = brine.getPhase("solid").getNumberOfMolesInPhase();
    double solidMassKg = solidMoles * brine.getPhase("solid").getMolarMass();
  }
  

• Screening (from SI): ScaleMassCalculator estimates the mg/L precipitated to reach equilibrium from ion concentrations and SI (per mineral):

  ScaleMassCalculator mass = new ScaleMassCalculator(predictionCalc);
  mass.setWaterVolume(1.0); // litres
  double caco3MgL = mass.calcCaCO3Mass(cCaMolL, cCO3MolL, siCaCO3);  // mg/L
  // also calcBaSO4Mass, calcCaSO4Mass, calcSrSO4Mass, calcFeCO3Mass
  

  ⚠️ ScaleMassCalculator is decoupled — each mineral is treated independently, so it overpredicts competing minerals (e.g. it reports celestite scaling even when barite has already consumed the shared sulphate).

• Coupled screening (recommended for shared-ion brines): MultiMineralScaleEquilibrium drops all supersaturated minerals simultaneously and re-equilibrates the shared ion pools (sulphate: BaSO4 / SrSO4 / CaSO4; carbonate: CaCO3 / FeCO3; calcium shared by anhydrite and calcite), giving OLI/ScaleChem-style precipitated amounts plus the residual brine. It reuses the same Ksp/ion-pairing/activity chemistry as ScalePredictionCalculator, and can upgrade the activity model from Davies (I ≤ 0.5 m) to an extended B-dot (Helgeson) model for high-salinity brines:

  MultiMineralScaleEquilibrium eq = new MultiMineralScaleEquilibrium(predictionCalc)
      .setActivityModel(MultiMineralScaleEquilibrium.ActivityModel.BDOT) // high-I option
      .setWaterVolume(1.0);
  eq.solve();
  double baso4MgL = eq.getPrecipitatedMassMgPerL("BaSO4");
  double totalMgL = eq.getTotalScaleMassMgPerL();
  double srResid  = eq.getResidualFreeIonMolPerL("SO4--");
  // each precipitated mineral ends at SI≈0; suppressed minerals stay SI≤0
  

  For high-pressure brines, enable the second-order (compressibility) Ksp term on the predictor first: predictionCalc.setSecondOrderPressureCorrection(true).

• In a process flowsheet (scaling mass RATE): StreamScaleAnalyzer extracts the aqueous ion chemistry + produced-water throughput from a run Stream and turns the coupled equilibrium into a kg/day scaling rate — the number that matters for operations and deposition:

  StreamScaleAnalyzer a = StreamScaleAnalyzer.fromStream(producedWaterStream);
  a.analyze();
  double bariteKgPerDay = a.getScaleRateKgPerDay("BaSO4");
  double totalKgPerDay  = a.getTotalScaleRateKgPerDay();
  String dominant       = a.getDominantScale();
  

  (Needs an electrolyte-CPA fluid with an aqueous ion phase; a non-electrolyte stream carries no ions and reports no scale.)

• Root cause analysis: RootCauseAnalyser.setWaterChemistryFromStream(stream) (or the ion setters) makes a scale-deposit symptom self-quantifying — the MINERAL_SCALE candidate then carries the coupled dominant mineral and its predicted mg/L, and is down-weighted when the brine is thermodynamically undersaturated (deposit is more likely wax/asphaltene/corrosion product).

• Agentic (MCP): the runChemistry tool exposes analysis:"multiMineralScale" (ions in mg/L + T/P + optional waterFlow_LPerDay, activityModel:"BDOT", secondOrderPressure) returning per-mineral amounts and a kg/day rate.

Model validation: the Ksp correlations are pinned to published log10(Ksp) at 25 °C in ScaleKspLiteratureBenchmarkTest — calcite −8.48, barite −9.97, celestite −6.63, anhydrite −4.36, siderite −10.89 (Greenberg & Tomson 1992). Keep the siderite form -59.3498 - 0.041377*T - 2.1963/T + 24.5724*log10(T) consistent across ScalePredictionCalculator, CheckScalePotential and CalcSaltSatauration (no T^2 term).

Kinetics — SI says if, not how fast. ScaleKinetics adds induction time and growth regime on top of the SI:

ScaleKinetics k = new ScaleKinetics().setSaturationIndex(si).evaluate();
double tInd = k.getInductionTimeHours();      // classical nucleation
String regime = k.getLimitingRegime();        // REACTION / TRANSPORT / NONE
double growth = k.getEffectiveGrowthRateMmYr();

Brine mixing (seawater + formation water incompatibility) — sulphate scale often peaks at an intermediate mixing ratio, not either end member:

BrineMixingScaleEvaluator mix = new BrineMixingScaleEvaluator(formationWater, seawater);
mix.setConditions(60, 100).setPH(6.0).setSteps(21).evaluate();
double worstFrac = mix.getWorstFractionA();    // fraction of formation water
String worst = mix.getWorstMineral();          // e.g. BaSO4

Deposition along a line: ScaleDepositionAccumulator (walks a PipeBeggsAndBrills), or the coupled corrosion+scale PipeSegmentIntegrity above.

Valve scale drift (plugging → Cv loss → opening drift → RCA)

Scale/solids on a control-valve trim do not fail the valve suddenly: the deposit shrinks the open flow area, the effective Kv/Cv drops, and the level/pressure controller opens the valve further to hold setpoint — an observable upward drift of the opening until it pins at 100% and control is lost. This is the classic level-valve-plugging signature. ThrottlingValve now carries a fouling term and ValveScaleDrift turns a deposit growth rate into that drift:

valve.setFoulingFraction(0.5);        // effectiveKv = Kv*(1-f); 0=clean
double effCv = valve.getEffectiveCv();

ValveScaleDrift drift = new ValveScaleDrift(valve)
    .setPortDiameter(0.05)            // trim port diameter [m]
    .setKinetics(scaleKinetics);      // or setGrowthRateMmPerYear(20.0)
for (int day = 0; day < 60; day++) {
  drift.advance(1.0);                 // grow deposit, update valve foulingFraction
  process.run();                      // controller compensates -> opening drifts up
}
double daysToPin = drift.predictTimeToPinDays(cleanOpeningPercent); // loss of control
double tPlug     = drift.getTimeToPlugDays();                       // full closure

Uniform radial deposit model: foulingFraction = 1 - ((d0 - 2t)/d0)^2; the opening required to hold flow rises as cleanOpening / (1 - fouling). Use it in an RCA to reproduce the "both LVs → 100% open, level still rising, no inflow surge" trend and to bound time-to-plug.

Scale / precipitation remediation (dissolvers & washing for RCA solutions)

Once a deposit is identified, ScaleRemediationAdvisor recommends the dissolver / solvent / wash to clean the already-fouled equipment (valves, trim, bridles, lines, exchangers) — the proposed-solution counterpart to inhibitor prevention. It is backed by a knowledge-base CSV (/data/scale_remediation.csv) and maps common mineral names to canonical formulae:

ScaleRemediationAdvisor advisor = new ScaleRemediationAdvisor();
List<ScaleRemediationAdvisor.RemediationOption> opts = advisor.recommendFor("calcite");
ScaleRemediationAdvisor.RemediationOption best = opts.get(0); // ordered by effectiveness
best.getDissolver();   // e.g. "Hydrochloric acid (HCl) with corrosion inhibitor ..."
best.getMethod();      // "Acid soak or circulation, then neutralise and flush"
best.getCautions();    // safety/operational cautions
String json = advisor.toJson("BaSO4");

Coverage & the key gotcha it encodes: acid dissolves carbonate/sulfide scale (CaCO3, FeCO3, FeS) but NOT sulfate scale — barite/celestite (BaSO4, SrSO4) need a high-pH chelant (DTPA/EDTA); dithiazine (H2S-scavenger solids) needs a proprietary dissolver plus restoring water pH control to prevent recurrence; asphaltene/wax/naphthenate use aromatic solvent / hot oil. Unknown deposits return a "sample-and-identify first (XRD/SEM-EDX)" guard.

RootCauseAnalyser uses this automatically: every deposit candidate (MINERAL_SCALE, WAX_DEPOSITION, ASPHALTENE, FES_DEPOSITION) now appends a concrete "to clean fouled equipment: …" hint to its recommendation, targeting the coupled dominant mineral when ion chemistry is available.

Deposits on compressors (fouling → performance loss → washing)

When deposit-forming species (elemental sulfur S8, salt from entrained produced water, mineral scale, wax) reach a compressor, they foul the impeller and reduce head/efficiency. The neqsim.process.equipment.compressor deposit model bridges these flow-assurance calculations to compressor performance and washing:

CompressorDeposit dep = CompressorDeposit.fromCompressor(comp);
dep.accumulate(new SolidFlashDepositSource(feed, "S8", DepositMechanism.SULFUR_S8, 0.3), 500.0);
dep.accumulate(new EntrainedSaltDepositSource(10.0, 0.05), 500.0); // entrained brine -> salt
comp.setDepositModel(dep);            // run() degrades efficiency/power
CompressorChart chart500 = comp.buildDegradedChart();               // map after 500 h
WashFluid wash = CompressorDepositWash.recommend(dep);              // salt->WATER, S8->XYLENE
comp.washOnline(wash, 300.0, 3.0);    // online wash removes matching deposits

Per-impeller deposit location: CompressorDepositProfile.compute(...) / computeFromPropertyProfile(...). Full API in the neqsim-api-patterns skill and the compressor_deposit_degradation doc. Use this to plan sulfur/salt washing of fouled compressors and to recommend a wash fluid.

6. Thermal Analysis

Cooldown Calculation

// Simple cooldown time estimate for insulated pipeline
// Use Newton's cooling law: T(t) = Tsea + (Tin - Tsea) * exp(-t/tau)
// where tau = m * Cp / (U * A)
// U = overall heat transfer coefficient (W/m2K)
// A = pipe surface area per unit length (m2/m)
// m = fluid mass per unit length (kg/m)
// Cp = fluid heat capacity (J/kgK)

Arrival Temperature Check

// Critical checks:
// 1. Arrival T > hydrate formation T + margin (typically 5°C subcooling)
// 2. Arrival T > WAT (if waxy crude) + margin
// 3. Arrival T > pour point (for restart)
// If not met: increase insulation, add DEH, reduce flow, or inject inhibitor

7. Flow Assurance Decision Matrix

8. CO2 Injection Well Analysis

For CCS/injection well safety analysis, use the dedicated analyzer:

CO2InjectionWellAnalyzer analyzer = new CO2InjectionWellAnalyzer("InjWell-1");
analyzer.setFluid(co2Fluid);
analyzer.setWellGeometry(measuredDepth, innerDiameter, roughness);
analyzer.setOperatingConditions(inletP_bara, inletT_C, massFlowRate_kg_h);
analyzer.setFormationTemperature(surfaceT_C, bottomholeT_C);
analyzer.addTrackedComponent("hydrogen", maxMolFracLimit);
analyzer.runFullAnalysis();
boolean safe = analyzer.isSafeToOperate();

Impurity Monitoring

ImpurityMonitor monitor = new ImpurityMonitor("H2-Mon", stream);
monitor.addTrackedComponent("hydrogen", 0.10);  // 10 mol% limit
monitor.addTrackedComponent("H2S", 0.001);      // 0.1% limit
double enrichment = monitor.getEnrichmentFactor("hydrogen");
boolean exceeds = monitor.exceedsLimit("hydrogen");

9. Elemental Sulfur (S8) Deposition

Elemental sulfur (cyclo-octasulfur, S8) drops out of natural gas as a yellow/grey solid at points of pressure and/or temperature letdown — compressor inlets, control and letdown valves, dry-gas seals, pressure regulators, and filters. The gas dissolves a small amount of S8 (solubility rises with pressure and with H2S/CO2 content, and is high in condensate, MEG, methanol and TEG); when pressure or temperature drops, the gas becomes supersaturated and S8 deposits. The root cause is almost always an oxygen source oxidising H2S: 8 H2S + 4 O2 -> S8 + 8 H2O (O2 typically enters via preservation fluids, platform-nitrogen purge/seal gas, or injected chemicals).

Screening workflow: assume the gas is S8-saturated at a baseline (e.g. scrubber / separator conditions such as 100 bara, 45 C), then check whether S8 drops to a solid as the stream follows the compressor / valve pressure-temperature path. S8 is a database component (fluid.addComponent("S8", ...)), and deposition is detected with a solid flash.

S8 solid-drop-out (solubility) check

// Gas saturated with a trace of S8 at baseline, then evaluate a P/T letdown point
SystemInterface fluid = new SystemSrkEos(273.15 + 45.0, 100.0);
fluid.addComponent("methane", 0.90);
fluid.addComponent("CO2", 0.02);
fluid.addComponent("H2S", 0.001);
fluid.addComponent("S8", 1.0e-6);   // trace S8 (mole basis)
fluid.setMixingRule("classic");
fluid.setMultiPhaseCheck(true);
fluid.setSolidPhaseCheck("S8");      // enable S8 solid-phase modelling

// Evaluate at the letdown condition (e.g. valve / compressor-inlet P and T)
fluid.setPressure(60.0, "bara");
fluid.setTemperature(30.0, "C");
ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.TPSolidflash();                 // NOT TPflash — needed to form the solid S8 phase
fluid.initProperties();

boolean s8Deposits = fluid.hasPhaseType(PhaseType.SOLID);
if (s8Deposits) {
  double solidMass = fluid.getPhaseOfType("solid").getMass("kg");  // S8 solid inventory
}

Rule of thumb: deposition risk exists wherever a saturated (or nearly saturated) stream sees a pressure or temperature drop. Heating the gas before letdown or reducing the dP moves the point back into the single-phase (dissolved) region.

SulfurFilter equipment (solid S8 removal + change interval)

SulfurFilter runs a TPSolidflash internally, removes the solid S8, and tracks the kg/hr loading for element sizing and change-out interval:

SulfurFilter filter = new SulfurFilter("S8 Filter", valveOutletStream);
filter.setRemovalEfficiency(0.99);       // 99% solid removal
filter.setFilterElementCapacity(50.0);   // kg S8 per element set
filter.setDeltaP(0.5);                    // clean pressure drop [bar]
filter.run();

boolean s8Present = filter.isSolidS8Detected();
double s8Rate = filter.getSolidSulfurRemovalRate();  // kg/hr captured
double interval = filter.getChangeIntervalHours();   // hours until element change

Particle nucleation (optional, for particle-size / filter rating)

ClassicalNucleationTheory cnt = ClassicalNucleationTheory.sulfurS8();
// pair with PopulationBalanceModel for particle-size distribution across the letdown

Mitigations (in order of effectiveness): eliminate the O2 source (nitrogen/ preservation-fluid quality and routing); heat the gas upstream of pressure reduction; reduce the dP / minimum landing pressure; install a SulfurFilter (watch the dP); keep dry-gas-seal gas warm and well above its S8 saturation point.

Related: SourServiceAssessment (setElementalSulfurPresent(true)) for materials impact; neqsim-electrolyte-systems for the aqueous-phase corrosivity of wet sulfur.

10. Common Pitfalls