Fluid Characterization in NeqSim: Mathematical Foundations

This document provides detailed mathematical documentation of the fluid characterization methods implemented in NeqSim, with emphasis on plus fraction (C7+) modeling, TBP fraction property correlations, and pseudo-component generation.

Table of Contents

  1. Overview
  2. Plus Fraction Distribution Models
  3. Critical Property Correlations
  4. Density Correlations
  5. Boiling Point Correlations
  6. Lumping Methods
  7. Usage Examples
  8. Common Fluid Characterization

Overview

Petroleum fluids contain thousands of hydrocarbon compounds. For practical thermodynamic calculations, heavy fractions (typically C7+) must be characterized using continuous distribution functions or discrete pseudo-components. NeqSim implements industry-standard characterization methods from Whitson (1983), Pedersen et al. (1984), and others.

The characterization workflow consists of three main steps:

  1. Plus Fraction Splitting: Distribute the C7+ fraction into discrete Single Carbon Number (SCN) groups
  2. Property Estimation: Calculate critical properties (Tc, Pc, ω) for each pseudo-component
  3. Lumping: Combine pseudo-components into computationally efficient groups

Plus Fraction Distribution Models

Whitson Gamma Distribution Model

The Whitson gamma distribution (SPE 12233, 1983) is the most widely used continuous distribution for petroleum C7+ fractions.

Probability Density Function

The molar distribution of the plus fraction is described by a three-parameter gamma distribution:

\[p(M) = \frac{(M - \eta)^{\alpha - 1}}{\beta^\alpha \cdot \Gamma(\alpha)} \exp\left(-\frac{M - \eta}{\beta}\right)\]

where:

Gamma Function Approximation

NeqSim uses a polynomial approximation for the gamma function:

\[\Gamma(\alpha) \approx \alpha \cdot \left(1 + \sum_{i=1}^{8} b_i \cdot (\alpha - 1)^i\right)\]

with coefficients: | $i$ | $b_i$ | |—–|——-| | 1 | -0.577191652 | | 2 | 0.988205891 | | 3 | -0.897056937 | | 4 | 0.918206857 | | 5 | -0.756704078 | | 6 | 0.482199394 | | 7 | -0.193527818 | | 8 | 0.035868343 |

Parameter Relationships

The mean of the gamma distribution equals the plus fraction average molecular weight:

\[\bar{M}_{C7+} = \eta + \alpha \cdot \beta\]

Therefore, the scale parameter is calculated as:

\[\beta = \frac{\bar{M}_{C7+} - \eta}{\alpha}\]

Cumulative Distribution Functions

For splitting the plus fraction into SCN groups, NeqSim computes:

\[P_0(M_b) = \int_\eta^{M_b} p(M) \, dM = Q(M_b) \cdot S(M_b)\]

where:

\[Q = \frac{\exp(-Y) \cdot Y^\alpha}{\Gamma(\alpha)}, \quad Y = \frac{M_b - \eta}{\beta}\] \[S = \sum_{j=0}^{\infty} \frac{Y^j}{\prod_{k=0}^{j}(\alpha + k)} = \frac{1}{\alpha} + \frac{Y}{\alpha(\alpha+1)} + \frac{Y^2}{\alpha(\alpha+1)(\alpha+2)} + \cdots\]

The first moment $P_1$ is used for average molecular weight calculation:

\[P_1(M_b) = Q(M_b) \cdot \left(S(M_b) - \frac{1}{\alpha}\right)\]

Mole Fraction of SCN Group

For an SCN group bounded by molecular weights $M_L$ (lower) and $M_U$ (upper):

\[z_i = z_{C7+} \cdot \left[P_0(M_U) - P_0(M_L)\right]\]

Average Molecular Weight of SCN Group

\[\bar{M}_i = \eta + \alpha \cdot \beta \cdot \frac{P_1(M_U) - P_1(M_L)}{P_0(M_U) - P_0(M_L)}\]

Shape Parameter Guidelines

The shape parameter $\alpha$ characterizes the fluid type:

Fluid Type Typical $\alpha$ Range Watson $K_w$
Gas condensates 0.5 - 1.0 > 12.5
Black oils 1.0 - 2.0 11.5 - 12.5
Heavy/naphthenic oils 2.0 - 4.0 < 11.5

Automatic Alpha Estimation

NeqSim can automatically estimate $\alpha$ using the Watson characterization factor:

\[K_w = 4.5579 \cdot M_{C7+}^{0.15178} \cdot \rho_{C7+}^{-1.18241}\]

where $\rho$ is specific gravity (g/cm³). The estimated alpha is:

\[\alpha = \begin{cases} 0.5 + 0.1(K_w - 12.5) & K_w \geq 12.5 \text{ (paraffinic)} \\ 1.0 + 0.5(K_w - 11.5) & 11.5 \leq K_w < 12.5 \text{ (mixed)} \\ 1.5 + 0.5(K_w - 10.5) & 10.5 \leq K_w < 11.5 \text{ (naphthenic)} \\ 2.0 + 0.5(10.5 - K_w) & K_w < 10.5 \text{ (aromatic)} \end{cases}\]

Pedersen Exponential Model

The Pedersen model (Pedersen et al., 1984) uses an exponential distribution:

\[\ln(z_n) = A + B \cdot n\]

where:

The coefficients are determined by matching:

  1. The plus fraction mole fraction $z_{C7+}$
  2. The plus fraction molecular weight $M_{C7+}$
  3. The plus fraction density $\rho_{C7+}$

Density Correlation (Pedersen)

\[\rho_n = C_1 + C_2 \cdot \ln(n)\]

where $C_1$ and $C_2$ are fitted to match the measured plus fraction density.

Critical Property Correlations

Pedersen Correlations

For SRK equation of state, critical properties are calculated using Pedersen et al. correlations:

Critical Temperature (K)

\[T_c = a_0 \cdot \rho + a_1 \cdot \ln(M) + a_2 \cdot M + \frac{a_3}{M}\]
Coefficient Light Oil (M < 1120 g/mol) Heavy Oil
$a_0$ 163.12 830.63
$a_1$ 86.052 17.5228
$a_2$ 0.43475 0.0455911
$a_3$ -1877.4 -11348.4

Critical Pressure (bar)

\[P_c = \exp\left(0.01325 + b_0 + b_1 \cdot \rho^{b_4} + \frac{b_2}{M} + \frac{b_3}{M^2}\right)\]
Coefficient Light Oil Heavy Oil
$b_0$ -0.13408 0.802988
$b_1$ 2.5019 1.78396
$b_2$ 208.46 156.740
$b_3$ -3987.2 -6965.59
$b_4$ 1.0 0.25

Acentric Factor (SRK m-parameter)

\[m = c_0 + c_1 \cdot M + c_2 \cdot \rho + c_3 \cdot M^2\]

For Peng-Robinson, different coefficients are used.

Riazi-Daubert Correlations

Alternative correlations based on molecular weight and specific gravity:

Critical Temperature (K)

\[T_c = \frac{5}{9} \cdot 554.4 \cdot \exp(-1.3478 \times 10^{-4} M - 0.61641 \cdot \rho) \cdot M^{0.2998} \cdot \rho^{1.0555}\]

Critical Pressure (bar)

\[P_c = 0.068947 \cdot 4.5203 \times 10^4 \cdot \exp(-1.8078 \times 10^{-3} M - 0.3084 \cdot \rho) \cdot M^{-0.8063} \cdot \rho^{1.6015}\]

Kesler-Lee Correlations

Acentric Factor

For $T_{br} = T_b/T_c < 0.8$:

\[\omega = \frac{\ln(P_{br}) - 5.92714 + \frac{6.09649}{T_{br}} + 1.28862 \ln(T_{br}) - 0.169347 T_{br}^6}{15.2518 - \frac{15.6875}{T_{br}} - 13.4721 \ln(T_{br}) + 0.43577 T_{br}^6}\]

For $T_{br} \geq 0.8$:

\[\omega = -7.904 + 0.1352 K_w - 0.007465 K_w^2 + 8.359 T_{br} + \frac{1.408 - 0.01063 K_w}{T_{br}}\]

where $P_{br} = 1.01325 / P_c$ and $K_w = T_b^{1/3} / \rho$ is the Watson K-factor.

Density Correlations

Watson K-Factor Method (UOP)

The Universal Oil Products (UOP) characterization assumes constant Watson K for all SCN groups:

\[K_w = 4.5579 \cdot M_{C7+}^{0.15178} \cdot \rho_{C7+}^{-1.18241}\]

Individual SCN densities:

\[\rho_i = 6.0108 \cdot M_i^{0.17947} \cdot K_w^{-1.18241}\]

Limitation: Less accurate for heavy fractions (C20+).

Søreide Correlation

Søreide (1989) developed a more accurate correlation for heavy fractions:

\[SG = 0.2855 + C_f \cdot (M - 66)^{0.13}\]

where $SG$ is specific gravity (same as $\rho$ in g/cm³) and $C_f$ is calculated from the plus fraction:

\[C_f = \frac{SG_{C7+} - 0.2855}{(M_{C7+} - 66)^{0.13}}\]

Individual SCN densities:

\[\rho_i = 0.2855 + C_f \cdot (M_i - 66)^{0.13}\]

Constrained to physical limits: $0.6 \leq \rho_i \leq 1.2$ g/cm³.

Boiling Point Correlations

Søreide Correlation

\[T_b = \frac{1}{1.8}\left(1928.3 - 1.695 \times 10^5 \cdot M^{-0.03522} \cdot \rho^{3.266} \cdot \exp\left(-4.922 \times 10^{-3} M - 4.7685 \rho + 3.462 \times 10^{-3} M \cdot \rho\right)\right)\]

Riazi-Daubert Correlation

\[T_b = \left(\frac{M}{5.805 \times 10^{-5} \cdot \rho^{0.9371}}\right)^{1/2.3776}\]

Polynomial Correlation (Light Components)

For $M < 540$ g/mol:

\[T_b = 2 \times 10^{-6} M^3 - 0.0035 M^2 + 2.4003 M + 171.74\]

Lumping Methods

Lumping reduces computational cost by grouping many SCN (Single Carbon Number) pseudo-components into fewer lumped components while preserving bulk properties.

Available Lumping Models in NeqSim

Model Name API Name Description
PVT Lumping Model "PVTlumpingModel" Default. Preserves TBP fractions (C6-C9) as separate pseudo-components, only lumps the plus fraction (C10+)
Standard Lumping Model "standard" Lumps all TBP fractions and plus fractions together starting from C6
No Lumping "no lumping" Keeps all individual SCN components (C6, C7, … C80)

NeqSim provides a fluent builder API for configuring lumping, which makes the intent clearer and avoids confusion between parameters:

// PVTlumpingModel: keep C6-C9 separate, lump C10+ into 5 groups
fluid.getCharacterization().configureLumping()
    .model("PVTlumpingModel")
    .plusFractionGroups(5)
    .build();

// Standard model: create exactly 6 total pseudo-components from C6+
fluid.getCharacterization().configureLumping()
    .model("standard")
    .totalPseudoComponents(6)
    .build();

// No lumping: keep all individual SCN components
fluid.getCharacterization().configureLumping()
    .noLumping()
    .build();

// Custom boundaries: match specific PVT lab report groupings
// Creates groups: C6, C7-C9, C10-C14, C15-C19, C20+
fluid.getCharacterization().configureLumping()
    .customBoundaries(6, 7, 10, 15, 20)
    .build();

Legacy Configuration Parameters

The lumping model has two key parameters:

Parameter Method What it Controls
numberOfPseudoComponents setNumberOfPseudoComponents(n) Total number of pseudo-components (TBP + lumped)
numberOfLumpedComponents setNumberOfLumpedComponents(n) Number of groups created from the plus fraction only

⚠️ Deprecation Notice: For PVTlumpingModel, the method setNumberOfPseudoComponents() is deprecated. Use setNumberOfLumpedComponents() or the fluent API plusFractionGroups() instead.

Model Fluent API Method Legacy Method
"standard" totalPseudoComponents(n) setNumberOfPseudoComponents(n)
"PVTlumpingModel" plusFractionGroups(n) setNumberOfLumpedComponents(n)
Custom grouping customBoundaries(...) setCustomBoundaries(int[])

Quick Reference

I want to… Model Fluent API
Get exactly N total pseudo-components (lumping from C6) "standard" .model("standard").totalPseudoComponents(N)
Keep C6-C9 separate, lump C10+ into N groups "PVTlumpingModel" .model("PVTlumpingModel").plusFractionGroups(N)
Match PVT lab report groupings (e.g., C6, C7-C9, C10+) Any .customBoundaries(6, 7, 10)
Keep all SCN components (C6-C80) "no lumping" .noLumping()

Custom Carbon Number Boundaries

When matching specific PVT lab report groupings, use custom boundaries to specify exactly which carbon numbers start each group:

// Match a PVT report with groups: C6, C7-C9, C10-C14, C15-C19, C20+
fluid.getCharacterization().configureLumping()
    .customBoundaries(6, 7, 10, 15, 20)
    .build();

Each value represents the starting carbon number for a group. The final group extends to the heaviest component (typically C80).

Boundary Array Resulting Groups
[6, 10, 20] C6-C9, C10-C19, C20+
[6, 7, 10, 15, 20] C6, C7-C9, C10-C14, C15-C19, C20+
[7, 12, 20, 30] C7-C11, C12-C19, C20-C29, C30+

PVTlumpingModel Behavior

The PVT lumping model keeps TBP fractions (e.g., C6, C7, C8, C9) as individual pseudo-components and only lumps the characterized plus fraction (C10 through C80).

The relationship between parameters:

\[n_{\text{lumped}} = n_{\text{pseudo}} - n_{\text{TBP}}\]

where:

⚠️ Override Behavior: If the calculated numberOfLumpedComponents is less than the current value (default 7), the model overrides your setting to ensure sufficient lumping groups. A warning is logged when this occurs.

Example with 4 TBP fractions (C6-C9):

You Set Calculation Final Result
plusFractionGroups(8) Direct: 8 lumped 4 TBP + 8 lumped = 12 total
totalPseudoComponents(12) 12 - 4 = 8 lumped 4 TBP + 8 lumped = 12 total
totalPseudoComponents(5) 5 - 4 = 1 < 7 (default) → override 4 + 7 = 11 total (not 5!)

Recommendation: Use plusFractionGroups() or setNumberOfLumpedComponents() for PVTlumpingModel to avoid unexpected overrides.

Standard Lumping Model Behavior

The standard model lumps all heavy components (TBP fractions + plus fraction) into equal-weight groups:

\[w_{\text{target}} = \frac{\sum_{i=C_6}^{C_{80}} z_i \cdot M_i}{N}\]

where $N$ is the total number of pseudo-components.

Weight-Based Lumping Algorithm

SCN pseudo-components are grouped into $N$ lumps with approximately equal weight fractions:

\[w_{\text{target}} = \frac{\sum_i z_i \cdot M_i}{N}\]

For each lump $k$, the properties are averaged:

Mole Fraction

\[z_k = \sum_{i \in k} z_i\]

Molecular Weight

\(M_k = \frac{\sum_{i \in k} z_i \cdot M_i}{z_k}\)

Density (Volume-Weighted)

\(\rho_k = \frac{\sum_{i \in k} z_i \cdot M_i}{\sum_{i \in k} \frac{z_i \cdot M_i}{\rho_i}}\)

Critical Properties

Mixing rules (typically mole-fraction weighted):

\[T_{c,k} = \sum_{i \in k} \frac{z_i}{z_k} \cdot T_{c,i}\] \[P_{c,k} = \sum_{i \in k} \frac{z_i}{z_k} \cdot P_{c,i}\] \[\omega_k = \sum_{i \in k} \frac{z_i}{z_k} \cdot \omega_i\]

Usage Examples

Basic Whitson Gamma Characterization

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

// Create fluid with plus fraction
SystemInterface fluid = new SystemSrkEos(350.0, 150.0);
fluid.addComponent("methane", 0.70);
fluid.addComponent("ethane", 0.10);
fluid.addComponent("propane", 0.05);
fluid.addPlusFraction("C7+", 0.15, 150.0 / 1000.0, 0.82);  // MW in kg/mol, density in g/cm³
fluid.setMixingRule(2);

// Configure Whitson Gamma Model
fluid.getCharacterization()
    .setPlusFractionModel("Whitson Gamma Model")
    .setGammaShapeParameter(1.5)      // alpha for black oil
    .setGammaMinMW(84.0)              // eta for C7+
    .setGammaDensityModel("Soreide"); // Use Søreide density correlation

// Run characterization
fluid.getCharacterization().characterisePlusFraction();

Automatic Alpha Estimation

// Let NeqSim estimate alpha based on Watson K-factor
fluid.getCharacterization()
    .setPlusFractionModel("Whitson Gamma Model")
    .setAutoEstimateGammaAlpha(true)
    .setGammaDensityModel("Soreide");

fluid.getCharacterization().characterisePlusFraction();

// Check the estimated alpha
double alpha = ((PlusFractionModel.WhitsonGammaModel)
    fluid.getCharacterization().getPlusFractionModel()).getAlpha();
double watsonK = ((PlusFractionModel.WhitsonGammaModel)
    fluid.getCharacterization().getPlusFractionModel()).getWatsonKFactor();

System.out.println("Estimated alpha: " + alpha);
System.out.println("Watson K-factor: " + watsonK);

PVTlumpingModel - Preserve TBP Fractions

// Fluid with C6-C9 TBP fractions and C10+ plus fraction
fluid.getCharacterization().setTBPModel("PedersenSRK");
fluid.addTBPfraction("C6", 1.0, 0.086, 0.66);
fluid.addTBPfraction("C7", 2.0, 0.092, 0.73);
fluid.addTBPfraction("C8", 2.0, 0.104, 0.76);
fluid.addTBPfraction("C9", 1.0, 0.118, 0.78);
fluid.addPlusFraction("C10+", 15.0, 0.280, 0.84);

fluid.getCharacterization().setPlusFractionModel("Pedersen");

// Fluent API (Recommended): Control number of groups from C10+
fluid.getCharacterization().configureLumping()
    .model("PVTlumpingModel")
    .plusFractionGroups(5)  // C6-C9 remain separate
    .build();

fluid.getCharacterization().characterisePlusFraction();
// Result: C6_PC, C7_PC, C8_PC, C9_PC + 5 lumped groups = 9 total pseudo-components

Standard Lumping - Lump Everything from C6

// Use standard model to lump ALL heavy fractions together
fluid.getCharacterization().setPlusFractionModel("Pedersen");

// Fluent API (Recommended): Total 6 pseudo-components covering C6 through C80
fluid.getCharacterization().configureLumping()
    .model("standard")
    .totalPseudoComponents(6)
    .build();

fluid.getCharacterization().characterisePlusFraction();
// Result: PC1, PC2, PC3, PC4, PC5, PC6 covering entire C6-C80 range

Custom Boundaries - Match PVT Lab Report

// Match specific PVT lab report groupings
fluid.getCharacterization().setPlusFractionModel("Pedersen");

// Creates groups: C6, C7-C9, C10-C14, C15-C19, C20+
fluid.getCharacterization().configureLumping()
    .customBoundaries(6, 7, 10, 15, 20)
    .build();

fluid.getCharacterization().characterisePlusFraction();

Accessing Characterized Data

PlusFractionModelInterface model = fluid.getCharacterization().getPlusFractionModel();

double[] moleFractions = model.getZ();
double[] molecularWeights = model.getM();
double[] densities = model.getDens();

for (int i = model.getFirstPlusFractionNumber(); i < model.getLastPlusFractionNumber(); i++) {
    if (moleFractions[i] > 1e-10) {
        System.out.printf("SCN %d: z=%.6f, M=%.1f g/mol, rho=%.4f g/cm³%n",
            i, moleFractions[i], molecularWeights[i] * 1000, densities[i]);
    }
}

PVT Regression Framework

This section documents the PVT regression framework for automatic EOS model tuning based on experimental PVT report data. The framework is implemented in the neqsim.pvtsimulation.regression package.

Package Overview

Class Description
PVTRegression Main regression framework class
RegressionParameter Enum defining tunable parameters (BIPs, volume shifts, critical properties, CSP viscosity factors)
ExperimentType Enum for experiment types (CCE, CVD, DLE, SEPARATOR, etc.)
CCEDataPoint, CVDDataPoint, DLEDataPoint, SeparatorDataPoint, ViscosityDataPoint Data point classes for each experiment type
RegressionParameterConfig Configuration for each regression parameter with bounds
PVTRegressionFunction Objective function extending LevenbergMarquardtFunction
RegressionResult Result container with tuned fluid and uncertainty analysis
UncertaintyAnalysis Statistical uncertainty quantification

Overview

PVT regression involves adjusting equation of state parameters to minimize the deviation between calculated and experimental properties. The framework handles multiple experiment types simultaneously while maintaining physical consistency.

Multi-Objective Regression

Objective Function

The total objective function combines weighted contributions from different PVT experiments:

\[F_{obj} = \sum_{k} w_k \cdot F_k\]

where $w_k$ are user-defined weights and $F_k$ are individual experiment objective functions.

Constant Composition Expansion (CCE)

\[F_{CCE} = \sum_{i=1}^{N_{CCE}} \left[ \left(\frac{V_{rel,i}^{calc} - V_{rel,i}^{exp}}{V_{rel,i}^{exp}}\right)^2 + \lambda_Y \left(\frac{Y_i^{calc} - Y_i^{exp}}{Y_i^{exp}}\right)^2 \right]\]

where:

Key match points:

Constant Volume Depletion (CVD)

\[F_{CVD} = \sum_{i=1}^{N_{CVD}} \left[ \left(\frac{L_i^{calc} - L_i^{exp}}{L_i^{exp}}\right)^2 + \left(\frac{Z_i^{calc} - Z_i^{exp}}{Z_i^{exp}}\right)^2 + \sum_{j} \left(\frac{y_{j,i}^{calc} - y_{j,i}^{exp}}{y_{j,i}^{exp}}\right)^2 \right]\]

where:

Key match points:

Differential Liberation Expansion (DLE)

\[F_{DLE} = \sum_{i=1}^{N_{DLE}} \left[ \left(\frac{R_{s,i}^{calc} - R_{s,i}^{exp}}{R_{s,i}^{exp}}\right)^2 + \left(\frac{B_{o,i}^{calc} - B_{o,i}^{exp}}{B_{o,i}^{exp}}\right)^2 + \left(\frac{\rho_{o,i}^{calc} - \rho_{o,i}^{exp}}{\rho_{o,i}^{exp}}\right)^2 \right]\]

where:

Key match points:

Separator Test

\[F_{SEP} = \left(\frac{GOR^{calc} - GOR^{exp}}{GOR^{exp}}\right)^2 + \left(\frac{B_o^{calc} - B_o^{exp}}{B_o^{exp}}\right)^2 + \left(\frac{API^{calc} - API^{exp}}{API^{exp}}\right)^2\]

Viscosity (Optional)

\[F_{\mu} = \sum_{i} \left(\frac{\mu_i^{calc} - \mu_i^{exp}}{\mu_i^{exp}}\right)^2\]

Viscosity data are added with PVTRegression.addViscosityData(...). The framework flashes the fluid at each pressure and temperature, initializes physical properties, and compares the selected phase viscosity in Pa s. Supported phase names are gas, vapor, oil, liquid, aqueous, and water.

When the regression parameters include VISCOSITY_CSP_1 through VISCOSITY_CSP_4, NeqSim activates the PFCT/CSP viscosity model on each phase and fits the four CSP viscosity correction factors. The helper addCspViscosityRegressionParameters() registers all four factors with their default bounds.

Automatic Binary Interaction Parameter (BIP) Optimization

BIP Matrix Structure

For a system with $N_c$ components, the symmetric BIP matrix has $N_c(N_c-1)/2$ independent parameters. To reduce dimensionality, group-based BIPs are used:

Group Pairs Typical Starting $k_{ij}$ Regression Range
CH₄ - C₂-C₆ 0.00 - 0.02 Fixed or narrow
CH₄ - C7+ 0.02 - 0.05 Primary target
CO₂ - HC 0.10 - 0.15 If CO₂ present
N₂ - HC 0.04 - 0.08 If N₂ present
H₂S - HC 0.05 - 0.10 If H₂S present
C7+ - C7+ 0.00 Usually fixed

BIP Correlation (Carbon Number Based)

For C7+ pseudo-components, BIPs can be correlated:

\[k_{ij} = k_{CH_4-C7+} \cdot \left(1 - \left(\frac{2\sqrt{T_{c,i} \cdot T_{c,j}}}{T_{c,i} + T_{c,j}}\right)^n\right)\]

where $n$ is a tunable exponent (typically 0.5-2.0).

Optimization Strategy

  1. Level 1: Tune $k_{CH_4-C7+}$ to match saturation pressure
  2. Level 2: Tune $k_{C_2-C_6, C7+}$ to improve phase envelope shape
  3. Level 3: Tune correlation exponent $n$ for density matching

Volume Translation Parameter Tuning

Peneloux Volume Shift

The Peneloux (1982) volume translation corrects liquid density without affecting VLE:

\[V_{corrected} = V_{EOS} - \sum_i x_i \cdot c_i\]

where $c_i$ is the component volume shift parameter.

Regression Approach

For pseudo-components, express volume shift as:

\[c_i = c_0 + c_1 \cdot M_i + c_2 \cdot M_i^2\]

Objective: Minimize density deviation in single-phase liquid region:

\[F_c = \sum_{i} \left(\frac{\rho_i^{calc} - \rho_i^{exp}}{\rho_i^{exp}}\right)^2\]

Rackett Equation Alternative

\[Z_{RA,i} = Z_{RA}^{ref} + a \cdot (M_i - M_{ref})\]

where $Z_{RA}$ is the Rackett compressibility factor and $a$ is a tunable coefficient.

Critical Property Regression for Pseudo-Components

Parameterization

Instead of regressing individual $T_c$, $P_c$, $\omega$ values, tune the correlation coefficients:

Critical Temperature: \(T_c = (a_0 + \Delta a_0) \cdot \rho + (a_1 + \Delta a_1) \cdot \ln(M) + a_2 \cdot M + \frac{a_3}{M}\)

Critical Pressure: \(\ln(P_c) = (b_0 + \Delta b_0) + b_1 \cdot \rho^{b_4} + \frac{b_2}{M} + \frac{b_3}{M^2}\)

Acentric Factor: \(\omega = (\omega_{base} + \Delta\omega) \cdot f(M, \rho)\)

where $\Delta a_0$, $\Delta a_1$, $\Delta b_0$, $\Delta\omega$ are regression parameters.

Constraints

Physical bounds must be enforced:

\(T_c > T_b > 0\) \(P_c > 0\) \(0 < \omega < 2\) \(\frac{\partial T_c}{\partial M} > 0 \text{ (monotonic increase)}\)

Uncertainty Quantification

Parameter Sensitivity Analysis

Compute the Jacobian matrix at the optimum:

\[J_{ij} = \frac{\partial F_i}{\partial \theta_j}\]

where $F_i$ are individual residuals and $\theta_j$ are parameters.

Covariance Matrix

\[\text{Cov}(\theta) = s^2 \cdot (J^T J)^{-1}\]

where $s^2$ is the residual variance:

\[s^2 = \frac{F_{obj}}{N_{data} - N_{params}}\]

Confidence Intervals

95% confidence interval for parameter $\theta_j$:

\[\theta_j \pm t_{0.975, N-p} \cdot \sqrt{\text{Cov}(\theta)_{jj}}\]

Monte Carlo Uncertainty Propagation

  1. Sample parameters from multivariate normal: $\theta \sim N(\hat{\theta}, \text{Cov}(\theta))$
  2. Run flash calculations for each sample
  3. Compute prediction intervals for properties:
\[P_{95\%} = \left[\mu - 1.96\sigma, \mu + 1.96\sigma\right]\]

API Usage Example

import neqsim.pvtsimulation.regression.*;
import neqsim.thermo.system.SystemInterface;
import neqsim.thermo.system.SystemSrkEos;

// Create base fluid
SystemInterface fluid = new SystemSrkEos(373.15, 200.0);
fluid.addComponent("methane", 0.70);
fluid.addComponent("ethane", 0.10);
fluid.addComponent("propane", 0.05);
fluid.addComponent("n-pentane", 0.05);
fluid.addPlusFraction("C7+", 0.10, 0.150, 0.82);
fluid.setMixingRule(2);

// Create PVT regression framework
PVTRegression regression = new PVTRegression(fluid);

// Add experimental CCE data
double[] pressures = {300.0, 250.0, 200.0, 150.0, 100.0};
double[] relativeVolumes = {0.98, 1.00, 1.08, 1.25, 1.55};
regression.addCCEData(pressures, relativeVolumes, 373.15);

// Add DLE data
double[] dlePressures = {250.0, 200.0, 150.0, 100.0};
double[] rs = {150.0, 120.0, 85.0, 50.0};
double[] bo = {1.45, 1.38, 1.30, 1.20};
double[] oilDensity = {720.0, 740.0, 760.0, 780.0};
regression.addDLEData(dlePressures, rs, bo, oilDensity, 373.15);

// Configure regression parameters (with custom bounds or use defaults)
regression.addRegressionParameter(RegressionParameter.BIP_METHANE_C7PLUS, 0.0, 0.10, 0.03);
regression.addRegressionParameter(RegressionParameter.VOLUME_SHIFT_C7PLUS);  // Uses defaults

// Set weights for multi-objective optimization
regression.setExperimentWeight(ExperimentType.CCE, 1.0);
regression.setExperimentWeight(ExperimentType.DLE, 1.5);  // Prioritize DLE matching

// Configure optimization
regression.setMaxIterations(100);
regression.setVerbose(true);

// Run regression
RegressionResult result = regression.runRegression();

// Get tuned fluid
SystemInterface tunedFluid = result.getTunedFluid();

// Get optimized parameter values
double optimizedBIP = result.getOptimizedValue(RegressionParameter.BIP_METHANE_C7PLUS);
System.out.println("Optimized BIP (CH4-C7+): " + optimizedBIP);

// Uncertainty analysis
UncertaintyAnalysis uncertainty = result.getUncertainty();
double[] ci = uncertainty.getConfidenceIntervalBounds(0);
System.out.println("95% CI for BIP: [" + ci[0] + ", " + ci[1] + "]");

// Generate summary report
String summary = result.generateSummary();
System.out.println(summary);

Available Regression Parameters

Parameter Description Default Bounds
BIP_METHANE_C7PLUS BIP between methane and C7+ fractions [0.0, 0.10, 0.03]
BIP_C2C6_C7PLUS BIP between C2-C6 and C7+ fractions [0.0, 0.05, 0.01]
BIP_CO2_HC BIP between CO₂ and hydrocarbons [0.08, 0.18, 0.12]
BIP_N2_HC BIP between N₂ and hydrocarbons [0.02, 0.12, 0.05]
VOLUME_SHIFT_C7PLUS Volume shift multiplier for C7+ [0.8, 1.2, 1.0]
TC_MULTIPLIER_C7PLUS Critical temperature multiplier [0.95, 1.05, 1.0]
PC_MULTIPLIER_C7PLUS Critical pressure multiplier [0.95, 1.05, 1.0]
OMEGA_MULTIPLIER_C7PLUS Acentric factor multiplier [0.90, 1.10, 1.0]
PLUS_MOLAR_MASS_MULTIPLIER Plus fraction MW multiplier [0.90, 1.10, 1.0]
GAMMA_ALPHA Gamma distribution shape parameter [0.5, 4.0, 1.0]
GAMMA_ETA Gamma distribution minimum MW [75.0, 95.0, 84.0]
VISCOSITY_CSP_1 PFCT/CSP critical-temperature scaling correction [0.2, 2.0, 1.0]
VISCOSITY_CSP_2 PFCT/CSP critical-pressure scaling correction [0.2, 2.0, 1.0]
VISCOSITY_CSP_3 PFCT/CSP molar-mass scaling correction [0.2, 2.0, 1.0]
VISCOSITY_CSP_4 PFCT/CSP alpha correction scaling [0.2, 2.0, 1.0]

Data Input Format (CSV)

CCE Data:

Pressure(bara),RelativeVolume,YFactor
350.0,0.9850,
300.0,0.9912,
250.0,1.0000,  # Saturation point
200.0,1.0523,1.0234
150.0,1.1876,1.0456

DLE Data:

Pressure(bara),Rs(Sm3/Sm3),Bo(m3/Sm3),OilDensity(kg/m3),GasGravity
250.0,150.5,1.425,725.3,0.85
200.0,120.2,1.380,742.1,0.82
150.0,85.6,1.312,761.5,0.79

Viscosity Data:

Pressure(bara),Temperature(K),Viscosity(Pa s),Phase
1.0,320.0,0.000312,oil
5.0,320.0,0.000331,oil
15.0,320.0,0.000370,oil

Implementation Status

Phase Feature Status
1 CCE/DLE simulation with objective function ✅ Implemented
2 BIP regression for saturation pressure ✅ Implemented
3 Volume translation optimization ✅ Implemented
4 CVD simulation and regression ✅ Implemented
5 Critical property correlation tuning ✅ Implemented
6 Multi-objective optimization framework ✅ Implemented
7 Uncertainty quantification ✅ Implemented
8 Viscosity data and CSP/PFCT parameter regression ✅ Implemented
9 GUI/Report generation 🔲 Future work

Common Fluid Characterization (Matching Pseudo-Components)

When working with multiple reservoir fluids in a simulation model (e.g., compositional reservoir simulation, commingled production), all fluids must share the same pseudo-component (PC) structure. NeqSim provides utilities for this workflow based on Pedersen et al. (Chapter 5.5-5.6).

PseudoComponentCombiner Utility

The PseudoComponentCombiner class provides three complementary workflows, all based on Pedersen et al. (Chapter 5.5–5.6):

import neqsim.thermo.characterization.PseudoComponentCombiner;

// 1. Match a source fluid to an existing reference PC structure (Chapter 5.6,
//    "Common EoS" slate). The source inherits the reference lump properties.
SystemInterface matched = PseudoComponentCombiner.characterizeToReference(
    sourceFluid, referenceFluid);

// 2. Blend several fluids into ONE merged fluid sharing a common PC structure
//    (Chapter 5.5).
SystemInterface combined = PseudoComponentCombiner.combineReservoirFluids(
    6,                       // target number of shared pseudo-components
    fluid1, fluid2, fluid3);

// 3. Re-characterize several fluids onto a SHARED common slate while keeping
//    them as SEPARATE fluids (Chapter 5.6, Eqs. 5.55–5.60). Boundaries come
//    from an equal-mass cut of the weighted "imaginary" composition.
List<SystemInterface> commonSlate = PseudoComponentCombiner.characterizeToCommonSlate(
    Arrays.asList(fluidA, fluidB),
    new double[] {0.5, 0.5});  // per-fluid weights
Method Pedersen Output Use case
characterizeToReference(src, ref[, opts]) Ch. 5.6 one fluid on the reference slate match a field to a trusted master fluid
combineReservoirFluids(n, fluids…) Ch. 5.5 one merged fluid commingled / blended stream
characterizeToCommonSlate(fluids, weights[, n]) Ch. 5.6 (Eqs. 5.55–5.60) a list of fluids on a shared slate keep fields separate but EoS-compatible

CharacterizationOptions

For advanced control, use the CharacterizationOptions builder:

import neqsim.thermo.characterization.CharacterizationOptions;
import neqsim.thermo.characterization.CharacterizationOptions.NamingScheme;

CharacterizationOptions options = CharacterizationOptions.builder()
    .transferBinaryInteractionParameters(true)  // Copy BIPs from reference
    .normalizeComposition(true)                 // Ensure mole fractions sum to 1.0
    .namingScheme(NamingScheme.REFERENCE)       // Use reference component names
    .generateValidationReport(true)             // Create before/after comparison
    .inheritReferenceProperties(true)           // Inherit reference lump properties
    .delumpBeforeRecharacterization(false)      // Split source lumps before re-binning
    .delumpResolution(12)                        // Sub-fractions per source lump
    .sharedImaginaryBoundaries(false)           // Equal-mass reference cut points
    .build();

SystemInterface matched = PseudoComponentCombiner.characterizeToReference(
    sourceFluid, referenceFluid, options);
Option Description Default
transferBinaryInteractionParameters Copy BIPs from reference fluid false
normalizeComposition Normalize mole fractions to sum to 1.0 false
namingScheme Use SOURCE, REFERENCE, or MERGED names REFERENCE
generateValidationReport Generate validation report false
inheritReferenceProperties Inherit the reference lump properties (molar mass, density, critical constants). When false, lump properties are recomputed from the source mass on the reference cut grid true
delumpBeforeRecharacterization Split each coarse source lump into finer single-carbon-number sub-fractions before re-binning onto the reference cuts (Pedersen Ch. 5 delumping). Conserves source moles and mass exactly false
delumpResolution Number of sub-fractions per source lump when delumping (also the fine grid for equal-mass reference boundaries). Values ≤ 1 disable splitting 12
sharedImaginaryBoundaries Place reference cut boundaries as carbon-number equal-mass cut points on the reference’s imaginary (delumped) composition (Pedersen Ch. 5.6, Eqs. 5.58–5.59) instead of boiling-point midpoints false

Property inheritance vs. recomputation (inheritReferenceProperties)

The Common EoS slate of Pedersen Ch. 5.6 requires every fluid characterized to the same reference to share an identical pseudo-component property set, differing only in the mole fractions. With inheritReferenceProperties(true) (default) each output lump copies the reference molar mass, density, and critical constants — so all fluids are directly EoS-compatible. Set it to false to keep grid-only behaviour: the cut boundaries still come from the reference, but each lump’s molar mass and density are recomputed from the source mass landing in that cut.

Delumping before re-characterization (delumpBeforeRecharacterization)

When a field’s native lumps already sit close to the reference grid, the naive source→reference mapping is effectively the identity: lump mole fractions are frozen and per-cut mass is not conserved against the reference molar masses. Enabling delumpBeforeRecharacterization(true) first splits each parent lump into delumpResolution single-carbon-number sub-fractions whose moles and mass exactly reproduce the parent (a single linear molar-mass rescale enforces $\sum_k n_k M_k = n_\text{parent} M_\text{parent}$). The boiling point spreads monotonically with molar mass so sub-fractions can cross reference cut boundaries; density and critical constants stay at the parent values. The sub-fractions are then re-lumped onto the reference grid, so each cut’s molar mass = mass / moles is recomputed self-consistently. This is most effective with inheritReferenceProperties(false); combining it with inheritReferenceProperties(true) logs a warning because the reference values still overwrite the redistributed lump properties.

Equal-mass reference boundaries (sharedImaginaryBoundaries)

By default the reference cut boundaries are the arithmetic boiling-point midpoints between adjacent reference lumps, which ignores how much mass each lump represents. With sharedImaginaryBoundaries(true) the reference is first rebuilt into a fine single-carbon-number “imaginary” composition (using delumpResolution) and the cut points are placed so each cut carries an equal mass fraction — the reference-only (NFLUID = 1) form of the Pedersen Ch. 5.6 common-slate rule (Eqs. 5.58–5.59 with §5.3 cut criterion). Each equal-mass cut is then clamped into the gap between the two adjacent reference pseudo-components, guaranteeing every reference lump stays inside its own cut so the one-to-one property inheritance is preserved even when the reference lumps carry unequal mass. With delumpResolution ≤ 1 it falls back to boiling-point midpoints.

Why reference-only here? characterizeToReference inherits the trusted reference’s lump properties one-to-one, so the reference fluid is the correct authority for the grid. The pooled multi-fluid imaginary composition of Eqs. 5.58–5.59 remains used by characterizeToCommonSlate (free new slate); applying it to the inherit-from-fixed-reference path would break the inherit alignment.

Common Slate for Multiple Fluids (characterizeToCommonSlate)

characterizeToCommonSlate implements the Pedersen et al. (Chapter 5.6, Eqs. 5.55–5.60) common-slate procedure for keeping several fluids as separate systems while forcing them onto an identical pseudo-component slate (so they are mutually EoS-compatible without merging compositions).

import neqsim.thermo.characterization.PseudoComponentCombiner;

List<SystemInterface> aligned = PseudoComponentCombiner.characterizeToCommonSlate(
    Arrays.asList(fluidA, fluidB, fluidC),
    new double[] {0.5, 0.3, 0.2});   // per-fluid weights (need not sum to 1)

// Optional: set the shared lump count explicitly instead of inferring it
List<SystemInterface> aligned6 = PseudoComponentCombiner.characterizeToCommonSlate(
    Arrays.asList(fluidA, fluidB), new double[] {0.5, 0.5}, 6);

Procedure:

  1. Imaginary composition — each input fluid is normalized to unit pseudo-component mass and pooled with its weight $w_f$ into a single weighted “imaginary” composition (Eqs. 5.55–5.57). This composition is not a real fluid; it only sets the cut grid, so no single fluid is privileged.

    \[\bar{m}(T_b) = \sum_f w_f \, \frac{m_f(T_b)}{\sum_{T_b} m_f(T_b)}\]

  2. Equal-mass cut points — the imaginary composition is lumped into $N$ pseudo-components using equal-mass cut points on the carbon-number axis (Eqs. 5.58–5.59, §5.3 criterion): each cut $i$ carries $1/N$ of the total imaginary mass.

    \[\sum_{T_b \le T_b^{(i)}} \bar{m}(T_b) = \frac{i}{N} \sum_{T_b} \bar{m}(T_b), \qquad i = 1 \dots N-1\]

  3. Apply to every fluid — each input fluid is independently binned onto the shared boundaries and its lump properties are computed from its own mass in each cut (Eq. 5.60). All returned fluids therefore share the same $N$ lump boiling-point ranges and remain individually mass-conserving.

If the shared lump count is not supplied it is inferred from the input fluids (inferCommonSlateSize). Weights let a more representative or higher-rate fluid pull the cut grid toward its own carbon-number distribution.

BIP Transfer

Binary Interaction Parameters (BIPs) can be transferred between fluids:

// Transfer BIPs during characterization
PseudoComponentCombiner.transferBinaryInteractionParameters(
    sourceFluid, referenceFluid);

// Fluent API on Characterise class
SystemInterface fluid = new SystemSrkEos(298, 50);
fluid.addComponent("methane", 0.7);
fluid.addPlusFraction("C7+", 0.3, 0.200, 0.85);
fluid.getCharacterization()
    .setTBPModel("PedersenSRK")
    .characterize()
    .transferBipsFrom(tunedReferenceFluid);

Validation Reports

The CharacterizationValidationReport provides before/after comparison:

CharacterizationValidationReport report =
    PseudoComponentCombiner.generateValidationReport(sourceFluid, matchedFluid);

System.out.println("Mass conserved: " + report.isMassConserved());
System.out.println("Moles conserved: " + report.isMolesConserved());
System.out.println("PC count before: " + report.getSourcePseudoComponentCount());
System.out.println("PC count after: " + report.getResultPseudoComponentCount());
System.out.println(report.toReportString());

Mathematical Background

When matching a source fluid to a reference PC structure:

  1. Component Mapping: Discrete components (C1, C2, CO₂, etc.) are mapped directly
  2. Cut boundaries: The reference lumps define cut boundaries on the carbon-number (boiling-point) axis — boiling-point midpoints by default, or carbon-number equal-mass cut points with sharedImaginaryBoundaries(true)
  3. PC binning: Each source pseudo-component (optionally first delumped into single-carbon-number sub-fractions, see delumpBeforeRecharacterization) is assigned to the reference cut whose boundary range contains its boiling point; the lump moles and mass are accumulated per cut
  4. Lump properties: Either inherited from the reference (inheritReferenceProperties(true)) or recomputed as molar mass = mass / moles per cut (inheritReferenceProperties(false))
  5. Mass Conservation: Total source pseudo-fraction moles and mass are preserved through the binning
  6. BIP Transfer: For EOS phases, BIPs are copied element-by-element from the reference

References

  1. Whitson, C.H. (1983). “Characterizing Hydrocarbon Plus Fractions.” SPE Journal, 23(4), 683-694. SPE-12233-PA.

  2. Pedersen, K.S., Thomassen, P., and Fredenslund, A. (1984). “Thermodynamics of Petroleum Mixtures Containing Heavy Hydrocarbons. 1. Phase Envelope Calculations by Use of the Soave-Redlich-Kwong Equation of State.” Industrial & Engineering Chemistry Process Design and Development, 23(1), 163-170.

  3. Pedersen, K.S., Christensen, P.L., and Shaikh, J.A. (2014). “Phase Behavior of Petroleum Reservoir Fluids,” 2nd ed., Chapter 5 (lumping/delumping, §5.5 fluid combination, §5.6 common-slate Eqs. 5.55–5.60). CRC Press.

  4. Søreide, I. (1989). “Improved Phase Behavior Predictions of Petroleum Reservoir Fluids from a Cubic Equation of State.” Dr.Ing. Thesis, Norwegian Institute of Technology (NTH), Trondheim.

  5. Riazi, M.R. and Daubert, T.E. (1980). “Simplify Property Predictions.” Hydrocarbon Processing, 59(3), 115-116.

  6. Kesler, M.G. and Lee, B.I. (1976). “Improve Prediction of Enthalpy of Fractions.” Hydrocarbon Processing, 55(3), 153-158.

  7. Whitson, C.H. and Brulé, M.R. (2000). “Phase Behavior.” SPE Monograph Series, Vol. 20. Society of Petroleum Engineers.

Appendix A: Unit Conventions in NeqSim

Property Internal Unit Common Input Unit
Molecular weight kg/mol g/mol (÷1000)
Density g/cm³ g/cm³ or kg/m³ (÷1000)
Temperature K °C (+273.15)
Pressure bar bara
Critical temperature K K
Critical pressure bar bar

Note: When using addPlusFraction(), molecular weight should be in kg/mol and density in g/cm³ (specific gravity).

Appendix B: Mathematical Notation Summary

Symbol Description Typical Range
$\alpha$ Gamma shape parameter 0.5 - 4.0
$\beta$ Gamma scale parameter Calculated
$\eta$ Minimum molecular weight 84 - 90 g/mol
$M$ Molecular weight g/mol
$\rho$ Density (specific gravity) 0.6 - 1.0 g/cm³
$K_w$ Watson characterization factor 10 - 13
$T_c$ Critical temperature K
$P_c$ Critical pressure bar
$\omega$ Acentric factor 0.2 - 1.5
$T_b$ Normal boiling point K
$z$ Mole fraction 0 - 1