Phase Envelope and Critical Points Guide

The phase envelope defines the boundary between single-phase and two-phase regions. This guide explains how to calculate phase envelopes and extract key points like cricondenbar and cricondentherm using NeqSim.

Key Concepts

Term	Definition	Importance
Bubble Point	Pressure/temperature where first vapor bubble forms	Oil systems, pipeline pressure
Dew Point	Pressure/temperature where first liquid droplet forms	Gas systems, condensation
Cricondenbar	Maximum pressure on phase envelope	Gas processing design
Cricondentherm	Maximum temperature on phase envelope	Pipeline operations
Critical Point	Where bubble and dew point curves meet	Phase behavior understanding
Quality Lines	Constant liquid fraction within two-phase region	Flash calculations

Algorithm Overview

NeqSim uses a continuation method (Newton-Raphson based) to trace the phase envelope. Understanding the algorithm helps interpret results and troubleshoot issues.

How the Algorithm Works

Starting Point Selection: Based on the bubblePointFirst parameter:
- bubblePointFirst=true: Starts at low pressure, finds bubble point using the most volatile component (lowest critical temperature)
- bubblePointFirst=false (default): Starts at low pressure, finds dew point using the least volatile component (highest critical temperature)
Initial Temperature Estimation: Uses Wilson K-value correlation to estimate starting temperature at the specified low pressure
Envelope Tracing: Newton-Raphson solver steps along the envelope, adjusting T and P to maintain phase equilibrium
Critical Point Detection: When K-values for both light and heavy components approach 1.0 (K ≈ 1), the algorithm:
- Records the critical point position
- Inverts phase types (switches from bubble to dew branch or vice versa)
- Continues tracing on the other side
Automatic Restart: If the algorithm fails before completing, it automatically restarts from the opposite end to capture the full envelope
Cricondenbar/Cricondentherm Refinement: After tracing, the discrete maximum-pressure and maximum-temperature estimates are refined using the Michelsen simultaneous Newton method (see Direct Cricondenbar and Cricondentherm Refinement)

Data Structure

The Michelsen continuation method traces the envelope as two logical branches separated by the critical point:

Branch	Keys	Description
Dew curve	`dewT`, `dewP`	Temperatures/pressures along the dew point curve
Bubble curve	`bubT`, `bubP`	Temperatures/pressures along the bubble point curve
Critical	`criticalPoint1`	`[T(K), P(bara)]` where K-values converge to 1
Cricondenbar	`cricondenbar`	`[T(K), P(bara)]` at maximum pressure
Cricondentherm	`cricondentherm`	`[T(K), P(bara)]` at maximum temperature

Quick Phase Envelope Calculation

Java

import neqsim.thermo.system.SystemSrkEos;
import neqsim.thermodynamicoperations.ThermodynamicOperations;

// Create natural gas
SystemSrkEos fluid = new SystemSrkEos(273.15, 50.0);
fluid.addComponent("nitrogen", 0.02);
fluid.addComponent("CO2", 0.03);
fluid.addComponent("methane", 0.80);
fluid.addComponent("ethane", 0.08);
fluid.addComponent("propane", 0.04);
fluid.addComponent("i-butane", 0.015);
fluid.addComponent("n-butane", 0.015);
fluid.setMixingRule("classic");

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);

// Calculate phase envelope (default: starts from dew point)
ops.calcPTphaseEnvelope();

// Get cricondenbar [T(K), P(bara)]
double[] cricondenbar = ops.get("cricondenbar");
double cricondenbarP = cricondenbar[1];    // Pressure (bara)
double cricondenbarT = cricondenbar[0];    // Temperature at cricondenbar (K)

// Get cricondentherm [T(K), P(bara)]
double[] cricondentherm = ops.get("cricondentherm");
double criconderthermT = cricondentherm[0];  // Temperature (K)
double criconderthermP = cricondentherm[1];  // Pressure at cricondentherm (bara)

// Get critical point [T(K), P(bara)]
double[] criticalPoint = ops.get("criticalPoint1");
double criticalT = criticalPoint[0];
double criticalP = criticalPoint[1];

System.out.println("=== Phase Envelope Results ===");
System.out.printf("Cricondenbar: %.2f bara at %.1f °C%n",
    cricondenbarP, cricondenbarT - 273.15);
System.out.printf("Cricondentherm: %.1f °C at %.2f bara%n",
    criconderthermT - 273.15, criconderthermP);
System.out.printf("Critical point: %.1f °C, %.2f bara%n",
    criticalT - 273.15, criticalP);

// Get all envelope data for plotting (2D array)
double[][] envelope = ops.getData();
// envelope contains multiple branches - see "Retrieving Envelope Data" section

Python

from neqsim import jneqsim

SystemSrkEos = jneqsim.thermo.system.SystemSrkEos
ThermodynamicOperations = jneqsim.thermodynamicoperations.ThermodynamicOperations

# Create gas
fluid = SystemSrkEos(273.15 + 20.0, 50.0)
fluid.addComponent("nitrogen", 0.02)
fluid.addComponent("CO2", 0.03)
fluid.addComponent("methane", 0.80)
fluid.addComponent("ethane", 0.08)
fluid.addComponent("propane", 0.04)
fluid.addComponent("n-butane", 0.03)
fluid.setMixingRule("classic")

ops = ThermodynamicOperations(fluid)

# Calculate phase envelope
ops.calcPTphaseEnvelope()

# Get critical points using the get() method
cricondenbar = list(ops.get("cricondenbar"))      # [T(K), P(bara)]
cricondentherm = list(ops.get("cricondentherm"))  # [T(K), P(bara)]
critical = list(ops.get("criticalPoint1"))        # [T(K), P(bara)]

print(f"Cricondenbar: {cricondenbar[1]:.2f} bara at {cricondenbar[0] - 273.15:.1f} °C")
print(f"Cricondentherm: {cricondentherm[0] - 273.15:.1f} °C at {cricondentherm[1]:.2f} bara")
print(f"Critical Point: {critical[0] - 273.15:.1f} °C, {critical[1]:.2f} bara")

Controlling the Starting Point

You can control whether the algorithm starts from the bubble point or dew point side:

// Start from BUBBLE point (trace bubble curve first, then dew after critical)
ops.calcPTphaseEnvelope(true);  // bubblePointFirst = true

// Start from DEW point (default behavior)
ops.calcPTphaseEnvelope(false); // bubblePointFirst = false

// With custom low pressure starting point
ops.calcPTphaseEnvelope(true, 0.5);  // bubbleFirst=true, lowPres=0.5 bar

// With custom phase fraction (0 = bubble, 1 = dew)
ops.calcPTphaseEnvelope(1.0, 1e-10);  // lowPres=1 bar, phaseFraction≈0 (bubble)

When to use each:

Bubble first: Better for oil-rich systems where bubble point is more important
Dew first (default): Better for gas systems where dew point/HCDP is critical

Retrieving Envelope Data

Available Data Keys

Use ops.get(key) to retrieve specific envelope data:

Key	Returns	Description
`"dewT"`	`double[]`	Dew point temperatures (K) - first branch
`"dewP"`	`double[]`	Dew point pressures (bara) - first branch
`"bubT"`	`double[]`	Bubble point temperatures (K) - second branch
`"bubP"`	`double[]`	Bubble point pressures (bara) - second branch
`"dewT2"`	`null`	Legacy key, returns `null` (Michelsen merges all dew data into `dewT`)
`"dewP2"`	`null`	Legacy key, returns `null` (Michelsen merges all dew data into `dewP`)
`"bubT2"`	`null`	Legacy key, returns `null` (Michelsen merges all bubble data into `bubT`)
`"bubP2"`	`null`	Legacy key, returns `null` (Michelsen merges all bubble data into `bubP`)
`"cricondenbar"`	`double[]`	[T(K), P(bara)] at maximum pressure
`"cricondentherm"`	`double[]`	[T(K), P(bara)] at maximum temperature
`"criticalPoint1"`	`double[]`	[Tc(K), Pc(bara)] critical point
`"dewH"`	`double[]`	Dew point enthalpies
`"bubH"`	`double[]`	Bubble point enthalpies
`"dewDens"`	`double[]`	Dew point densities
`"bubDens"`	`double[]`	Bubble point densities
`"dewS"`	`double[]`	Dew point entropies
`"bubS"`	`double[]`	Bubble point entropies

Getting All Data

// Get all envelope points as 2D array
double[][] allData = ops.getData();

// allData structure depends on whether envelope completed or restarted:
// allData[0], allData[1] = First branch (T, P)
// allData[2], allData[3] = Second branch (T, P) after critical point
// allData[4], allData[5] = Third branch if restarted
// allData[6], allData[7] = Fourth branch if restarted

Complete Example: Extracting Both Curves

// Calculate envelope
ops.calcPTphaseEnvelope();

// Get temperatures and pressures for plotting
double[] dewT = ops.get("dewT");   // First branch temperatures
double[] dewP = ops.get("dewP");   // First branch pressures
double[] bubT = ops.get("bubT");   // Second branch temperatures (after critical)
double[] bubP = ops.get("bubP");   // Second branch pressures

// Note: "dew" and "bub" naming depends on bubblePointFirst setting
// With default (bubblePointFirst=false):
//   dewT/dewP = dew point curve (before critical)
//   bubT/bubP = bubble point curve (after critical)

// Critical point
double[] crit = ops.get("criticalPoint1");
System.out.printf("Critical: %.1f K, %.2f bar%n", crit[0], crit[1]);

// Combine for full envelope
System.out.println("Dew curve points: " + dewT.length);
System.out.println("Bubble curve points: " + bubT.length);

Plotting Phase Envelopes

With Matplotlib (Python)

import matplotlib.pyplot as plt
import numpy as np
from neqsim import jneqsim

# Create fluid
SystemSrkEos = jneqsim.thermo.system.SystemSrkEos
ThermodynamicOperations = jneqsim.thermodynamicoperations.ThermodynamicOperations

fluid = SystemSrkEos(273.15 + 20.0, 50.0)
fluid.addComponent("methane", 0.80)
fluid.addComponent("ethane", 0.10)
fluid.addComponent("propane", 0.05)
fluid.addComponent("n-butane", 0.05)
fluid.setMixingRule("classic")

ops = ThermodynamicOperations(fluid)
ops.calcPTphaseEnvelope()

# Get envelope data
dewT = np.array(list(ops.get("dewT"))) - 273.15  # Convert to °C
dewP = np.array(list(ops.get("dewP")))
bubT = np.array(list(ops.get("bubT"))) - 273.15
bubP = np.array(list(ops.get("bubP")))

# Get critical points
cricondenbar = list(ops.get("cricondenbar"))
cricondentherm = list(ops.get("cricondentherm"))
critical = list(ops.get("criticalPoint1"))

# Plot
plt.figure(figsize=(10, 7))
plt.plot(dewT, dewP, 'b-', label='Dew Point', linewidth=2)
plt.plot(bubT, bubP, 'r-', label='Bubble Point', linewidth=2)

# Mark critical points
plt.plot(cricondenbar[0] - 273.15, cricondenbar[1], 'ko', markersize=10, label='Cricondenbar')
plt.plot(cricondentherm[0] - 273.15, cricondentherm[1], 'g^', markersize=10, label='Cricondentherm')
plt.plot(critical[0] - 273.15, critical[1], 'rs', markersize=10, label='Critical Point')

plt.xlabel('Temperature (°C)', fontsize=12)
plt.ylabel('Pressure (bara)', fontsize=12)
plt.title('Phase Envelope', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

With JFreeChart (Java - built-in)

// NeqSim has built-in JFreeChart display
ops.calcPTphaseEnvelope();
ops.displayResult();  // Opens a window with the phase envelope plot

Hydrocarbon Dew Point (HCDP)

For gas pipeline specifications, the hydrocarbon dew point is critical:

// Calculate hydrocarbon dew point at specific pressure
double pipelinePressure = 70.0;  // bara
fluid.setPressure(pipelinePressure);

ops.dewPointTemperatureFlash();
double hcdp = fluid.getTemperature("C");

System.out.printf("HCDP at %.0f bara: %.1f °C%n", pipelinePressure, hcdp);

// Check against spec (e.g., -2°C at delivery pressure)
double maxHCDP = -2.0;
if (hcdp > maxHCDP) {
    System.out.printf("⚠️ HCDP %.1f °C exceeds spec %.1f °C%n", hcdp, maxHCDP);
    System.out.println("   Need NGL extraction or C3+ removal");
} else {
    System.out.printf("✅ HCDP %.1f °C meets spec%n", hcdp);
}

Water Dew Point

For pipeline water specifications:

// Add water to gas
fluid.addComponent("water", 0.001);  // Small amount

// Calculate water dew point at pipeline pressure
fluid.setPressure(70.0);
ops.waterDewPointTemperatureFlash();
double waterDP = fluid.getTemperature("C");

System.out.printf("Water dew point at 70 bara: %.1f °C%n", waterDP);

// Pipeline spec typically requires water DP below minimum ground temp
double minGroundTemp = 5.0;
if (waterDP > minGroundTemp) {
    System.out.println("⚠️ Need dehydration (TEG)");
}

Quality Lines (Iso-Vapor Fraction)

Quality lines show constant vapor/liquid fraction within the two-phase region. NeqSim calculates these by running phase envelope calculations at different phase fractions:

// Calculate quality lines by varying phaseFraction parameter
// phaseFraction = 0 -> bubble point (0% vapor)
// phaseFraction = 1 -> dew point (100% vapor)

double[] qualities = {0.1, 0.25, 0.5, 0.75, 0.9};  // vapor fractions

for (double q : qualities) {
    ThermodynamicOperations opsQ = new ThermodynamicOperations(fluid.clone());
    // lowPres=1.0, phaseFraction=q
    opsQ.calcPTphaseEnvelope(1.0, q);

    double[][] qualityLine = opsQ.getData();
    System.out.printf("Quality %.0f%% vapor - %d points%n", q * 100, qualityLine[0].length);
}

Note: The calcPTphaseEnvelope(double lowPres, double phaseFraction) method allows specifying a phase fraction between 0 (bubble) and 1 (dew) to trace quality lines.

Phase Envelope for Different Fluid Types

Lean Gas (Dry Gas)

SystemSrkEos leanGas = new SystemSrkEos(273.15, 50.0);
leanGas.addComponent("methane", 0.95);
leanGas.addComponent("ethane", 0.03);
leanGas.addComponent("propane", 0.01);
leanGas.addComponent("nitrogen", 0.01);
leanGas.setMixingRule("classic");

// Lean gas has small phase envelope, low cricondentherm
// Typically single-phase at pipeline conditions

Rich Gas (Wet Gas / Condensate)

SystemSrkEos richGas = new SystemSrkEos(273.15, 50.0);
richGas.addComponent("methane", 0.70);
richGas.addComponent("ethane", 0.10);
richGas.addComponent("propane", 0.08);
richGas.addComponent("i-butane", 0.03);
richGas.addComponent("n-butane", 0.04);
richGas.addComponent("n-pentane", 0.02);
richGas.addTBPfraction("C6+", 0.03, 90.0, 0.70);
richGas.setMixingRule("classic");

// Rich gas has larger phase envelope
// May condense liquids in pipeline - flow assurance concern

Black Oil

SystemSrkEos oil = new SystemSrkEos(373.15, 200.0);
oil.addComponent("methane", 0.30);
oil.addComponent("ethane", 0.05);
oil.addComponent("propane", 0.04);
oil.addComponent("n-butane", 0.03);
oil.addComponent("n-pentane", 0.03);
oil.addTBPfraction("C7", 0.10, 100.0, 0.74);
oil.addTBPfraction("C15", 0.20, 210.0, 0.82);
oil.addTBPfraction("C30+", 0.25, 420.0, 0.90);
oil.setMixingRule("classic");

// Oil systems: bubble point line is most important
// Critical point at high temperature (often > 400°C)

Bubble Point Pressure

For oils, the bubble point at reservoir temperature is key:

// Set reservoir temperature
double reservoirT = 100.0;  // °C
fluid.setTemperature(reservoirT + 273.15);

// Calculate bubble point pressure
ops.bubblePointPressureFlash();
double bubblePointP = fluid.getPressure("bara");

System.out.printf("Bubble point at %.0f°C: %.1f bara%n", reservoirT, bubblePointP);

// Compare to reservoir pressure for undersaturation check
double reservoirP = 300.0;
double undersaturation = reservoirP - bubblePointP;
System.out.printf("Undersaturation: %.1f bar%n", undersaturation);

Effect of Composition on Phase Envelope

Adding Heavy Components

// Base case
SystemSrkEos gas1 = createBaseGas();
ThermodynamicOperations ops1 = new ThermodynamicOperations(gas1);
ops1.calcPTphaseEnvelope();
double[] cricoT1 = ops1.get("cricondentherm");
double cricondentherm1 = cricoT1[0] - 273.15;

// With added C6+
SystemSrkEos gas2 = createBaseGas();
gas2.addComponent("n-hexane", 0.02);
ThermodynamicOperations ops2 = new ThermodynamicOperations(gas2);
ops2.calcPTphaseEnvelope();
double[] cricoT2 = ops2.get("cricondentherm");
double cricondentherm2 = cricoT2[0] - 273.15;

System.out.printf("Cricondentherm without C6+: %.1f °C%n", cricondentherm1);
System.out.printf("Cricondentherm with C6+: %.1f °C%n", cricondentherm2);
System.out.println("Adding heavy components increases cricondentherm!");

Adding CO2

// CO2 typically reduces cricondentherm but increases cricondenbar
// Important for CCS applications

Pipeline Operating Point Analysis

Check if pipeline conditions are in single-phase region:

// Get cricondentherm from envelope calculation
double[] cricoTherm = ops.get("cricondentherm");
double cricondenthermC = cricoTherm[0] - 273.15;

// Pipeline operating conditions
double[] pipelineT = {10, 5, 0, -5};  // °C along pipeline
double pipelineP = 70.0;  // bara (constant pressure assumption)

System.out.println("Pipeline Operating Point Analysis:");
System.out.printf("Cricondentherm: %.1f °C%n", cricondenthermC);
System.out.printf("Operating pressure: %.0f bara%n%n", pipelineP);

for (double t : pipelineT) {
    fluid.setTemperature(t + 273.15);
    fluid.setPressure(pipelineP);
    ops.TPflash();

    int nPhases = fluid.getNumberOfPhases();
    String status = (nPhases == 1) ? "Single phase ✅" : "Two phase ⚠️";
    double liquidFrac = (nPhases > 1) ? fluid.getBeta(1) : 0.0;

    System.out.printf("T=%3.0f°C: %s (liquid frac: %.3f)%n", t, status, liquidFrac);
}

Troubleshooting Phase Envelope Calculations

Issue	Possible Cause	Solution
Envelope doesn’t close	Algorithm failed before reaching critical point	Try starting from the other side with `calcPTphaseEnvelope(true)`
Missing branch	Algorithm failed before completing	Try starting from the other side with `calcPTphaseEnvelope(true)` or adjust starting pressure
Cricondentherm too high	Heavy components or characterization	Review plus fraction properties
No dew point found	Fluid too light (dry gas)	Normal for methane-rich gas
Calculation fails	Near-critical region or bad initial guess	Try different starting pressure with `calcPTphaseEnvelope(true, 0.5)`
NaN values in results	Non-convergence	Check fluid composition, try simpler EoS

Understanding Algorithm Behavior

The Michelsen continuation method traces the envelope through the critical point, producing two branches:

ops.calcPTphaseEnvelope();

// Check what data is available
double[] dewT = ops.get("dewT");
double[] bubT = ops.get("bubT");
double[] cricondenbar = ops.get("cricondenbar");
double[] criticalPoint = ops.get("criticalPoint1");

System.out.println("Dew curve: " + (dewT != null ? dewT.length : 0) + " points");
System.out.println("Bubble curve: " + (bubT != null ? bubT.length : 0) + " points");
System.out.println("Cricondenbar: " + cricondenbar[1] + " bara");
System.out.println("Critical point: " + criticalPoint[0] + " K, " + criticalPoint[1] + " bara");

Quality Lines

For constant vapor fraction curves inside the two-phase region:

// Calculate envelope with quality lines at specified vapor fractions
double[] betaValues = {0.1, 0.25, 0.5, 0.75, 0.9};
ops.calcPTphaseEnvelopeWithQualityLines(betaValues);

// Access quality line data
double[] qualT = ops.get("qualityT_0.5");      // Temperatures at 50% vapor
double[] qualP = ops.get("qualityP_0.5");      // Pressures at 50% vapor
double[] volFrac = ops.get("qualityVolFrac_0.5");   // Volume fractions
double[] massFrac = ops.get("qualityMassFrac_0.5"); // Mass fractions

The phase envelope tracer finds the cricondenbar and cricondentherm as discrete maximum-pressure and maximum-temperature points on the traced curve. For high-precision applications, these estimates can be refined to machine precision using thethe Michelsen simultaneous Newton method.

How It Works

Instead of the envelope tracer’s discrete tracking (which is limited by the step size between traced points), the refinement solves a full $(n_c + 2)$-dimensional Newton system simultaneously:

Unknowns: $\ln K_1, \dots, \ln K_{n_c}, \ln T, \ln P$
Equations: $n_c$ equilibrium equations (fugacity equality), 1 Rachford-Rice summation, and 1 extremum specification

The extremum specification differs for each type:

Extremum	Condition	Specification Equation
Cricondenbar	$dT/dP = 0$	$S_P = \sum_i s_i \cdot T(\partial \ln \hat{\varphi}_i^V/\partial T - \partial \ln \hat{\varphi}_i^L/\partial T) = 0$
Cricondentherm	$dP/dT = 0$	$S_T = \sum_i s_i \cdot P(\partial \ln \hat{\varphi}_i^V/\partial P - \partial \ln \hat{\varphi}_i^L/\partial P) = 0$

where $s_i = (y_i - x_i) / |\mathbf{y} - \mathbf{x}|$ is a normalized composition direction vector.

Convergence

The simultaneous Newton system provides quadratic convergence, typically reaching $|\mathbf{g}|_2 < 10^{-10}$ in 3-8 iterations from the discrete envelope estimate. Key features:

Analytical Jacobian using fugacity coefficient derivatives from system.init(3)
Damped Newton steps to handle large initial corrections
Gaussian elimination with partial pivoting (no external matrix library)
Graceful fallback to the discrete estimate if Newton does not converge

Usage (Java)

The refinement is invoked through calcCricoP() and calcCricoT():

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.calcPTphaseEnvelope();

// Get refined cricondenbar [T(K), P(bara)]
double[] cricondenbar = ops.get("cricondenbar");
System.out.printf("Cricondenbar: %.4f bara at %.2f K%n", cricondenbar[1], cricondenbar[0]);

// Get refined cricondentherm [T(K), P(bara)]
double[] cricondentherm = ops.get("cricondentherm");
System.out.printf("Cricondentherm: %.2f K at %.4f bara%n", cricondentherm[0], cricondentherm[1]);

Implementation Classes

Class	Purpose
`CricondenBarFlash`	Michelsen $(n_c+2)$ Newton refinement for maximum pressure ($S_P = 0$)
`CricondenThermFlash`	Michelsen $(n_c+2)$ Newton refinement for maximum temperature ($S_T = 0$)

For the full mathematical formulation, see the Phase Envelope Algorithm documentation.

API Reference Summary

calcPTphaseEnvelope Overloads

Method	Description
`calcPTphaseEnvelope()`	Default: dew first, lowPres=1.0 bar
`calcPTphaseEnvelope(boolean bubbleFirst)`	Control start side
`calcPTphaseEnvelope(double lowPres)`	Custom starting pressure
`calcPTphaseEnvelope(boolean bubbleFirst, double lowPres)`	Both options
`calcPTphaseEnvelope(double lowPres, double phaseFraction)`	Quality lines (0=bubble, 1=dew)
`calcPTphaseEnvelopeWithQualityLines(double[] betaValues)`	Envelope + quality lines at specified vapor fractions
`calcPTphaseEnvelopeWithQualityLines(boolean bubFirst, double[] betaValues)`	With start side control
`calcPTphaseEnvelope2()`	Deprecated — delegates to `calcPTphaseEnvelope()`
`calcPTphaseEnvelopeNew()`	Deprecated — delegates to `calcPTphaseEnvelope()`

Method	Description
`calcCricoP(...)`	Refine cricondenbar using Michelsen simultaneous Newton ($S_P = 0$)
`calcCricoT(...)`	Refine cricondentherm using Michelsen simultaneous Newton ($S_T = 0$)