This guide documents the density correction models available in NeqSim for improving volumetric predictions.

Table of Contents

Overview

Density predictions from cubic equations of state (SRK, PR) often have systematic errors:

NeqSim provides volume translation and correlation-based methods to improve liquid density predictions.

Basic density access:

fluid.init(3);  // Initialize with derivatives
fluid.initPhysicalProperties();

double density = fluid.getPhase(1).getDensity("kg/m3");  // Liquid phase
double molarVolume = fluid.getPhase(1).getMolarVolume();  // m³/mol

Equation of State Density

Direct EoS Calculation

Cubic equations of state calculate compressibility factor $Z$:

\[PV = ZnRT\]

Molar volume is then: \(V_m = \frac{ZRT}{P}\)

Density: \(\rho = \frac{PM_w}{ZRT}\)

Issue with cubic EoS: At the critical point, $Z_c^{SRK} = 0.333$ and $Z_c^{PR} = 0.307$, while real hydrocarbons have $Z_c \approx 0.26$. This causes systematic liquid volume overprediction.

Volume Translation Methods

Peneloux Volume Shift

The Peneloux correction adds a constant shift to the EoS molar volume:

\[V_{corrected} = V_{EoS} - c\]

where $c$ is the volume shift parameter.

Class: Peneloux

Mixture shift: \(c_{mix} = \sum_i x_i c_i\)

Component shift correlation: \(c_i = 0.40768 \frac{RT_{c,i}}{P_{c,i}} \left( 0.29441 - Z_{RA,i} \right)\)

where $Z_{RA}$ is the Rackett compressibility factor (from COMP database).

Setting shift parameters:

// Enable Peneloux correction (default for SRK)
fluid.setDensityModel("Peneloux");

// Or set component-specific shifts
fluid.getPhase(0).getComponent("methane").setVolumeCorrectionConst(0.0);
fluid.getPhase(1).getComponent("n-heptane").setVolumeCorrectionConst(-0.0105);

Advantages:

Limitations:

Component-Specific Corrections

NeqSim stores volume correction constants in the COMP database. For heavy hydrocarbons or polar compounds, these may need tuning.

Accessing correction constants:

// Get current volume correction
double vc = fluid.getPhase(1).getComponent("n-decane").getVolumeCorrectionConst();

// Modify correction
fluid.getPhase(1).getComponent("n-decane").setVolumeCorrectionConst(-0.015);

Temperature-dependent shift (Jhaveri-Youngren): Some systems require temperature-dependent corrections:

\[c(T) = c_0 + c_1 (T - T_{ref})\]

This is implemented in specific component models.

Liquid Density Correlations

COSTALD

The COSTALD (COrreSponding STAtes Liquid Density) method is a generalized corresponding-states correlation for predicting liquid densities of pure compounds and mixtures. It is equivalent to the implementation in commercial simulators such as UniSim/HYSYS, Aspen Plus, and PRO/II.

Class: Costald (in neqsim.physicalproperties.methods.liquidphysicalproperties.density)

Method Overview

COSTALD has two parts:

  1. Saturated liquid volume — Hankinson and Thomson (1979)
  2. Compressed liquid correction — Aalto et al. (1996) modified Tait equation

The method automatically applies the compressed liquid correction when the system pressure exceeds the estimated saturation pressure.

Quick Start (Java)

// 1. Create fluid and run flash
SystemInterface fluid = new SystemSrkEos(293.15, 50.0);
fluid.addComponent("methane", 70.0);
fluid.addComponent("n-hexane", 20.0);
fluid.addComponent("water", 10.0);
fluid.setMixingRule("classic");
fluid.setMultiPhaseCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.TPflash();
fluid.initPhysicalProperties();

// 2. Switch liquid phases to COSTALD
fluid.setLiquidDensityModel("COSTALD");

// 3. Read density (applies to oil, liquid, and aqueous phases)
double oilDensity = fluid.getPhase("oil").getPhysicalProperties().getDensity();
double waterDensity = fluid.getPhase("aqueous").getPhysicalProperties().getDensity();

Quick Start (Python / Jupyter)

from neqsim import jneqsim

fluid = jneqsim.thermo.system.SystemSrkEos(293.15, 50.0)
fluid.addComponent("methane", 70.0)
fluid.addComponent("n-hexane", 20.0)
fluid.addComponent("water", 10.0)
fluid.setMixingRule("classic")
fluid.setMultiPhaseCheck(True)

ops = jneqsim.thermodynamicoperations.ThermodynamicOperations(fluid)
ops.TPflash()
fluid.initPhysicalProperties()

fluid.setLiquidDensityModel("COSTALD")

oil_density = fluid.getPhase("oil").getPhysicalProperties().getDensity()
water_density = fluid.getPhase("aqueous").getPhysicalProperties().getDensity()
print(f"Oil:   {oil_density:.1f} kg/m3")
print(f"Water: {water_density:.1f} kg/m3")

API Methods

Method Level Description
fluid.setLiquidDensityModel("COSTALD") System Applies COSTALD to all liquid/oil/aqueous phases at once
fluid.setLiquidDensityModel("Peneloux") System Switches back to the default Peneloux volume shift
phase.getPhysicalProperties().setDensityModel("Costald") Phase Applies COSTALD to a single phase
component.setCostaldCharacteristicVolume(v) Component Sets an explicit V* value (cm³/mol) for fine-tuning
component.getCostaldCharacteristicVolume() Component Returns the explicit V* if set, otherwise 0

Important: Call setLiquidDensityModel() after initPhysicalProperties(). Once set, the density model persists for subsequent calcDensity() calls until changed.

Saturated Liquid Volume (Hankinson-Thomson 1979)

The saturated liquid molar volume is:

\[V_s = V^* \cdot V_R^{(0)}(T_r) \cdot \left[1 - \omega \cdot V_R^{(\delta)}(T_r)\right]\]

where $T_r = T / T_{c,m}$ is the reduced temperature and $\omega$ is the acentric factor.

The dimensionless volume functions are:

\[V_R^{(0)} = 1 - 1.52816\,\tau^{1/3} + 1.43907\,\tau^{2/3} - 0.81446\,\tau + 0.190454\,\tau^{4/3}\] \[V_R^{(\delta)} = \frac{-0.296123 + 0.386914\,T_r - 0.0427258\,T_r^2 - 0.0480645\,T_r^3}{T_r - 1.00001}\]

where $\tau = 1 - T_r$.

Compressed Liquid Correction (Aalto et al. 1996)

When $P > P_{sat}$, a modified Tait correction shrinks the molar volume:

\[V = V_s^{sat} \cdot \frac{A + e^{(d - T_r)^B} \cdot (P_r - P_r^{sat})}{A + e \cdot (P_r - P_r^{sat})}\]

where:

\[A = a_0 + a_1 T_r + a_2 T_r^3 + a_3 T_r^6 + a_4 / T_r\] \[B = b_0 + b_1 \cdot \omega\]
Constant Value
$a_0$ -170.335
$a_1$ -28.578
$a_2$ 124.809
$a_3$ -55.5393
$a_4$ 130.01
$b_0$ 0.164813
$b_1$ -0.0914427
$c$ $e$ (Euler’s number)
$d$ 1.00588

The pseudocritical pressure is:

\[P_{c,m} = \frac{(0.291 - 0.080\,\omega_m) \cdot R \cdot T_{c,m}}{V^*_m}\]

The saturation pressure uses the Lee-Kesler correlation.

Mixing Rules

For mixtures, pseudocritical properties are calculated using the Hankinson-Thomson (1979) mixing rules (Poling, Table 4-12):

Characteristic volume (quadratic):

\[V^*_m = \frac{1}{4}\left[\sum_i x_i V^*_i + 3\left(\sum_i x_i {V^*_i}^{2/3}\right)\left(\sum_i x_i {V^*_i}^{1/3}\right)\right]\]

Acentric factor (linear):

\[\omega_m = \sum_i x_i \omega_i\]

Pseudocritical temperature:

\[T_{c,m} = \frac{\left[\sum_i x_i \sqrt{T_{c,i} \cdot V^*_i}\right]^2}{V^*_m}\]

These are the same mixing rules used in all major commercial simulators.

Characteristic Volume (V*) Estimation

The characteristic volume $V^$ is the most important parameter for COSTALD accuracy. It is not the same as the critical volume $V_c$ — for polar compounds such as water or glycols, $V^$ can be 10–20% lower than $V_c$.

NeqSim determines $V^*$ using a three-tier priority:

Priority Source When Used
1 Explicit V* (setCostaldCharacteristicVolume) User has literature or fitted V* value
2 Back-calculated from normal liquid density at 60 °F (288.71 K) Component has normalLiquidDensity > 0 and $T_r < 0.9$ at 288.71 K
3 Critical volume $V_c$ Fallback when no density data is available

The back-calculation (priority 2) solves:

\[V^* = \frac{M / \rho_{std}}{V_R^{(0)}(T_{r,std}) \cdot \left[1 - \omega \cdot V_R^{(\delta)}(T_{r,std})\right]}\]

where $\rho_{std}$ is the normal liquid density at 288.71 K from the component database. This is the same approach used by UniSim/HYSYS and PRO/II. It automatically handles:

Supported Compound Types

Compound Type V* Method Example Components Expected Accuracy
Light hydrocarbons ($T_c < 320$ K) Critical volume $V_c$ methane, ethane, propane 1–3%
Heavier hydrocarbons V* from density n-hexane, n-decane, nC16, nC20 1–2%
Polar / associating V* from density water, methanol, ethanol 3–5%
Glycols V* from density MEG, TEG, DEG 3–5%
TBP fractions V* from density C7, C10, C20 pseudo-components 2–5%
Plus fractions (C20+) V* from density Characterized heavy end 3–8%

Aqueous and Multi-Phase Systems

setLiquidDensityModel("COSTALD") applies to all liquid-type phases simultaneously: LIQUID, OIL, and AQUEOUS. This means a single call enables COSTALD for three-phase gas-oil-water systems.

// Gas-oil-water with MEG inhibitor
SystemInterface fluid = new SystemSrkEos(303.15, 80.0);
fluid.addComponent("methane", 70.0);
fluid.addComponent("n-hexane", 5.0);
fluid.addComponent("water", 10.0);
fluid.addComponent("MEG", 7.0);
fluid.setMixingRule("classic");
fluid.setMultiPhaseCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.TPflash();
fluid.initPhysicalProperties();

// Enable COSTALD for both oil and aqueous phases
fluid.setLiquidDensityModel("COSTALD");

if (fluid.hasPhaseType("oil")) {
    double oilRho = fluid.getPhase("oil").getPhysicalProperties().getDensity();
}
if (fluid.hasPhaseType("aqueous")) {
    double aqRho = fluid.getPhase("aqueous").getPhysicalProperties().getDensity();
}

Oil Characterization (TBP / Plus Fractions)

COSTALD works directly with the NeqSim TBP characterization. Each pseudo-component’s $V^*$ is automatically estimated from its specified density:

SystemInterface fluid = new SystemSrkEos(313.15, 100.0);
fluid.addComponent("methane", 50.0);
fluid.addComponent("ethane", 10.0);
fluid.addTBPfraction("C7", 5.0, 95.0 / 1000.0, 0.738);
fluid.addTBPfraction("C10", 4.0, 134.0 / 1000.0, 0.792);
fluid.addPlusFraction("C20", 3.0, 350.0 / 1000.0, 0.895);
fluid.getCharacterization().setTBPModel("PedersenSRK");
fluid.getCharacterization().setLumpingModel("PVTlumpingModel");
fluid.getCharacterization().characterisePlusFraction();
fluid.setMixingRule("classic");

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.TPflash();
fluid.initPhysicalProperties();

// COSTALD will use the density-based V* for each TBP fraction
fluid.setLiquidDensityModel("COSTALD");
double density = fluid.getPhase("oil").getPhysicalProperties().getDensity();

Tuning with Custom V*

If you have published V* values (from Hankinson-Thomson 1979 Table I, DIPPR, or the API Technical Data Book), you can set them explicitly for better accuracy:

// Set published V* for n-hexane (371.0 cm3/mol from H-T 1979)
fluid.getPhase("oil").getComponent("n-hexane").setCostaldCharacteristicVolume(371.0);

// Set V* for water (46.4 cm3/mol, significantly less than Vc = 56 cm3/mol)
fluid.getPhase("aqueous").getComponent("water").setCostaldCharacteristicVolume(46.4);

Explicit V* values override the density-based estimation. To revert to automatic estimation, set it back to 0:

fluid.getPhase("oil").getComponent("n-hexane").setCostaldCharacteristicVolume(0.0);

Switching Back to Default

// Switch back to Peneloux volume shift (default)
fluid.setLiquidDensityModel("Peneloux");

// Or set per-phase
fluid.getPhase("oil").getPhysicalProperties().setDensityModel("Peneloux volume shift");

Valid Range and Limitations

Parameter Valid Range Notes
Reduced temperature $T_r$ 0.25 – 0.95 Warns below 0.25; falls back to EOS above 1.0
Pressure Up to ~700 bar Compressed liquid correction from Aalto (1996)
Component types Hydrocarbons, polar, associating, pseudo-components Best for hydrocarbons; 3–5% for polar
Phases Liquid, oil, aqueous Does not apply to gas or solid phases

Limitations:

Comparison with Other Commercial Simulators

Feature NeqSim COSTALD NeqSim NASTALD UniSim/HYSYS Aspen Plus PRO/II
Saturated volume (H-T 1979) Yes Yes Yes Yes Yes
NBS polar correction (1982) No Yes Yes Yes Yes
Compressed liquid Aalto (1996) Aalto (1996) Thomson (1982) Thomson (1982) Thomson (1982)
V* from density back-calc Yes Yes Yes DIPPR database Yes
V* explicit override Yes Yes Yes Yes Yes
Polar compound support Via V* from density Via V* + polar term Via V* + polar Via DIPPR V* Via V* + polar
Aqueous phase support Yes Yes Yes Yes Yes
TBP fraction support Yes Yes Yes Yes Yes
Mixing rules Standard H-T Standard H-T Standard H-T Standard H-T Standard H-T

NeqSim uses the Aalto et al. (1996) modified Tait equation for compressed liquid correction rather than the Thomson et al. (1982) original Tait form used by most commercial simulators. The Aalto form works in reduced pressure units for better numerical behavior at high pressures.

The NASTALD variant adds the Thomson-Brobst-Hankinson (1982) polar correction term, which is the same approach used by UniSim/HYSYS, Aspen Plus, and PRO/II for polar compounds.

References

  1. Hankinson, R.W. and Thomson, G.H. (1979). “A New Correlation for Saturated Densities of Liquids and Their Mixtures.” AIChE J. 25, 653–663.
  2. Thomson, G.H., Brobst, K.R. and Hankinson, R.W. (1982). “An Improved Correlation for Densities of Compressed Liquids and Liquid Mixtures.” AIChE J. 28, 671–676.
  3. Aalto, M., Keskinen, K.I., Aittamaa, J. and Liukkonen, S. (1996). “An Improved Correlation for Compressed Liquid Densities of Hydrocarbons. Part 2. Mixtures.” Fluid Phase Equil. 114, 1–19.
  4. Poling, B.E., Prausnitz, J.M. and O’Connell, J.P. (2001). The Properties of Gases and Liquids, 5th ed. McGraw-Hill, Chapters 4.
  5. API Technical Data Book — Petroleum Refining, Chapter 6.

NASTALD (COSTALD with Polar Correction)

NASTALD extends the standard COSTALD method by adding the Thomson-Brobst-Hankinson (1982) polar correction term $V_R^{(p)}$. This improves accuracy for strongly polar and hydrogen-bonding compounds such as water, alcohols, and glycols.

Equation

The saturated volume for a polar substance becomes:

\[V_s = V^* \left[ V_R^{(0)} \left(1 - \omega_{SRK} \, V_R^{(\delta)}\right) + \omega_p \, V_R^{(p)} \right]\]

where $\omega_p$ is the substance-specific polar parameter and $V_R^{(p)}$ is the polar correction function:

\[V_R^{(p)} = \frac{e + f \, T_r + g \, T_r^2 + h \, T_r^3}{T_r - 1.00001}\]

with the universal constants:

Constant Value
$e$ -1.52816
$f$ 1.43907
$g$ -0.81446
$h$ 0.190454

For mixtures, the polar parameter is mixed via:

\[\omega_{p,mix} = \sum_i x_i \, \omega_{p,i}\]

Polar Parameters ($\omega_p$)

The following polar parameters are built into NeqSim:

Compound $\omega_p$ Source
Water 0.3478 Thomson et al. (1982)
Ammonia 0.2872 Thomson et al. (1982)
HF 0.3750 Thomson et al. (1982)
Methanol 0.1907 Thomson et al. (1982)
Ethanol 0.1471 Thomson et al. (1982)
1-propanol 0.1100 Thomson et al. (1982)
2-propanol 0.1050 Thomson et al. (1982)
1-butanol 0.0922 Thomson et al. (1982)
Acetic acid 0.1530 Thomson et al. (1982)
Acetone 0.0547 Thomson et al. (1982)
MEG 0.2213 Estimated
DEG 0.2000 Estimated
TEG 0.1800 Estimated

For compounds not in this table, $\omega_p = 0$ and NASTALD reduces to standard COSTALD.

Important Caveat: V* Back-Calculation

When V* is back-calculated from the component’s normalLiquidDensity (which is the default for most components), the polarity is already partially captured in the V* value. Adding the polar correction on top can lead to over-correction (typically 5-15% too high for alcohols and glycols).

The NASTALD polar correction is most beneficial when:

For most practical purposes, standard COSTALD with V* from density gives excellent results for polar compounds without needing the explicit polar term.

Usage

// Enable NASTALD (COSTALD with polar correction) for all liquid phases
fluid.setLiquidDensityModel("NASTALD");

// Run flash and get density
ops.TPflash();
fluid.initPhysicalProperties();
double rho = fluid.getPhase("aqueous").getDensity("kg/m3");

References

  1. Thomson, G.H., Brobst, K.R. and Hankinson, R.W. (1982). “An Improved Correlation for Densities of Compressed Liquids and Liquid Mixtures.” AIChE J. 28, 671–676.

Rackett Equation

The Spencer-Danner (1972) modified Rackett equation provides a simple corresponding-states method for saturated liquid density. It requires only critical properties and the Rackett compressibility factor $Z_{RA}$.

Pure Component Equation

\[V_s = \frac{R T_c}{P_c} \; Z_{RA}^{\left[1 + (1 - T_r)^{2/7}\right]}\]

where $T_r = T / T_c$ is the reduced temperature and $Z_{RA}$ is the Rackett compressibility factor (an empirically fitted parameter, not the true critical compressibility $Z_c$).

$Z_{RA}$ Selection Priority

  1. Database value — component.getRacketZ() from the NeqSim component database (fitted to experimental data)
  2. Yamada-Gunn estimate (1973) — If no database value is available:
\[Z_{RA} = 0.29056 - 0.08775 \, \omega\]

where $\omega$ is the acentric factor.

Mixture Rules (Li, 1971)

For mixtures, pseudo-critical properties are calculated via mole-fraction weighted mixing:

\[V_{c,mix} = \sum_i x_i V_{c,i}, \quad T_{c,mix} = \frac{\sum_i x_i V_{c,i} T_{c,i}}{V_{c,mix}}, \quad Z_{RA,mix} = \sum_i x_i Z_{RA,i}\] \[P_{c,mix} = \frac{Z_{RA,mix} \, R \, T_{c,mix}}{V_{c,mix}}\]

The mixture saturated volume is then:

\[V_{s,mix} = \frac{R \, T_{c,mix}}{P_{c,mix}} \; Z_{RA,mix}^{\left[1 + (1 - T_{r,mix})^{2/7}\right]}\]

Usage

// Enable Rackett density model for all liquid phases
fluid.setLiquidDensityModel("Rackett");

// Run flash and get density
ops.TPflash();
fluid.initPhysicalProperties();
double rho = fluid.getPhase("oil").getDensity("kg/m3");

// Access Z_RA for a component
double Zra = fluid.getPhase(1).getComponent("n-pentane").getRacketZ();

Strengths and Limitations

References

  1. Rackett, H.G. (1970). “Equation of State for Saturated Liquids.” J. Chem. Eng. Data 15, 514–517.
  2. Spencer, C.F. and Danner, R.P. (1972). “Improved Equation for Prediction of Saturated Liquid Density.” J. Chem. Eng. Data 17, 236–241.
  3. Li, C.C. (1971). “Critical Temperature Estimation for Simple Mixtures.” Can. J. Chem. Eng. 49, 709–710.
  4. Yamada, T. and Gunn, R.D. (1973). “Saturated Liquid Molar Volumes. The Rackett Equation.” J. Chem. Eng. Data 18, 234–236.

Binary Electrolyte Volumetric Pitzer Kernel

PitzerBinaryVolumetricModel evaluates the standard pressure derivative of the binary-electrolyte Pitzer excess Gibbs energy. It is a parameter-neutral thermodynamic kernel: it does not contain a built-in coefficient table and it is not selected by setLiquidDensityModel(...).

For a binary salt $M_{\nu_M}X_{\nu_X}$ with formula-unit molality $m$,

\[\begin{aligned} \phi_V={}&V^\circ+\nu|z_Mz_X|\frac{A_V}{2b}\ln(1+b\sqrt{I})\\ &+\nu_M\nu_XRT\left[2mB^V_{MX}+m^2\sqrt{\nu_M\nu_X}C^{\phi V}_{MX}\right],\\ B^V_{MX}={}&\beta^{(0)V}_{MX}+\beta^{(1)V}_{MX}g(\alpha\sqrt{I}),\\ g(x)={}&\frac{2[1-(1+x)e^{-x}]}{x^2},\\ I={}&\frac{1}{2}\nu|z_Mz_X|m. \end{aligned}\]

The implementation uses the standard $b=1.2$ and $\alpha=2.0$ values and an analytical small-$x$ expansion for $g(x)$ to prevent dilute-limit cancellation. All inputs use SI units. The pressure derivatives $\beta^{(0)V}$, $\beta^{(1)V}$, and $C^{\phi V}$ must already be evaluated at the temperature and pressure stored in StateParameters.

The non-zero reference form is also available:

\[\phi_V(m)=\phi_V(m_r)+F(m)-F(m_r),\]

where $F$ is the Debye–Hückel plus binary-interaction contribution in the equation above. This form avoids deriving $V^\circ$ from small differences of dilute-solution densities.

Density conversion uses an exact one-kilogram-solvent balance:

\[\rho=\frac{1+mM}{1/\rho_w+m\phi_V}.\]

The inverse conversion is provided for auditable data preparation.

PitzerBinaryVolumetricModel calciumChloride =
    new PitzerBinaryVolumetricModel(1, 2, 2, -1);

PitzerBinaryVolumetricModel.StateParameters state =
    new PitzerBinaryVolumetricModel.StateParameters(
        temperatureK,
        pressurePa,
        debyeHuckelVolumeSlope,
        beta0PressureDerivative,
        beta1PressureDerivative,
        cphiPressureDerivative);

double apparentMolarVolume = calciumChloride.calculateApparentMolarVolumeFromReference(
    molality,
    referenceMolality,
    referenceApparentMolarVolume,
    state);

double density = PitzerBinaryVolumetricModel.calculateDensity(
    molality,
    calciumChlorideMolarMass,
    pureWaterDensity,
    apparentMolarVolume);

Parameter regression and identifiability

PitzerBinaryVolumetricRegression performs weighted linear regression for one common temperature-pressure state. The caller supplies every apparent-molar- volume observation, its absolute one-sigma uncertainty, a laboratory or source-lineage identifier, and the Debye-Hückel volume slope. The regression fits V°, β⁽⁰⁾_V, β⁽¹⁾_V, and Cφ_V without installing the result in a phase or changing a default model.

The weighted design matrix is column-scaled and solved by singular-value decomposition. A fit fails closed when the concentration design is rank deficient or excessively ill-conditioned. FitResult reports coefficient covariance and standard uncertainties, chi-square, reduced chi-square, weighted RMS residual, maximum standardized residual, and deterministic per-source-group residual diagnostics. Covariance assumes that input uncertainties are absolute one-sigma values.

Regression capability does not make an input dataset acceptable. Callers must separately establish row provenance, redistribution rights, uncertainty meaning, species and composition basis, validity range, and an independent laboratory-grouped validation split. The reserved Al Ghafri/NIST ThermoML pressure series must not be supplied as calibration input when it is retained as the campaign hold-out.

Independent source-group holdout validation

PitzerBinaryVolumetricGroupedValidation fits only an explicit calibration list and evaluates a separate untouched holdout list. It rejects any source-lineage identifier present in both lists, reports holdout chi-square, weighted RMS and maximum standardized residual, and returns deterministic per-lineage diagnostics. The class deliberately defines no universal pass threshold.

The string identifiers are an auditable software boundary, not proof of scientific independence. The caller must map each row to its real laboratory, apparatus, publication lineage and reuse history before assigning groups. A validation result is admissible only when those lineages are genuinely independent, all rows and uncertainties have qualified provenance, and the acceptance limits were fixed before the holdout was evaluated.

Repository dataset provenance

PitzerBinaryVolumetricDatasetProvenance records a machine-auditable manifest for observations distributed with NeqSim. Each source-lineage record fixes its calibration or validation role, full citation and stable URL, license and redistribution decision, SHA-256 source-file checksum, uncertainty basis and explicit qualification, exact row count, and molality, temperature and pressure envelope.

Rows are repository-admissible only when the manifest explicitly marks their uncertainty mapping as QUALIFIED_ABSOLUTE_ONE_SIGMA. Instrument accuracy, precision, or repeatability specifications and unresolved empirical error envelopes are non-qualified states and fail before fitting or holdout evaluation.

validateRepositoryDatasets(...) checks both observation lists against that manifest before fitting. It fails closed for an undeclared lineage, role or row count mismatch, a state outside the declared envelope, or any source whose redistribution status is restricted or unknown. The original validate(...) method remains available for caller-owned private or in-memory data that are not bundled with NeqSim.

The manifest verifies declared metadata against the supplied observations. It does not hash source bytes, prove that a license interpretation is correct, or prove experimental independence. Repository maintainers must separately compare the recorded checksum with the distributed file, audit the permission decision, map source groups to real laboratories and apparatus, and pre-register acceptance limits before evaluating a holdout.

Acoustic compressibility conversion

PitzerBinaryVolumetricAcousticConversion provides the missing thermodynamic bridge for assessing sound-speed evidence without equating isentropic and isothermal response:

\[\kappa_S=\frac{1}{\rho c^2},\qquad \kappa_T=\kappa_S+\frac{T\alpha^2}{\rho c_p}.\]

Inputs use SI units: temperature in K, density in kg/m3, speed of sound in m/s, volumetric thermal expansivity in 1/K, and mass-specific isobaric heat capacity in J/(kg K). Compressibilities are returned in 1/Pa. The uncertainty method accepts a full 5-by-5 covariance matrix in that input order, checks that it is finite, symmetric, and positive semidefinite, and evaluates $u^2(\kappa_T)=J\Sigma J^T$ with an analytic sensitivity Jacobian. This retains correlations among measurements instead of silently assuming independence.

The conversion does not fill missing density, expansivity, or heat-capacity inputs, infer their covariance, fit any coefficient, or activate a density model. Acoustic rows therefore remain inadmissible for volumetric-Pitzer calibration until every matched property, uncertainty, source lineage, and redistribution right passes the campaign provenance gates.

The SI conversion is checked against the CC BY 4.0 compressed-water results of El Hawary and Meier (2023) at 303.15 K and 50 MPa. Their published sound-speed correlation and Table 5 density and heat-capacity values give $\kappa_T=3.96032\times10^{-10}$ 1/Pa. This agrees within 0.018% with the alternative $\kappa_Sc_p/c_v$ identity and within 0.030% with a symmetric 45/55 MPa density derivative. The benchmark validates units and thermodynamic consistency of the conversion. Because the source derived density and heat capacity through thermodynamic integration of its acoustic measurements, this is not an independent experimental validation. It is pure-water evidence and does not qualify CaCl2 acoustic data or any volumetric-Pitzer coefficient.

Reference: El Hawary and Meier (2023), doi:10.1007/s10765-023-03276-1.

Qualification boundary

Reference: Rowland and May (2013), doi:10.1016/j.fluid.2012.10.021.

Usage Example

The following complete Java 8 program compares the maintained liquid-density models without changing the equation of state. Constructor inputs are kelvin and bara. Assertions are deliberately broad physical checks: model selection still requires independent data representative of the fluid and operating envelope.

import neqsim.thermo.system.SystemInterface;
import neqsim.thermo.system.SystemSrkEos;
import neqsim.thermodynamicoperations.ThermodynamicOperations;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;

public final class DensityModelComparison {
  private static final Logger logger = LogManager.getLogger(DensityModelComparison.class);

  private DensityModelComparison() {}

  public static void main(String[] args) {
    SystemInterface hydrocarbon = createLiquid("n-hexane", 298.15, 10.0);
    String oilPhase = "oil";

    double penelouxDensity = density(hydrocarbon, oilPhase);
    hydrocarbon.setLiquidDensityModel("COSTALD");
    double costaldDensity = density(hydrocarbon, oilPhase);
    hydrocarbon.setLiquidDensityModel("Rackett");
    double rackettDensity = density(hydrocarbon, oilPhase);
    hydrocarbon.setLiquidDensityModel("Peneloux");
    double restoredDensity = density(hydrocarbon, oilPhase);

    assertPhysicalDensity("Peneloux", penelouxDensity, 500.0, 800.0);
    assertPhysicalDensity("COSTALD", costaldDensity, 500.0, 800.0);
    assertPhysicalDensity("Rackett", rackettDensity, 500.0, 800.0);
    assert Math.abs(restoredDensity - penelouxDensity) < 0.01
        : "Restoring Peneloux must restore the original density";

    double originalTemperatureK = hydrocarbon.getTemperature("K");
    double originalPressureBara = hydrocarbon.getPressure("bara");
    double referenceDensity =
        hydrocarbon.getDensityAtReferenceConditions(15.0, "C", 1.01325, "bara");
    assertPhysicalDensity("15 C reference", referenceDensity, 600.0, 720.0);
    assert Math.abs(hydrocarbon.getTemperature("K") - originalTemperatureK) < 1.0e-10
        : "Reference-condition calculation changed the live temperature";
    assert Math.abs(hydrocarbon.getPressure("bara") - originalPressureBara) < 1.0e-10
        : "Reference-condition calculation changed the live pressure";

    SystemInterface water = createLiquid("water", 293.15, 10.0);
    String waterPhase = water.hasPhaseType("aqueous") ? "aqueous" : "oil";
    water.setLiquidDensityModel("COSTALD");
    double waterCostaldDensity = density(water, waterPhase);
    water.setLiquidDensityModel("NASTALD");
    double waterNastaldDensity = density(water, waterPhase);

    assertPhysicalDensity("water COSTALD", waterCostaldDensity, 900.0, 1100.0);
    assertPhysicalDensity("water NASTALD", waterNastaldDensity, 900.0, 1100.0);
    assert Math.abs(waterCostaldDensity - waterNastaldDensity) > 0.1
        : "The polar correction should change the water result";

    logger.info(
        "n-hexane density kg/m3: Peneloux={}, COSTALD={}, Rackett={}, 15 C reference={}",
        penelouxDensity,
        costaldDensity,
        rackettDensity,
        referenceDensity);
    logger.info(
        "water density kg/m3: COSTALD={}, NASTALD={}",
        waterCostaldDensity,
        waterNastaldDensity);
  }

  private static SystemInterface createLiquid(
      String component, double temperatureK, double pressureBara) {
    SystemInterface fluid = new SystemSrkEos(temperatureK, pressureBara);
    fluid.addComponent(component, 1.0);
    fluid.setMixingRule("classic");
    new ThermodynamicOperations(fluid).TPflash();
    fluid.initPhysicalProperties();
    return fluid;
  }

  private static double density(SystemInterface fluid, String phaseName) {
    double value = fluid.getPhase(phaseName).getPhysicalProperties().calcDensity();
    assert Double.isFinite(value) : phaseName + " density is not finite";
    return value;
  }

  private static void assertPhysicalDensity(
      String label, double value, double lowerBound, double upperBound) {
    assert value > lowerBound && value < upperBound
        : label + " density outside the demonstration bounds: " + value;
  }
}

Run with assertions enabled so the physical and state-preservation checks are active. The exact program is extracted, compiled with Java 8 source/target settings, and executed by the documentation contract.

Interpretation and calibration boundary

Model Selection Guide

Situation Recommended Model Notes
General hydrocarbons Peneloux Default, good accuracy
Near saturation COSTALD Better for saturated liquids
Polar compounds (water, glycols) COSTALD V* from density handles polarity
Polar with database V* NASTALD Adds explicit polar correction term
Binary electrolyte with qualified volumetric parameters PitzerBinaryVolumetricModel Low-level opt-in kernel; no bundled coefficient set
Multi-phase (gas-oil-water) COSTALD Single call applies to oil + aqueous
Different model per phase Per-phase API e.g. COSTALD on oil, Peneloux on aqueous
TBP/plus fractions COSTALD V* from density input, no tuning needed
High pressure (> 200 bar) COSTALD Aalto compressed liquid correction
Quick simple estimate Rackett No high-pressure correction
Saturated hydrocarbons only Rackett Fast, 2–5% accuracy
Critical region GERG-2008 If available
Quick estimate EoS only 5–15% error typical

Expected Accuracy

Method Liquid Density Error Vapor Density Error
SRK (no correction) 5–15% 1–3%
SRK + Peneloux 1–3% 1–3%
PR (no correction) 3–10% 1–3%
PR + Peneloux 1–3% 1–3%
COSTALD (hydrocarbons) 1–2% N/A
COSTALD (polar/aqueous) 3–5% N/A
COSTALD (TBP fractions) 2–5% N/A
NASTALD (polar with database V*) 1–3% N/A
Rackett (hydrocarbons) 2–5% N/A
GERG-2008 0.1–0.5% 0.1–0.5%
Binary volumetric Pitzer kernel Dataset-dependent N/A

API Reference

Setting Density Model

// Set COSTALD for all liquid phases (oil, liquid, aqueous)
fluid.setLiquidDensityModel("COSTALD");

// Set NASTALD (polar-corrected COSTALD) for all liquid phases
fluid.setLiquidDensityModel("NASTALD");

// Set Rackett equation for all liquid phases
fluid.setLiquidDensityModel("Rackett");

// Switch back to default Peneloux
fluid.setLiquidDensityModel("Peneloux");

// Set model for a specific phase type only
fluid.setLiquidDensityModel("COSTALD", "oil");      // Oil phase uses COSTALD
fluid.setLiquidDensityModel("Peneloux", "aqueous");  // Aqueous keeps Peneloux

// Valid phase type names: "oil", "aqueous", "liquid"

// Low-level: set per-phase directly
fluid.getPhase("oil").getPhysicalProperties().setDensityModel("Costald");
fluid.getPhase("aqueous").getPhysicalProperties().setDensityModel("Rackett");

Accepted Model Strings

setLiquidDensityModel(...) setDensityModel(...) Description
"COSTALD" "Costald" Standard COSTALD (Hankinson-Thomson)
"NASTALD" or "COSTALD-polar" "Costald polar" COSTALD with polar correction
"Rackett" "Rackett" Spencer-Danner modified Rackett
"Peneloux" "Peneloux volume shift" Default Peneloux volume translation

Accessing Density

// Mass density
double rhoMass = phase.getDensity("kg/m3");
double rhoMass2 = phase.getDensity("lb/ft3");

// Molar density
double rhoMolar = phase.getDensity("mol/m3");

// Molar volume
double Vm = phase.getMolarVolume();  // m³/mol

Density at Reference Conditions

// Calculate density at standard conditions without modifying the fluid
double rhoStd = fluid.getDensityAtReferenceConditions(15.0, "C", 1.01325, "bara");
// Returns density in kg/m3 at 15°C / 1.01325 bara

// API standard temperature (60°F)
double rhoAPI = fluid.getDensityAtReferenceConditions(15.56, "C", 1.01325, "bara");

// The original fluid state is NOT modified by this call

Volume Correction Parameters

// Get/set volume correction constant
double c = component.getVolumeCorrectionConst();
component.setVolumeCorrectionConst(newValue);

// Get Rackett parameter
double Zra = component.getRacketZ();

// Get/set COSTALD characteristic volume (V*)
double Vstar = component.getCostaldCharacteristicVolume();
component.setCostaldCharacteristicVolume(newValue);  // cm³/mol

Binary Volumetric Pitzer API

Method Purpose
calculateApparentMolarVolume(...) Evaluate the standard infinite-dilution form
calculateApparentMolarVolumeFromReference(...) Evaluate the numerically stable non-zero reference form
calculateDensity(...) Convert apparent molar volume to density on a 1 kg solvent basis
calculateApparentMolarVolumeFromDensity(...) Invert a density observation to apparent molar volume
calculateIonicStrength(...) Return the binary salt’s stoichiometric ionic strength
PitzerBinaryVolumetricRegression.fit(...) Fit caller-supplied one-state observations with SVD rank checks
FitResult.getGroupStatistics() Inspect residuals by laboratory or source lineage
PitzerBinaryVolumetricGroupedValidation.validate(...) Fit one lineage set and evaluate a disjoint untouched holdout
PitzerBinaryVolumetricDatasetProvenance.validateRepositoryObservations(...) Check role, license, checksum manifest, row count and state envelope
PitzerBinaryVolumetricGroupedValidation.validateRepositoryDatasets(...) Enforce repository provenance before grouped fitting and holdout evaluation
PitzerBinaryVolumetricAcousticConversion.calculateIsentropicCompressibility(...) Convert density and sound speed to isentropic compressibility
PitzerBinaryVolumetricAcousticConversion.calculateIsothermalCompressibility(...) Apply the expansivity and heat-capacity correction in SI units
PitzerBinaryVolumetricAcousticConversion.convertWithUncertainty(...) Propagate a full input covariance through the acoustic conversion

These methods do not install a parameter dataset or alter a phase. A caller must explicitly supply a provenance-qualified StateParameters instance.

References

  1. Peneloux, A., Rauzy, E., Freze, R. (1982). A Consistent Correction for Redlich-Kwong-Soave Volumes. Fluid Phase Equilib. 8, 7–23.
  2. Hankinson, R.W. and Thomson, G.H. (1979). A New Correlation for Saturated Densities of Liquids and Their Mixtures. AIChE J. 25, 653–663.
  3. Thomson, G.H., Brobst, K.R. and Hankinson, R.W. (1982). An Improved Correlation for Densities of Compressed Liquids and Liquid Mixtures. AIChE J. 28, 671–676.
  4. Aalto, M., Keskinen, K.I., Aittamaa, J. and Liukkonen, S. (1996). An Improved Correlation for Compressed Liquid Densities of Hydrocarbons. Part 2. Mixtures. Fluid Phase Equil. 114, 1–19.
  5. Rackett, H.G. (1970). Equation of State for Saturated Liquids. J. Chem. Eng. Data 15, 514–517.
  6. Spencer, C.F. and Danner, R.P. (1972). Improved Equation for Prediction of Saturated Liquid Density. J. Chem. Eng. Data 17, 236–241.
  7. Yamada, T. and Gunn, R.D. (1973). Saturated Liquid Molar Volumes. The Rackett Equation. J. Chem. Eng. Data 18, 234–236.
  8. Li, C.C. (1971). Critical Temperature Estimation for Simple Mixtures. Can. J. Chem. Eng. 49, 709–710.
  9. Jhaveri, B.S. and Youngren, G.K. (1988). Three-Parameter Modification of the Peng-Robinson Equation of State. SPE Reservoir Eng. 3, 1033–1040.
  10. Poling, B.E., Prausnitz, J.M. and O’Connell, J.P. (2001). The Properties of Gases and Liquids, 5th ed. McGraw-Hill.
  11. Rowland, D. and May, P.M. (2013). A Pitzer-Based Characterization of Aqueous Magnesium Chloride, Calcium Chloride and Potassium Iodide Solution Densities to High Temperature and Pressure. Fluid Phase Equilib. 338, 54–62. doi:10.1016/j.fluid.2012.10.021.