Skip to the content.

ThermodynamicOperations is the public facade for equilibrium calculations on a NeqSim SystemInterface. It changes the attached fluid in place: set the known state variables on the fluid, call an operation, and then read the solved state from that same fluid.

Units and initialization

Verified state-flash example

The following Java 8 program exercises TP, PH, and PS flashes. The residual checks demonstrate the constraints imposed by the PH and PS operations.

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

public final class StateFlashExample {
  private StateFlashExample() {}

  public static void main(String[] args) {
    SystemInterface fluid = new SystemSrkEos(298.15, 50.0);
    fluid.addComponent("methane", 0.90);
    fluid.addComponent("ethane", 0.07);
    fluid.addComponent("propane", 0.03);
    fluid.setMixingRule("classic");

    ThermodynamicOperations operations = new ThermodynamicOperations(fluid);
    operations.TPflash();
    fluid.initProperties();

    double initialEnthalpy = fluid.getEnthalpy();
    fluid.setPressure(30.0, "bara");
    operations.PHflash(initialEnthalpy);
    if (Math.abs(fluid.getEnthalpy() - initialEnthalpy) > 1.0e-3) {
      throw new IllegalStateException("PH flash did not preserve enthalpy");
    }

    double initialEntropy = fluid.getEntropy();
    fluid.setPressure(70.0, "bara");
    operations.PSflash(initialEntropy);
    if (Math.abs(fluid.getEntropy() - initialEntropy) > 1.0e-6) {
      throw new IllegalStateException("PS flash did not preserve entropy");
    }
  }
}

Common state flashes exposed by the facade include:

Method Fixed variables Typical use
TPflash() Temperature, pressure Equilibrium state at specified conditions
PHflash(H) Pressure, total enthalpy Throttling and heat-duty calculations
PSflash(S) Pressure, total entropy Isentropic compression or expansion
PUflash(U) Pressure, total internal energy Pressure-constrained energy balance
TVflash(V) Temperature, volume Rigid-volume state calculation
VUflash(V, U) Volume, total internal energy Closed-vessel dynamic update

TVfractionFlash(fraction) constrains the volume fraction of the first, lightest phase. It is not a general molar vapour-fraction specification; inspect the resulting phase types before interpreting the fraction.

PS flash convergence

PSflash(S), the unit-qualified PSflash(S, unit) overload and PSflash2(S) return only when the total entropy residual satisfies abs(Sactual - Sspecified) <= max(1e-7 * n, 1e-9 * abs(Sspecified)) J/K, where n is the system amount in moles. Unit-qualified targets are converted to total J/K before applying this criterion. The amount-scaled absolute tolerance also applies when the target entropy is zero or negative. Pressure is preserved; temperature and pressure must be finite and positive, and phase fractions must be finite, bounded by zero and one, and normalized.

The temperature solver uses bounded Newton steps and a sign-changing bracket across mixture phase boundaries. If progress stalls, it rebuilds the bracket with cold TP flashes before checking the endpoint with a final cold TP flash. The bracket trials each use a fresh TP-flash operation: disabling K-value warm starts alone does not discard the operation’s convergence history. Reusing that history can make entropy evaluations with a trace liquid depend on earlier trials and collapse the bracket away from the smooth equilibrium root. The independent bracket uses secant interpolation for smooth roots and bisects at least every third trial to guarantee progress. Rounded or out-of-bracket interpolants fall back to bisection. This bounds the additional TP work without changing entropy acceptance or the normal warm-started Newton path. The independent endpoint check can reject a metastable homogeneous root even when its entropy residual is small. Recovery then uses a fresh inner TP flash so that its previous iteration history cannot reselect the rejected root. The iteration first seeks the tighter max(1e-8 * n, 1e-10 * abs(Sspecified)) J/K residual. Only when a cold sign-changing bracket narrows to four floating-point temperature spacings may the solver accept the best endpoint within the stated postcondition (ten times the iteration tolerance). This handles inner TP-flash resolution without stopping ordinary Newton iterations prematurely or accepting a larger unresolved entropy discontinuity. Pure-component two-phase states retain the tighter iteration tolerance and use saturation temperature and an entropy-based phase fraction, since temperature alone cannot span latent entropy at a fixed pressure. This covers dense CO2 and CO2-rich mixtures crossing into the two-phase region; it does not add solid CO2 equilibrium.

Non-finite entropy targets or invalid initial temperature, pressure or fluid amount raise IllegalArgumentException. Failure to converge raises IllegalStateException, including the entropy residual and state when available, instead of returning the last iterate as a solution. The fluid is modified in place: after an exception, restore a saved inlet or reinitialize before retrying. These semantics concern the standard EOS PS methods above; specialized GERG2008, Leachman and Vega PS methods have their own implementations.

Regression coverage is in PSFlashEntropyClosureTest: fresh and continuation starts for CO2/nitrogen and CO2/hydrogen, pure-CO2 phase entry, independent cold TP-root comparisons, component inventory and phase checks, entropy units and system amounts, and explicit failure behavior. PSFlashTraceLiquidTest adds the issue #4828 PR78 heavy-end reproducer, an independent TP-only root, the complete 216-point compression and amount sweep, and a gas-only control. It checks the tighter iteration tolerance, independent cold-TP agreement, pressure, phase fractions and component conservation. Deterministic evaluation-count checks also guard smooth trace-liquid and CO2 cold recovery and confirm that ordinary gas compression does not enter it.

Saturation and phase-aware operations

The same facade provides the following public operations:

Purpose Method or setup
Bubble pressure at current temperature bubblePointPressureFlash(false)
Bubble temperature at current pressure bubblePointTemperatureFlash()
Dew pressure at current temperature dewPointPressureFlash()
Dew temperature at current pressure dewPointTemperatureFlash()
Water dew temperature waterDewPointTemperatureFlash()
Water-rich liquid dew temperature including dissolved gas waterDewPointTemperatureFlash(true)
CPA liquid-water, ice and hydrate boundaries together cpaWaterPhaseBoundaryTemperatures(); see the setup and phase-reference guide
Hydrate formation temperature hydrateFormationTemperature()
Hydrate formation pressure hydrateFormationPressure()
Multiliquid equilibrium Enable setMultiPhaseCheck(true), then call TPflash()
Wax or other configured solid Call setSolidPhaseCheck(name), then TPflash()

The mixed-liquid water-dew-point operation holds the feed gas composition fixed and solves fugacity equilibrium for all components in an incipient water-rich liquid. It requires more than 50 mol% water in that liquid and excludes reactions, other liquids, solids and hydrates. Its subzero reference is metastable liquid water. The CPA water-dew-point analyser uses this operation; failed measurements return NaN. Use the full multiphase operation for hydrocarbon liquids or glycol-rich condensates. See the dew-point guide for the distinction from GERG-water’s ISO ice reference and validity range.

The direct water/MEG dew-point and Pitzer hydrate temperature/pressure searches reuse bracket-endpoint fugacity residuals for secant interpolation, with bisection at least every second refinement trial. Non-finite, rounded or out-of-bracket predictions use bisection. The logarithmic fugacity tolerances remain 1e-9 for water/MEG dew points and 1e-8 for Pitzer hydrate onset. Pitzer pressure searches retain the first crossing located by the low-pressure scan. BracketedSaturationPerformanceTest checks the previous equilibrium roots, strict residuals, physical state and bounded refinement work.

Solid selection preserves the fluid multiphase-check setting. Request additional liquid phases explicitly with setMultiPhaseCheck(true). The empirical solid model excludes methane and is not a methane-freezing model; see the solid phase flash limitations.

Solid checking is a fluid configuration used by TPflash(); there is no public TPsolidflash() method on ThermodynamicOperations. Hydrate calculations also require a fluid model and components suitable for hydrate equilibrium. See the thermodynamic model guide before selecting an equation of state.

PT phase envelope

calcPTphaseEnvelope() calculates the pressure-temperature phase boundary. Use the facade’s named result accessor for the extrema; the internal operation interface does not expose getCricondenbar() or getCricondentherm() methods.

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

public final class PhaseEnvelopeExample {
  private PhaseEnvelopeExample() {}

  public static double[] run() {
    SystemInterface fluid = new SystemSrkEos(280.0, 10.0);
    fluid.addComponent("methane", 0.75);
    fluid.addComponent("ethane", 0.12);
    fluid.addComponent("propane", 0.08);
    fluid.addComponent("n-butane", 0.05);
    fluid.setMixingRule("classic");

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

    double[] cricondenbar = operations.get("cricondenbar");
    double[] cricondentherm = operations.get("cricondentherm");
    double[] extrema = {
      cricondenbar[0], cricondenbar[1], cricondentherm[0], cricondentherm[1]
    };
    for (double value : extrema) {
      if (!Double.isFinite(value) || value <= 0.0) {
        throw new IllegalStateException("Phase-envelope extremum is not finite and positive");
      }
    }
    return extrema;
  }

  public static void main(String[] args) {
    run();
  }
}

Each named extremum is returned as [temperature K, pressure bara]. Phase envelopes are iterative and may be sensitive near critical conditions. Validate the returned curve and extrema rather than relying only on successful method completion. The facade does not currently expose a public PH-envelope method.

Reactive TP flash

Use reactiveTPflash() for a system containing components registered in NeqSim’s chemical-reaction data. The example below uses a single vapour phase for the water-gas-shift system.

SystemInterface reactive = new SystemSrkEos(600.0, 1.0);
reactive.addComponent("CO", 0.25);
reactive.addComponent("water", 0.25);
reactive.addComponent("CO2", 0.25);
reactive.addComponent("hydrogen", 0.25);
reactive.setMixingRule("classic");
reactive.setMaxNumberOfPhases(1);
reactive.setNumberOfPhases(1);
reactive.init(0);
reactive.init(1);

ThermodynamicOperations operations = new ThermodynamicOperations(reactive);
operations.reactiveTPflash();

double compositionSum = 0.0;
for (int i = 0; i < reactive.getPhase(0).getNumberOfComponents(); i++) {
  compositionSum += reactive.getPhase(0).getComponent(i).getx();
}
if (Math.abs(compositionSum - 1.0) > 1.0e-10) {
  throw new IllegalStateException("Reactive composition does not close");
}

There are no setChemicalReactions(true) or calcChemicalEquilibrium() methods on these public interfaces. For reaction selection, phase constraints, and reactive PH/PS operations, see the reactive flash guide.

Convergence and result checks

ThermodynamicOperations does not provide generic setMaxIterations(...) or setTolerance(...) methods. Solver controls are operation-specific and are not interchangeable across all flash types. Prefer physical validation of each result:

  1. Check that temperature, pressure, phase count, and phase types are plausible.
  2. Check that phase compositions and phase fractions close to one.
  3. For PH, PS, PU, or VU flashes, compare the solved extensive-property residual with a tolerance appropriate to the system amount and application.
  4. Test nearby states to detect branch switching or critical-region sensitivity.
  5. Treat non-convergence as a model, initialization, or state-specification issue; do not silently accept the last iterate.

The legacy OLGApropertyTableGenerator is not a portable general-purpose export API: it does not expose filename or water-cut setters and its current writer has platform-specific behavior. Do not build new workflows around the obsolete example previously shown on this page.