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
- Constructors use kelvin for temperature and bara for pressure.
- Prefer unit-qualified setters such as
setPressure(30.0, "bara"). - Add all components and set the mixing rule before constructing the operations facade.
- Call
initProperties()after an equilibrium calculation before reading transport properties such as viscosity or thermal conductivity. - Extensive properties such as enthalpy and entropy depend on the amount of fluid in the system. Preserve the same system amount when using them as PH or PS flash targets.
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:
- Check that temperature, pressure, phase count, and phase types are plausible.
- Check that phase compositions and phase fractions close to one.
- For PH, PS, PU, or VU flashes, compare the solved extensive-property residual with a tolerance appropriate to the system amount and application.
- Test nearby states to detect branch switching or critical-region sensitivity.
- 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.