Class PhaseSrkCPAfullyImplicit

java.lang.Object
neqsim.thermo.phase.Phase
neqsim.thermo.phase.PhaseEos
neqsim.thermo.phase.PhaseSrkEos
neqsim.thermo.phase.PhaseSrkCPA
neqsim.thermo.phase.PhaseSrkCPAs
neqsim.thermo.phase.PhaseSrkCPAfullyImplicit
All Implemented Interfaces:
Serializable, Cloneable, PhaseCPAInterface, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface
public class PhaseSrkCPAfullyImplicit extends PhaseSrkCPAs
Fully implicit CPA phase implementation based on Igben et al. (2026).

Implements the algorithm from: "Fully implicit algorithm for Cubic Plus Association equation of state", Fluid Phase Equilibria 608 (2026) 114734.

Key differences from the standard nested approach in PhaseSrkCPA:

  • Solves molar volume and association site fractions simultaneously using coupled Newton-Raphson
  • Eliminates inner iterations for XA solving during volume calculation
  • Includes restart criterion to suppress unnecessary root searches in supercritical regions
  • Reports 30-80% reduction in computational cost
Version:
1.0
Author:
Even Solbraa
See Also:

  • Field Details

    • serialVersionUID

      private static final long serialVersionUID
      Serialization version UID.
      See Also:

    • logger

      private static final org.apache.logging.log4j.Logger logger
      Logger object for class.

    • MAX_IMPLICIT_ITERATIONS

      private static final int MAX_IMPLICIT_ITERATIONS
      Maximum iterations for the fully implicit solver.
      See Also:

    • CONVERGENCE_TOL

      private static final double CONVERGENCE_TOL
      Convergence tolerance for the coupled system.
      See Also:

    • MAX_REL_STEP

      private static final double MAX_REL_STEP
      Maximum relative step size per iteration to prevent divergence.
      See Also:

    • RESTART_ALPHA

      private static final double RESTART_ALPHA
      Restart threshold parameter alpha (supercritical detection).
      See Also:

    • implicitCallCount

      private static volatile long implicitCallCount
      Total molarVolume calls via implicit solver.

    • totalImplicitIters

      private static volatile long totalImplicitIters
      Total coupled Newton iterations across all calls.

    • fallbackCount

      private static volatile long fallbackCount
      Total fallbacks to nested solver.

    • workKlk

      private transient double[][] workKlk

    • workHess

      private transient double[][] workHess

    • workKsi

      private transient double[] workKsi

    • workM

      private transient double[] workM

    • workKlkKsi

      private transient double[] workKlkKsi

    • workXV

      private transient double[] workXV

    • workNs

      private transient int workNs

  • Constructor Details

    • PhaseSrkCPAfullyImplicit

      public PhaseSrkCPAfullyImplicit()
      Constructor for PhaseSrkCPAfullyImplicit.

  • Method Details

    • resetProfileCounters

      public static void resetProfileCounters()
      Reset profiling counters.

    • getProfileSummary

      public static String getProfileSummary()
      Get profiling summary string.
      Returns:
      a summary of profiling data

    • clone

      public PhaseSrkCPAfullyImplicit clone()

      clone.

      Specified by:
      clone in interface PhaseInterface
      Overrides:
      clone in class PhaseSrkCPAs
      Returns:
      a PhaseInterface object

    • addComponent

      public void addComponent(String name, double moles, double molesInPhase, int compNumber)

      Add component to component array and update moles variables.

      Specified by:
      addComponent in interface PhaseInterface
      Overrides:
      addComponent in class PhaseSrkCPAs
      Parameters:
      name - Name of component.
      moles - Total number of moles of component.
      molesInPhase - Number of moles in phase.
      compNumber - Index number of component in phase object component array.

    • molarVolume

      public double molarVolume(double pressure, double temperature, double A, double B, PhaseType pt) throws IsNaNException, TooManyIterationsException

      molarVolume.

      Fully implicit molar volume calculation. Solves the normalized density zeta = b/v and association site fractions X_k simultaneously using a coupled Newton-Raphson method.

      Per Igben et al. (2026), the key speedup is eliminating inner XA iterations during the volume loop. The CPA pressure derivatives (dFCPAdV, dFCPAdVdV) are computed directly from the current X_k values using O(ns^2) sums, avoiding costly matrix inversions.

      Specified by:
      molarVolume in interface PhaseInterface
      Overrides:
      molarVolume in class PhaseSrkCPA
      Parameters:
      pressure - a double
      temperature - a double
      A - a double
      B - a double
      pt - the PhaseType of the phase
      Returns:
      a double
      Throws:
      IsNaNException - if any.
      TooManyIterationsException - if any.

    • ensureWorkArrays

      private void ensureWorkArrays(int ns)
      Ensure work arrays are allocated for the given association site count.
      Parameters:
      ns - number of association sites

    • initCPAMatrix

      public void initCPAMatrix(int type)
      Override CPA matrix initialization.

      For type == 1 (volume derivatives), computes FCPA, dFCPAdV, dFCPAdVdV, dFCPAdVdVdV from scratch using Gaussian elimination to solve H*XV = KlkV*ksi directly. This avoids the expensive EJML hessianMatrix.invert() in the parent's solveX().

      Key insight: KlkV[i][j] = fV * Klk[i][j] where fV = gcpav - 1/V. Similarly for KlkVV and KlkVVV. This means only Klk needs to be stored; volume derivatives are scalar multiples.

      For type >= 2, calls solveX() to populate hessianInvers and delegates to super.initCPAMatrix(type).

      Overrides:
      initCPAMatrix in class PhaseSrkCPA
      Parameters:
      type - 1 for volume derivatives, 2+ for temperature/composition derivatives

    • molarVolumeChangePhase

      public double molarVolumeChangePhase(double pressure, double temperature, double A, double B, PhaseType pt) throws IsNaNException, TooManyIterationsException

      molarVolumeChangePhase.

      Overrides:
      molarVolumeChangePhase in class PhaseSrkCPA
      Parameters:
      pressure - a double
      temperature - a double
      A - a double
      B - a double
      pt - the PhaseType of the phase
      Returns:
      a double
      Throws:
      IsNaNException - if any.
      TooManyIterationsException - if any.

    • setXsiteOnComponents

      private void setXsiteOnComponents(double[] xSite)
      Set XA site fraction values on component objects from flat array.
      Parameters:
      xSite - array of site fractions, indexed by total site number

    • readXsiteFromComponents

      private void readXsiteFromComponents(double[] xSite)
      Read XA site fraction values from component objects into flat array.
      Parameters:
      xSite - array to fill with site fractions

    • getDeltaij

      private double getDeltaij(int i, int j)
      Get delta[i][j] from the parent's delta array. This array is populated by calcDelta() and updateDeltaWithG().
      Parameters:
      i - site index i
      j - site index j
      Returns:
      delta value

    • getDeltaNogij

      private double getDeltaNogij(int i, int j)
      Get deltaNog[i][j] from the parent's deltaNog array.
      Parameters:
      i - site index i
      j - site index j
      Returns:
      deltaNog value

    • updateDeltaWithG

      private void updateDeltaWithG()
      Update delta = deltaNog * g where g is the radial distribution function.

    • solveLinearSystem

      private static boolean solveLinearSystem(double[][] a, double[] b, int n)
      Solve a linear system Ax = b in-place using Gaussian elimination with partial pivoting. Solution is stored in b.
      Parameters:
      a - coefficient matrix (modified in place)
      b - right-hand side (solution on return)
      n - system dimension
      Returns:
      true if solve succeeded, false if singular