Class PhaseSAFTVRMie

java.lang.Object
neqsim.thermo.phase.Phase
neqsim.thermo.phase.PhaseEos
neqsim.thermo.phase.PhaseSrkEos
neqsim.thermo.phase.PhaseSAFTVRMie
All Implemented Interfaces:
Serializable, Cloneable, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface
public class PhaseSAFTVRMie extends PhaseSrkEos
Phase class for the SAFT-VR Mie equation of state following Lafitte et al. (2013).

Implements the Helmholtz free energy as: A = A_ideal + A_mono + A_chain where A_mono comprises hard-sphere reference plus first-order (a1) and second-order (a2) Barker-Henderson perturbation terms for the Mie potential, and A_chain is the chain-connectivity correction. The perturbation terms use the Sutherland mean-field energy a1S with effective packing fraction mapping and the B correction, following the exact formulation of Lafitte et al. J. Chem. Phys. 139, 154504 (2013).

Version:
$Id: $Id
Author:
Even Solbraa
See Also:

  • Field Details

    • serialVersionUID

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

    • logger

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

    • cachedMolarVolume

      private transient double cachedMolarVolume
      Cached molar volume for successive initial guess.

    • gMieCorrectionEnabled

      private transient boolean gMieCorrectionEnabled
      Controls whether g_Mie perturbation corrections are used in chain term.

    • gMieBlendFraction

      private transient double gMieBlendFraction
      Blending fraction for g_Mie perturbation correction in chain term. 0.0 = g_HS(x0) only, 1.0 = full g_Mie. Used by the continuation/homotopy solver.

    • useHS

      int useHS
      Toggle hard-sphere plus chain contribution (1 = on, 0 = off).

    • useDISP

      int useDISP
      Toggle dispersion contribution (1 = on, 0 = off).

    • volumeSAFT

      double volumeSAFT

    • nSAFT

      double nSAFT

    • dSAFT

      double dSAFT

    • dmeanSAFT

      double dmeanSAFT

    • mSAFT

      double mSAFT

    • mmin1SAFT

      double mmin1SAFT

    • ghsSAFT

      double ghsSAFT

    • aHSSAFT

      double aHSSAFT

    • dnSAFTdV

      double dnSAFTdV

    • dnSAFTdVdV

      double dnSAFTdVdV

    • daHSSAFTdN

      double daHSSAFTdN

    • daHSSAFTdNdN

      double daHSSAFTdNdN

    • dgHSSAFTdN

      double dgHSSAFTdN

    • dgHSSAFTdNdN

      double dgHSSAFTdNdN

    • dNSAFTdT

      double dNSAFTdT

    • dNSAFTdTdT

      double dNSAFTdTdT

    • dNSAFTdTdV

      double dNSAFTdTdV

    • a1Disp

      double a1Disp

    • da1DispDeta

      double da1DispDeta

    • d2a1DispDeta2

      double d2a1DispDeta2

    • da1DispDT

      double da1DispDT

    • d2a1DispDT2

      double d2a1DispDT2

    • d2a1DispDetaDT

      double d2a1DispDetaDT

    • a2Disp

      double a2Disp

    • da2DispDeta

      double da2DispDeta

    • d2a2DispDeta2

      double d2a2DispDeta2

    • da2DispDT

      double da2DispDT

    • d2a2DispDT2

      double d2a2DispDT2

    • d2a2DispDetaDT

      double d2a2DispDetaDT

    • a3Disp

      double a3Disp

    • da3DispDeta

      double da3DispDeta

    • d2a3DispDeta2

      double d2a3DispDeta2

    • da3DispDT

      double da3DispDT

    • d2a3DispDT2

      double d2a3DispDT2

    • d2a3DispDetaDT

      double d2a3DispDetaDT

    • aDispPerComp

      double[] aDispPerComp
      Per-component weighted pair-sum of dispersion terms: aDispPerComp[i] = sum_l xs_l * (a1_il+a2_il+a3_il). Used for analytical dF_DISP/dNi computation. Only allocated for multi-component systems.

    • useASSOC

      int useASSOC
      Toggle association contribution (1 = on, 0 = off).

    • totalNumberOfAssociationSites

      int totalNumberOfAssociationSites
      Total number of association sites across all components.

    • siteOffset

      int[] siteOffset
      Site offset for component i: first site of comp i is siteOffset[i].

    • selfAssocScheme

      int[][][] selfAssocScheme
      Self-association scheme indicator. selfAssocScheme[comp][siteA][siteB] = 1 if site A can bond with site B on the same component (0 otherwise).

    • crossAssocScheme

      int[][][][] crossAssocScheme
      Cross-association scheme indicator. crossAssocScheme[compI][compJ][siteA][siteB] = 1 if site A on comp I can bond with site B on comp J.

    • deltaAssoc

      double[][] deltaAssoc
      Association strength delta[siteI][siteJ] between flattened sites (global indexing). Delta = scheme * (exp(eps_HB/RT) - 1) * sigma_ij^3 * NA * 1e5 * kappa_ij * g_HS.

    • deltadTAssoc

      double[][] deltadTAssoc
      Temperature derivative of delta.

    • hcpatot

      double hcpatot
      Sum of (1-XA)*ni over all sites: hcpatot = sum_i ni * sum_A (1 - X_Ai).

    • hcpatotdT

      double hcpatotdT
      dh/dT.

    • hcpatotdTdT

      double hcpatotdTdT
      d2h/dT2.

    • gcpaAssoc

      double gcpaAssoc
      Association g (hard-sphere RDF used in delta). For SAFT-VR Mie this equals ghsSAFT.

    • gcpavAssoc

      double gcpavAssoc
      d(ln g)/dV in SI units (m^3).

    • gcpavvAssoc

      double gcpavvAssoc
      d2(ln g)/dV2 in SI units (m^6).

    • gcpavvvAssoc

      double gcpavvvAssoc
      d3(ln g)/dV3 in SI units.

    • cachedFAssoc

      double cachedFAssoc
      Cached F_ASSOC value from volInit.

    • cacheddFAssocDV

      double cacheddFAssocDV
      Cached dF_ASSOC/dV_SI from volInit (m^-3).

    • cacheddFAssocDVDV

      double cacheddFAssocDVDV
      Cached d2F_ASSOC/dV_SI2 from volInit (m^-6), computed numerically.

    • assocDVDVValid

      boolean assocDVDVValid
      Whether cacheddFAssocDVDV is valid (computed after molar volume convergence).

    • etaEffCoeffs

      static final double[][] etaEffCoeffs
      Effective packing fraction parameterization coefficients from Lafitte 2013. c_i(lambda) = A[i][0] + A[i][1]/(lambda-3) + A[i][2]/(lambda-3)^2 + A[i][3]/(lambda-3)^3 Rows: c1, c2, c3, c4.

    • phiPade

      static final double[][] phiPade
      Padé coefficient matrix for f1-f6 functions (Lafitte 2013 Table 3).

      Each column m (0-5) corresponds to function f_{m+1}. Rows 0-3: numerator coefficients (alpha^0 to alpha^3). Rows 4-6: denominator coefficients (alpha^1 to alpha^3; alpha^0 = 1). f_m(alpha) = (phi[0][m] + phi[1][m]*alpha + phi[2][m]*alpha^2 + phi[3][m]*alpha^3) / (1 + phi[4][m]*alpha + phi[5][m]*alpha^2 + phi[6][m]*alpha^3).

    • DUFAL_C

      static final double[][] DUFAL_C
      Dufal 2015 association integral I(T*, rho*) coefficient matrix.

      From Dufal et al. (2015) Mol. Phys. 113(9-10), 948-984, Table A-1. I(Tr, rhoStar) = sum_n sum_m c[n][m] * Tr^m * rhoStar^n, where Tr = T / (epsilon/kB) and rhoStar = rhoS * sigma^3. The matrix is upper-triangular: for row n, m runs from 0 to 10-n. Row index n = power of rhoStar (0..10). Column index m = power of Tr (0..10).

  • Constructor Details

    • PhaseSAFTVRMie

      public PhaseSAFTVRMie()
      Constructor for PhaseSAFTVRMie.

  • Method Details

    • calcDufalI

      public static double calcDufalI(double Tr, double rhoStar)
      Evaluates the Dufal 2015 association integral I(Tr, rhoStar).

      I = sum_{n=0}^{10} sum_{m=0}^{10-n} c[n][m] * Tr^m * rhoStar^n

      Parameters:
      Tr - reduced temperature T / (epsilon/kB)
      rhoStar - reduced segment density rhoS * sigma^3
      Returns:
      value of the I integral

    • calcDufalIdRhoStar

      public static double calcDufalIdRhoStar(double Tr, double rhoStar)
      Evaluates the partial derivative dI/d(rhoStar) of the Dufal 2015 association integral.
      Parameters:
      Tr - reduced temperature T / (epsilon/kB)
      rhoStar - reduced segment density rhoS * sigma^3
      Returns:
      dI/d(rhoStar)

    • calcDufalIdTr

      public static double calcDufalIdTr(double Tr, double rhoStar)
      Evaluates the partial derivative dI/dTr of the Dufal 2015 association integral.
      Parameters:
      Tr - reduced temperature T / (epsilon/kB)
      rhoStar - reduced segment density rhoS * sigma^3
      Returns:
      dI/dTr

    • clone

      public PhaseSAFTVRMie clone()

      clone.

      Specified by:
      clone in interface PhaseInterface
      Overrides:
      clone in class PhaseSrkEos
      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 PhaseSrkEos
      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.

    • init

      public void init(double totalNumberOfMoles, int numberOfComponents, int initType, PhaseType pt, double beta)

      init.

      initType used in component.init()

      Calls component.init(initType)

      Specified by:
      init in interface PhaseInterface
      Overrides:
      init in class PhaseEos
      Parameters:
      totalNumberOfMoles - Total number of moles in all phases of Stream.
      numberOfComponents - Number of components in system.
      initType - a int. Use 0 to init, and 1 to reset.
      pt - Type of phase.
      beta - Mole fraction of this phase in system.

    • initAssociationSchemes

      private void initAssociationSchemes(int numberOfComponents)
      Detects associating components and builds the site-site bonding scheme arrays. Called once at initType=0. Sets useASSOC=1 if any component has association sites.
      Parameters:
      numberOfComponents - number of components

    • hasAssociationParams

      private boolean hasAssociationParams(int compIndex)
      Checks whether component i has valid SAFT-VR Mie association parameters (both VR Mie segment params and association energy/volume must be set).
      Parameters:
      compIndex - component index
      Returns:
      true if component has association params

    • setupAssociationScheme

      private void setupAssociationScheme(int[][] scheme, String schemeName, int nSites)
      Sets up self-association scheme matrix for a given scheme type. Following the standard CPA/SAFT convention: sites are split into electron donors and acceptors.
      Parameters:
      scheme - output matrix [nSites][nSites] to fill with 0/1
      schemeName - scheme name from database: "4C", "2B", "3B", "1A", "2A"
      nSites - number of sites

    • setupCrossAssociationScheme

      private void setupCrossAssociationScheme(int[][] scheme, String schemeI, int nSitesI, String schemeJ, int nSitesJ)
      Sets up cross-association scheme between two components using CR-1 combining rule. Donors on one component bond with acceptors on the other.
      Parameters:
      scheme - output matrix [nSitesI][nSitesJ]
      schemeI - scheme name of component I
      nSitesI - number of sites on I
      schemeJ - scheme name of component J
      nSitesJ - number of sites on J

    • getSiteDonors

      private boolean[] getSiteDonors(String schemeName, int nSites)
      Returns donor mask for sites of a given scheme type.
      Parameters:
      schemeName - scheme name
      nSites - number of sites
      Returns:
      boolean array where true = donor

    • getSiteAcceptors

      private boolean[] getSiteAcceptors(String schemeName, int nSites)
      Returns acceptor mask for sites of a given scheme type.
      Parameters:
      schemeName - scheme name
      nSites - number of sites
      Returns:
      boolean array where true = acceptor

    • volInit

      public void volInit()
      Initializes all SAFT variables for the current molar volume. Called at each iteration of the volume solver and after volume is converged.

    • recomputeSAFTBaseQuantities

      private void recomputeSAFTBaseQuantities()
      Recomputes the base SAFT quantities (nSAFT, dnSAFTdV, etc.) for the current molarVolume. This is used during numerical association derivative computation to update the SAFT state without running a full volInit (which would recurse).

    • calcLnGChainEffective

      private double calcLnGChainEffective(double etaVal, double alpha)
      Computes the effective chain ln(g) as a composition-weighted average over components. ln(g_eff) = sum_i [w_i * ln(g_Mie_ii)] where w_i = x_i*(m_i-1) / sum_j x_j*(m_j-1). For a single-component system this reduces to ln(g_Mie_ii). For mixtures, components with m_i = 1 do not contribute (their chain weight is zero).
      Parameters:
      etaVal - packing fraction
      alpha - g_Mie blend fraction (0 = g_HS only, 1 = full g_Mie)
      Returns:
      weighted average ln(g)

    • calcEtaEffCoeff

      static double calcEtaEffCoeff(int i, double lambda)
      Effective packing fraction coefficient c_i(lambda) (Lafitte 2013 Table 5). c_i(lambda) = A[i][0] + A[i][1]/lambda + A[i][2]/lambda^2 + A[i][3]/lambda^3.
      Parameters:
      i - coefficient index (0-3 for c1-c4)
      lambda - Mie exponent
      Returns:
      c_i(lambda)

    • calcEtaEff

      static double calcEtaEff(double eta, double lambda)
      Effective packing fraction eta_eff = c1*eta + c2*eta^2 + c3*eta^3 + c4*eta^4.
      Parameters:
      eta - actual packing fraction
      lambda - Mie exponent
      Returns:
      eta_eff

    • calcDEtaEffDeta

      static double calcDEtaEffDeta(double eta, double lambda)
      d(eta_eff)/d(eta).
      Parameters:
      eta - actual packing fraction
      lambda - Mie exponent
      Returns:
      derivative

    • calcA1Sutherland

      static double calcA1Sutherland(double eta, double lambda, double epsOverKT)
      First-order Sutherland mean-field energy a1S/(kT) (Lafitte 2013 Eq. 25).
      Parameters:
      eta - packing fraction
      lambda - Mie exponent
      epsOverKT - epsilon/(kT)
      Returns:
      a1S_bar (dimensionless)

    • calcBCorrection

      static double calcBCorrection(double eta, double lambda, double epsOverKT, double x0)
      B correction term (Lafitte 2013 Eq. 36) using I/J integrals and gHS contact value. B = 12 * eta * (eps/kT) * [gHS(eta)*I(x0,lambda) - 9*eta*(1+eta)/(2*(1-eta)^3)*J(x0,lambda)] where I = (1 - x0^(3-lambda))/(lambda-3) and J = (1 - (lambda-3)*x0^(4-lambda) + (lambda-4)*x0^(3-lambda)) / ((lambda-3)*(lambda-4)).
      Parameters:
      eta - packing fraction
      lambda - Mie exponent
      epsOverKT - epsilon/(kT)
      x0 - sigma/d ratio
      Returns:
      B (dimensionless, including 12*eta*eps/kT factor)

    • calcA1MieAtEta

      public static double calcA1MieAtEta(double eta, double lambdaR, double lambdaA, double epsOverKT, double cMie, double x0)
      Full first-order Mie perturbation a1/(NkT) at a given eta value (Lafitte 2013 Eq. 40).
      Parameters:
      eta - packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      cMie - Mie prefactor
      x0 - sigma/d ratio
      Returns:
      a1_Mie (dimensionless)

    • calcKHS

      public static double calcKHS(double eta)
      Hard-sphere isothermal compressibility KHS = (1-eta)^4 / (1+4eta+4eta^2-4eta^3+eta^4).
      Parameters:
      eta - packing fraction
      Returns:
      KHS

    • calcMieAlpha

      public static double calcMieAlpha(double lambdaR, double lambdaA)
      Alpha parameter for the Mie potential: alpha = C * [1/(lam_a-3) - 1/(lam_r-3)].
      Parameters:
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      Returns:
      alpha

    • calcPadeF

      public static double calcPadeF(int funcIndex, double alpha)
      Padé approximant function f_m(alpha) from Lafitte 2013 Table 3.
      Parameters:
      funcIndex - function index: 0=f1, 1=f2, 2=f3, 3=f4, 4=f5, 5=f6
      alpha - Mie alpha parameter
      Returns:
      f_m(alpha)

    • calcChi

      static double calcChi(double zetaSt, double lambdaR, double lambdaA)
      Chi correction for second-order perturbation (Lafitte 2013 Eq. 42). chi = f1(alpha) * zetaSt + f2(alpha) * zetaSt^5 + f3(alpha) * zetaSt^8, where zetaSt is the sigma-based packing fraction.
      Parameters:
      zetaSt - sigma-based packing fraction (pi/6 * rhoS * sigma^3)
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      Returns:
      chi

    • calcA2MieAtEta

      public static double calcA2MieAtEta(double eta, double zetaSt, double lambdaR, double lambdaA, double epsOverKT, double cMie, double x0)
      Full second-order Mie perturbation a2/(NkT) at a given eta value.
      Parameters:
      eta - packing fraction (d-based)
      zetaSt - sigma-based packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      cMie - Mie prefactor
      x0 - sigma/d ratio
      Returns:
      a2_Mie (dimensionless)

    • calcA3Mie

      public static double calcA3Mie(double zetaSt, double lambdaR, double lambdaA, double epsOverKT)
      Third-order Mie perturbation a3/(NkT) (Lafitte 2013 Eq. 50). a3 = -(eps/kT)^3 * f4(alpha) * zetaSt * exp(f5(alpha) * zetaSt + f6(alpha) * zetaSt^2).
      Parameters:
      zetaSt - sigma-based packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      Returns:
      a3 (dimensionless)

    • calcAS1Bare

      public static double calcAS1Bare(double eta, double lambda)
      Bare first-order Sutherland mean-field (without 12*eps*eta/kT prefactor). Following Clapeyron convention: aS1_bare = -1/(lambda-3) * gCS(zetaEff).
      Parameters:
      eta - packing fraction
      lambda - Mie exponent
      Returns:
      aS1 bare (dimensionless)

    • calcBBare

      public static double calcBBare(double eta, double lambda, double x0)
      Bare B correction (without 12*eps*eta/kT prefactor). Following Clapeyron convention.
      Parameters:
      eta - packing fraction
      lambda - Mie exponent
      x0 - sigma/d ratio
      Returns:
      B bare (dimensionless)

    • calcGHS_x0

      public static double calcGHS_x0(double eta, double x0)
      Hard-sphere RDF at separation x0 = sigma/d using the parametric form from Lafitte 2013. At x0=1 this approximates the CS contact value. At x0 > 1 it gives the correct RDF at sigma.
      Parameters:
      eta - packing fraction
      x0 - sigma/d ratio
      Returns:
      g_d^HS(sigma)

    • calcG1Chain

      public static double calcG1Chain(double eta, double lambdaR, double lambdaA, double cMie, double x0)
      First-order perturbation correction g1 for the chain RDF at sigma (Lafitte 2013 Eq. 37). Uses the bare (Clapeyron-convention) perturbation terms.
      Parameters:
      eta - packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      cMie - Mie prefactor
      x0 - sigma/d ratio
      Returns:
      g1 (dimensionless)

    • calcG2Chain

      public static double calcG2Chain(double eta, double zetaSt, double lambdaR, double lambdaA, double epsOverKT, double cMie, double x0)
      Second-order perturbation correction g2 for the chain RDF at sigma (Lafitte 2013 Eq. 38). Uses the bare perturbation terms and gamma_c correction.
      Parameters:
      eta - packing fraction
      zetaSt - sigma-based packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      cMie - Mie prefactor
      x0 - sigma/d ratio
      Returns:
      g2 (dimensionless)

    • calcGMie

      public static double calcGMie(double eta, double zetaSt, double lambdaR, double lambdaA, double epsOverKT, double cMie, double x0)
      Full Mie-corrected RDF at distance sigma for the chain contribution (Lafitte 2013 Eq. 35). g_Mie(sigma) = g_HS(x0; eta) * exp(tau * g1/g0 + tau^2 * g2/g0).
      Parameters:
      eta - packing fraction
      zetaSt - sigma-based packing fraction
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      cMie - Mie prefactor
      x0 - sigma/d ratio
      Returns:
      g_Mie(sigma)

    • calcGMieBlended

      static double calcGMieBlended(double eta, double x0, double lambdaR, double lambdaA, double epsOverKT, double cMie, double alpha)
      Blended g for the chain contribution: g_HS(x0) * exp(alpha * correction). When alpha=0, returns g_HS(x0). When alpha=1, returns full g_Mie.
      Parameters:
      eta - packing fraction
      x0 - sigma/d ratio
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      epsOverKT - epsilon/(kT)
      cMie - Mie prefactor
      alpha - blending fraction (0 to 1)
      Returns:
      blended g chain

    • computeDispersionTerms

      private void computeDispersionTerms()
      Computes first, second, and third-order perturbation terms and all their eta and T derivatives using pair summation per Lafitte 2013 Eqs. 37-40. For multi-component mixtures, a_k = sum_ij xs_i * xs_j * a_k^ij where a_k^ij uses pair-specific cross parameters (sigma_ij, eps_ij, lambda_ij). For single-component systems, falls back to direct evaluation.

    • computeDispTermsSingleFluid

      private void computeDispTermsSingleFluid(double eta, double epsOverKT, double cMie, double x0, double lambdaR, double lambdaA, double sigma, double epsk)
      Computes dispersion terms for a single effective fluid (pure component or one-fluid mapped).
      Parameters:
      eta - packing fraction
      epsOverKT - reduced energy parameter
      cMie - Mie prefactor
      x0 - sigma/d ratio
      lambdaR - repulsive exponent
      lambdaA - attractive exponent
      sigma - segment diameter in m
      epsk - energy parameter eps/k in K

    • computeDispTermsPairSum

      private void computeDispTermsPairSum(double eta)
      Computes pair-summed dispersion terms for multi-component mixtures. a_k = sum_ij xs_i * xs_j * a_k^ij per Lafitte 2013 Eqs. 37-40. Cross parameters follow Eq. 36: sigma_ij = arithmetic mean, eps_ij = geometric mean with sigma^3 correction, lambda_ij = 3 + geometric mean of (lambda-3), d_ij = BH diameter of cross potential.
      Parameters:
      eta - packing fraction

    • calcPairDispSum

      private double[] calcPairDispSum(double eta, double temp)
      Evaluates pair-summed a1, a2, a3 at given eta and temperature. For each pair (i,j), computes cross parameters per Lafitte 2013 Eq. 36 including the sigma^3-corrected epsilon and BH diameter from the cross-potential.
      Parameters:
      eta - packing fraction
      temp - temperature in K
      Returns:
      double[3] = {a1_sum, a2_sum, a3_sum}

    • calcPairDispSumPerComp

      private double[] calcPairDispSumPerComp(double eta, double temp)
      Computes per-component weighted dispersion sums: aDispPerComp[i] = sum_l xs_l * (a1_il+a2_il+a3_il). Used for the analytical dF_DISP/dNi formula.
      Parameters:
      eta - packing fraction
      temp - temperature in K
      Returns:
      array of per-component dispersion sums

    • initAssocGDerivatives

      private void initAssocGDerivatives()
      Initializes association-related quantities for the current volume. Computes the Dufal 2015 association integral I(Tr, rhoStar) and its volume derivatives using the polynomial from Dufal et al. (2015) Mol. Phys. 113(9-10), 948-984.

      The I integral replaces the simple g_HS contact RDF used in PC-SAFT. It captures the orientationally-averaged Mie potential Boltzmann factor integral and provides a more accurate association strength for SAFT-VR Mie. For the Michelsen-Hendriks dF/dV shortcut, gcpavAssoc is set to d(ln I)/dV.

    • calcAssocGMieEffective

      private double calcAssocGMieEffective(double etaVal)
      Computes the effective association g_Mie as a composition-weighted average over associating component pairs. For pure components, returns g_Mie(sigma; eta, T) directly. Uses the full Mie RDF from Lafitte 2013 Eq. 35 instead of the simple hard-sphere g_HS contact value.
      Parameters:
      etaVal - packing fraction
      Returns:
      effective g_Mie for association delta computation

    • calcD2LngDeta2

      private double calcD2LngDeta2(double eta)
      Second derivative of ln(g_HS) with respect to packing fraction eta.
      Parameters:
      eta - packing fraction
      Returns:
      d2(ln g)/d(eta)2

    • computeDeltaAssoc

      private void computeDeltaAssoc()
      Computes the association strength delta for all site-site pairs using Dufal 2015 I polynomial. Association cross parameters: epsilon_HB_ij = (eps_i + eps_j)/2 (CR-1), K_HB_ij = sqrt(K_HB_i * K_HB_j) (geometric mean).

      Uses delta = (exp(eps_HB/RT) - 1) * NA * K_HB * I(Tr, rho*) where I is the orientationally averaged association kernel from Dufal et al. (2015) Mol. Phys. 113(9-10), 948-984. The pair-specific reduced temperature Tr_ij = T / eps_disp_ij uses the dispersion cross epsilon.

    • solveAssociation

      private boolean solveAssociation()
      Solves the association fraction X_A for all sites using successive substitution. X_A = 1 / (1 + (1/V) * sum_j nj * sum_B delta_AB * X_B).
      Returns:
      true if converged

    • calcAssocHCPA

      private void calcAssocHCPA()
      Calculates hcpatot = sum_i ni * sum_A (1 - X_Ai). Also computes T derivatives.

    • calcHcpatotDerivatives

      private void calcHcpatotDerivatives()
      Calculates temperature derivatives of hcpatot using perturbation of delta. dh/dT = sum_i ni * sum_A (-dX_A/dT), where the XA T-derivative comes from implicit differentiation of the association equation.

    • solveAssociationFull

      private void solveAssociationFull()
      Performs the full association calculation sequence: compute g and its derivatives, compute delta, solve for XA, and compute hcpatot. Called inside the molarVolume solver at each iteration.

    • F_ASSOC_SAFT

      public double F_ASSOC_SAFT()
      Association Helmholtz free energy following PCSAFTa/CPA convention. F_ASSOC = sum_i ni * sum_A [ln(X_Ai) - X_Ai/2 + 1/2].
      Returns:
      F_ASSOC

    • dF_ASSOC_SAFTdV

      public double dF_ASSOC_SAFTdV()
      dF_ASSOC/dV (w.r.t. V_m3). Returns cached value computed during volInit.
      Returns:
      dF_ASSOC/dV_m3

    • dF_ASSOC_SAFTdVdV

      public double dF_ASSOC_SAFTdVdV()
      d2F_ASSOC/dV2 (w.r.t. V_m3). Computed via cloned-phase numerical differentiation to correctly capture the full implicit XA response (Michelsen-Hendriks Q-function correction). The analytical product-rule of dFdV fails for strong association (e.g., water) because hcpatot is nearly V-independent while F_ASSOC (through ln XA) varies steeply.
      Returns:
      second volume derivative in SI units (m^-6)

    • computeAssocDVDV

      private void computeAssocDVDV()
      Computes d2F_ASSOC/dV_SI2 numerically via cloned phases. Clones the phase, perturbs molar volume, runs volInit on the clone (which re-solves XA), and takes central difference of F_ASSOC. This avoids corrupting the original phase state.

    • dF_ASSOC_SAFTdT

      public double dF_ASSOC_SAFTdT()
      dF_ASSOC/dT.
      Returns:
      temperature derivative

    • dF_ASSOC_SAFTdTdT

      public double dF_ASSOC_SAFTdTdT()
      d2F_ASSOC/dT2.
      Returns:
      second temperature derivative

    • dF_ASSOC_SAFTdTdV

      public double dF_ASSOC_SAFTdTdV()
      d2F_ASSOC/dTdV (w.r.t. V_m3). Numerical via central difference.
      Returns:
      cross derivative

    • getTotalNumberOfAssociationSites

      public int getTotalNumberOfAssociationSites()
      Returns total number of association sites in this phase.
      Returns:
      total sites

    • getUseASSOC

      public int getUseASSOC()
      Returns association toggle.
      Returns:
      useASSOC

    • getHcpatot

      public double getHcpatot()
      Returns hcpatot (sum of (1-XA)*n over all sites).
      Returns:
      hcpatot

    • getGcpaAssoc

      public double getGcpaAssoc()
      Returns association g (hard-sphere contact RDF used in delta).
      Returns:
      gcpaAssoc

    • getCrossAssociationScheme

      public int getCrossAssociationScheme(int comp1, int comp2, int site1, int site2)
      Returns the cross-association scheme indicator between sites of two components.
      Parameters:
      comp1 - component 1 index
      comp2 - component 2 index
      site1 - site index on comp1
      site2 - site index on comp2
      Returns:
      1 if sites can bond, 0 otherwise

    • getDeltaAssoc

      public double getDeltaAssoc(int globalSiteA, int globalSiteB)
      Returns the global delta between two sites (global indexing).
      Parameters:
      globalSiteA - global site index A
      globalSiteB - global site index B
      Returns:
      delta value

    • getSiteOffset

      public int getSiteOffset(int compIndex)
      Returns the site offset for component i (first global site index of component i).
      Parameters:
      compIndex - component index
      Returns:
      site offset

    • getF

      public double getF()

      getF.

      Overrides:
      getF in class PhaseEos
      Returns:
      a double

    • dFdV

      public double dFdV()

      Calculate derivative of F per Volume, i.e., dF/dV.

      Specified by:
      dFdV in interface PhaseInterface
      Overrides:
      dFdV in class PhaseEos
      Returns:
      a double

    • dFdVdV

      public double dFdVdV()

      dFdVdV.

      Specified by:
      dFdVdV in interface PhaseInterface
      Overrides:
      dFdVdV in class PhaseEos
      Returns:
      a double

    • dFdVdVdV

      public double dFdVdVdV()

      dFdVdVdV.

      Overrides:
      dFdVdVdV in class PhaseEos
      Returns:
      a double

    • dFdT

      public double dFdT()

      Calculate derivative of F per Temperature, i.e., dF/dT.

      Specified by:
      dFdT in interface PhaseInterface
      Overrides:
      dFdT in class PhaseEos
      Returns:
      a double

    • dFdTdT

      public double dFdTdT()

      dFdTdT.

      Specified by:
      dFdTdT in interface PhaseInterface
      Overrides:
      dFdTdT in class PhaseEos
      Returns:
      a double

    • dFdTdV

      public double dFdTdV()

      Calculate derivative of F per Temperature and Volume, i.e., dF/dT * 1/dV.

      Specified by:
      dFdTdV in interface PhaseInterface
      Overrides:
      dFdTdV in class PhaseEos
      Returns:
      a double

    • F_HC_SAFT

      public double F_HC_SAFT()
      Hard-chain Helmholtz free energy. F_HC = n * [m * a_HS - (m-1) * ln(g_HS)]
      Returns:
      F_HC

    • dF_HC_SAFTdV

      public double dF_HC_SAFTdV()
      dF_HC/dV (w.r.t. V_m3).
      Returns:
      derivative

    • dF_HC_SAFTdVdV

      public double dF_HC_SAFTdVdV()
      d2F_HC/dV2 (w.r.t. V_m3).
      Returns:
      derivative

    • dF_HC_SAFTdT

      public double dF_HC_SAFTdT()
      dF_HC/dT.
      Returns:
      derivative

    • dF_HC_SAFTdTdT

      public double dF_HC_SAFTdTdT()
      d2F_HC/dT2.
      Returns:
      derivative

    • dF_HC_SAFTdTdV

      public double dF_HC_SAFTdTdV()
      d2F_HC/dTdV.
      Returns:
      derivative

    • F_DISP_SAFT

      public double F_DISP_SAFT()
      Total dispersion Helmholtz free energy. F_disp = n * m_bar * (a1 + a2 + a3) where m_bar is the mean segment number.
      Returns:
      F_DISP

    • dF_DISP_SAFTdV

      public double dF_DISP_SAFTdV()
      dF_DISP/dV (w.r.t. V_m3).
      Returns:
      derivative

    • dF_DISP_SAFTdVdV

      public double dF_DISP_SAFTdVdV()
      d2F_DISP/dV2 (w.r.t. V_m3).
      Returns:
      derivative

    • dF_DISP_SAFTdT

      public double dF_DISP_SAFTdT()
      dF_DISP/dT.
      Returns:
      derivative

    • dF_DISP_SAFTdTdT

      public double dF_DISP_SAFTdTdT()
      d2F_DISP/dT2.
      Returns:
      derivative

    • dF_DISP_SAFTdTdV

      public double dF_DISP_SAFTdTdV()
      d2F_DISP/dTdV.
      Returns:
      derivative

    • molarVolume

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

      molarVolume.

      Uses pressure-residual Newton-Raphson solver following PhasePCSAFTRahmat pattern.

      Specified by:
      molarVolume in interface PhaseInterface
      Overrides:
      molarVolume in class PhaseEos
      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.

    • calcmSAFT

      public double calcmSAFT()
      Calculates mean segment number m-bar = sum(xi * mi).
      Returns:
      mean segment number

    • calcmmin1SAFT

      public double calcmmin1SAFT()
      Calculates sum(xi * (mi - 1)).
      Returns:
      m-1 average

    • calcdSAFT

      public double calcdSAFT()
      Calculates sum(xi * mi * di^3).
      Returns:
      d^3 weighted sum

    • calcdmeanSAFT

      public double calcdmeanSAFT()
      Calculates mean diameter from (sum ni*mi*di^3 / sum ni*mi)^(1/3).
      Returns:
      mean diameter in m

    • getdDSAFTdT

      public double getdDSAFTdT()
      Temperature derivative of dSAFT = sum(xi * mi * 3*di^2 * ddi/dT).
      Returns:
      d(dSAFT)/dT

    • getd2DSAFTdTdT

      public double getd2DSAFTdTdT()
      Second temperature derivative of dSAFT (zero approx).
      Returns:
      d2(dSAFT)/dT2

    • getDSAFT

      public double getDSAFT()
      Returns the d^3 weighted sum for component access.
      Returns:
      dSAFT

    • getDaHSSAFTdN

      public double getDaHSSAFTdN()
      Returns daHS/deta.
      Returns:
      daHSSAFTdN

    • getDgHSSAFTdN

      public double getDgHSSAFTdN()
      Returns dgHS/deta.
      Returns:
      dgHSSAFTdN

    • getVolumeSAFT

      public double getVolumeSAFT()
      Returns the SAFT volume in m3.
      Returns:
      volumeSAFT

    • getNSAFT

      public double getNSAFT()
      Returns the packing fraction.
      Returns:
      nSAFT (eta)

    • getmSAFT

      public double getmSAFT()
      Returns mean segment number.
      Returns:
      m-bar

    • getADispPerComp

      public double getADispPerComp(int compIndex)
      Returns the per-component weighted dispersion sum: aDispPerComp[i] = sum_l xs_l * (a1_il+a2_il+a3_il). Used for analytical dF_DISP/dNi. Returns 0 if not computed (single component).
      Parameters:
      compIndex - component index
      Returns:
      per-component dispersion sum

    • getGhsSAFT

      public double getGhsSAFT()
      Returns hard-sphere RDF at contact.
      Returns:
      ghs

    • getAHSSAFT

      public double getAHSSAFT()
      Returns hard-sphere Helmholtz energy per mole.
      Returns:
      aHS

    • getMmin1SAFT

      public double getMmin1SAFT()
      Returns sum(xi*(mi-1)).
      Returns:
      mmin1SAFT

    • getUseHS

      public int getUseHS()
      Returns useHS toggle.
      Returns:
      useHS

    • getUseDISP

      public int getUseDISP()
      Returns useDISP toggle.
      Returns:
      useDISP

    • getA1Disp

      public double getA1Disp()
      Returns first-order dispersion perturbation a1/(NkT).
      Returns:
      a1Disp

    • getA2Disp

      public double getA2Disp()
      Returns second-order dispersion perturbation a2/(NkT).
      Returns:
      a2Disp

    • getA3Disp

      public double getA3Disp()
      Returns third-order dispersion perturbation a3/(NkT).
      Returns:
      a3Disp

    • getDa1DispDeta

      public double getDa1DispDeta()
      Returns da1Disp/deta derivative.
      Returns:
      da1Disp/deta

    • getDa2DispDeta

      public double getDa2DispDeta()
      Returns da2Disp/deta derivative.
      Returns:
      da2Disp/deta

    • getDa3DispDeta

      public double getDa3DispDeta()
      Returns da3Disp/deta derivative.
      Returns:
      da3Disp/deta

    • getdDSAFTdTprime

      public double getdDSAFTdTprime(double mi, double di, double nMoles)
      Returns the d(dSAFT)/dN_i helper for composition differentiation.
      Parameters:
      mi - segment number of component i
      di - diameter of component i
      nMoles - total moles
      Returns:
      d(dSAFT)/dN_i contribution (partial through d^3 sum)

    • getF1dispVolTerm

      public double getF1dispVolTerm()
      Returns the dispersion volume term (legacy, stub).
      Returns:
      0.0

    • getF1dispSumTerm

      public double getF1dispSumTerm()
      Returns F1 dispersion sum term (legacy, stub).
      Returns:
      0.0

    • getF2dispSumTerm

      public double getF2dispSumTerm()
      Returns F2 dispersion sum term (legacy, stub).
      Returns:
      0.0

    • calcF1dispI1dN

      public double calcF1dispI1dN()
      Returns I1 integral derivative w.r.t. eta (legacy, stub).
      Returns:
      0.0

    • calcF1dispI1dm

      public double calcF1dispI1dm()
      Returns I1 integral derivative w.r.t. m (legacy, stub).
      Returns:
      0.0

    • calcF2dispI2dN

      public double calcF2dispI2dN()
      Returns I2 integral derivative w.r.t. eta (legacy, stub).
      Returns:
      0.0

    • calcF2dispI2dm

      public double calcF2dispI2dm()
      Returns I2 integral derivative w.r.t. m (legacy, stub).
      Returns:
      0.0

    • getF2dispZHC

      public double getF2dispZHC()
      Returns the ZHC factor (legacy, stub returns 1.0 to avoid division by zero).
      Returns:
      1.0