Class PhaseElectrolyteCPAMM

All Implemented Interfaces:
Serializable, Cloneable, PhaseCPAInterface, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface

public class PhaseElectrolyteCPAMM extends PhaseSrkCPA
Electrolyte CPA (e-CPA) phase class implementing the Maribo-Mogensen model.

The residual Helmholtz free energy is computed as:

A^res = A^SRK + A^Association + A^Debye-Hückel + A^Born

Key features of this model:

  • SRK cubic equation of state with CPA association term
  • Debye-Hückel theory for long-range electrostatic interactions (simpler than MSA)
  • Born solvation term with empirical Born radius correlations
  • Temperature-dependent ion-solvent interaction parameters

References:

  • Maribo-Mogensen, B. (2014). PhD Thesis, DTU Chemical Engineering.
  • Maribo-Mogensen et al., Ind. Eng. Chem. Res. 2012, 51, 5353-5363
Version:
$Id: $Id
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.
    • VACUUM_PERMITTIVITY

      private static final double VACUUM_PERMITTIVITY
      Vacuum permittivity [F/m] = [C²/(J·m)].
      See Also:
    • ELECTRON_CHARGE

      private static final double ELECTRON_CHARGE
      Elementary charge [C].
      See Also:
    • AVOGADRO

      private static final double AVOGADRO
      Avogadro's number [1/mol].
      See Also:
    • BOLTZMANN

      private static final double BOLTZMANN
      Boltzmann constant [J/K].
      See Also:
    • solventPermittivity

      protected double solventPermittivity
      Solvent static permittivity (dielectric constant) [-].
    • solventPermittivitydT

      protected double solventPermittivitydT
      Temperature derivative of solvent permittivity [1/K].
    • solventPermittivitydTdT

      protected double solventPermittivitydTdT
      Second temperature derivative of solvent permittivity [1/K²].
    • mixturePermittivity

      protected double mixturePermittivity
      Mixture dielectric constant including ion contributions [-].
    • kappa

      protected double kappa
      Debye screening length inverse (kappa) [1/m].
    • kappadT

      protected double kappadT
      Temperature derivative of kappa [1/(m·K)].
    • bornX

      protected double bornX
      Born X factor: Σ(n_i * z_i² / R_Born,i) [mol/m].
    • ionSolventW

      protected double ionSolventW
      Ion-solvent short-range interaction parameter W [J/mol]. Calculated from wij mixing rule parameters which are populated from ΔU_iw values.
    • ionSolventWdT

      protected double ionSolventWdT
      Temperature derivative of W [J/(mol·K)].
    • packingFraction

      protected double packingFraction
      Packing fraction η = (π N_A / 6V) Σ n_i σ_i³ [-].
    • packingFractiondV

      protected double packingFractiondV
      Volume derivative of packing fraction [1/m³].
    • meanIonDiameter

      protected double meanIonDiameter
      Mean closest-approach diameter for extended Debye-Hückel [m]. Computed as the charge-weighted average of ion LJ diameters: d = Σ(n_i z_i² σ_i) / Σ(n_i z_i²).
    • tauDH

      protected double tauDH
      DH screening factor τ(χ) = 3ψ(χ)/χ³ where χ = κd [-]. Reduces the primitive DH (DHLL) contribution to account for finite ion size. τ = 1.0 recovers the DHLL; τ ≈ 0.6 at 1 molal NaCl.
    • tauDHprime

      protected double tauDHprime
      Derivative of screening factor dτ/dχ [-]. Used for composition derivatives where κ changes significantly with added ions.
    • wij

      protected double[][] wij
      Ion-solvent binary interaction parameters wij. These are populated from the ΔU_iw parameters from Maribo-Mogensen thesis Table 6.11. Format: wij[i][j] where i is ion index and j is solvent index. Units: Kelvin (energy/R).
    • wijT

      protected double[][] wijT
      Temperature derivatives of ion-solvent wij parameters. wijT[i][j] = d(wij)/dT. Units: K/K (dimensionless).
    • debyeHuckelOn

      protected boolean debyeHuckelOn
      Flag to enable/disable Debye-Hückel term.
    • bornOn

      protected boolean bornOn
      Flag to enable/disable Born term.
    • shortRangeOn

      protected boolean shortRangeOn
      Flag to enable/disable short-range ion-solvent term. Uses the Furst electrolyte mixing rule correlation for wij parameters, combined with the FSR2 Helmholtz framework. Enable this via setShortRangeOn(true).
    • electrolyteMixingRule

      protected transient ElectrolyteMixingRulesInterface electrolyteMixingRule
      Electrolyte mixing rule for computing W, Wi, WiT (Furst-style wij from Stokes diameter).
    • srW

      protected double srW
      Short-range W parameter from electrolyte mixing rule [m³·mol].
    • srWT

      protected double srWT
      Temperature derivative of W [m³·mol/K].
    • srWTT

      protected double srWTT
      Second temperature derivative of W [m³·mol/K²].
    • dielectricMixingRule

      protected PhaseElectrolyteCPAMM.DielectricMixingRule dielectricMixingRule
      Current dielectric mixing rule.
  • Constructor Details

    • PhaseElectrolyteCPAMM

      public PhaseElectrolyteCPAMM()
      Constructor for PhaseElectrolyteCPAMM.
  • Method Details

    • 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 PhaseSrkCPA
      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.
    • clone

      public PhaseElectrolyteCPAMM clone()

      clone.

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

      public void initElectrolyteProperties()
      Initialize volume-dependent electrolyte properties. Called during molar volume iteration.
    • calcMeanIonDiameter

      public double calcMeanIonDiameter()
      Calculate the mean closest-approach diameter for the extended Debye-Hückel. Uses the charge-squared weighted average of ionic LJ diameters.
      Returns:
      mean ion diameter [m]
    • psiDH

      static double psiDH(double u)
      Extended DH charging function ψ(u) = ln(1+u) - u + u²/2. In the limit u → 0, ψ → u³/3 (recovers primitive DH).
      Parameters:
      u - dimensionless argument κd
      Returns:
      ψ(u) [-]
    • psiDHprime

      static double psiDHprime(double u)
      First derivative of ψ: ψ'(u) = u²/(1+u). In the limit u → 0, ψ' → u².
      Parameters:
      u - dimensionless argument κd
      Returns:
      ψ'(u) [-]
    • calcTauDH

      static double calcTauDH(double chi)
      DH screening factor τ(χ) = 3ψ(χ)/χ³. Equals 1.0 in the primitive DH limit (χ → 0), with values less than 1 at finite concentration. For NaCl at 1 molal, τ ≈ 0.60.
      Parameters:
      chi - dimensionless screening parameter κd
      Returns:
      τ [-]
    • calcTauDHprime

      static double calcTauDHprime(double chi)
      Derivative of DH screening factor: dτ/dχ. Used for the composition derivative correction of the extended DH where κ changes with ion amount.
      Parameters:
      chi - dimensionless screening parameter κd
      Returns:
      dτ/dχ [-]
    • getTauDH

      public double getTauDH()
      Get the DH screening factor τ(κd) [-].
      Returns:
      τ value (1.0 for primitive DH, less than 1 for extended DH)
    • getTauDHprime

      public double getTauDHprime()
      Get the derivative of DH screening factor dτ/dχ [-].
      Returns:
      dτ/d(κd)
    • getMeanIonDiameter

      public double getMeanIonDiameter()
      Get the mean ion closest-approach diameter [m].
      Returns:
      mean ion diameter in meters
    • initMixingRuleWij

      public void initMixingRuleWij()
      Initialize the ion-solvent wij mixing rule parameters from component ΔU_iw values.

      Following the Maribo-Mogensen thesis, the ion-solvent short-range interaction is parameterized as: ΔU_iw(T) = u⁰_iw + uᵀ_iw × (T - 298.15)

      These are stored in the wij/wijT arrays which follow the standard mixing rule convention. The W parameter used in the short-range Helmholtz term is then calculated as: W = Σ_i Σ_j (n_i × n_j × wij(T)) for ion i and solvent j.

    • getWij

      public double getWij(int i, int j)
      Get the ion-solvent interaction parameter wij at current temperature.
      Parameters:
      i - first component index
      j - second component index
      Returns:
      wij value [K]
    • getWijT

      public double getWijT(int i, int j)
      Get the temperature derivative of ion-solvent interaction parameter wij.
      Parameters:
      i - first component index
      j - second component index
      Returns:
      dwij/dT [K/K]
    • setDielectricMixingRule

      public void setDielectricMixingRule(PhaseElectrolyteCPAMM.DielectricMixingRule rule)
      Set the dielectric constant mixing rule.
      Parameters:
      rule - the mixing rule to use
    • getDielectricMixingRule

      public PhaseElectrolyteCPAMM.DielectricMixingRule getDielectricMixingRule()
      Get the current dielectric constant mixing rule.
      Returns:
      the current mixing rule
    • calcSolventPermittivity

      public double calcSolventPermittivity(double T)
      Calculate solvent static permittivity using the selected mixing rule.

      This method dispatches to the appropriate mixing rule implementation based on the dielectricMixingRule setting.

      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivityMolarAvg

      private double calcSolventPermittivityMolarAvg(double T)
      Calculate solvent permittivity using molar average mixing rule.
      ε_mix = Σ(x_i × ε_i) / Σ(x_i)
      
      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivityVolumeAvg

      private double calcSolventPermittivityVolumeAvg(double T)
      Calculate solvent permittivity using volume average mixing rule.
      ε_mix = Σ(φ_i × ε_i)
      where φ_i = (n_i × V_i) / Σ(n_j × V_j) is the volume fraction
      

      This gives better accuracy for water-glycol mixtures (2.4% avg error vs 4.2% for molar average) based on Ma et al. (2010) data.

      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivityLooyenga

      private double calcSolventPermittivityLooyenga(double T)
      Calculate solvent permittivity using Looyenga mixing rule.
      ε_mix^(1/3) = Σ(φ_i × ε_i^(1/3))
      

      The Looyenga equation has theoretical basis for polar molecule mixtures and works well for water-organic systems.

      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivityOster

      private double calcSolventPermittivityOster(double T)
      Calculate solvent permittivity using Oster mixing rule.

      The Oster equation is specifically designed for water-alcohol mixtures. It accounts for the non-ideal mixing behavior of polar solvents.

      ε_mix = ε₁×φ₁ + ε₂×φ₂ + k×φ₁×φ₂×(ε₁ - ε₂)
      

      where k is an interaction parameter (typically -0.2 to -0.5 for water-alcohol).

      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivityLichtenecker

      private double calcSolventPermittivityLichtenecker(double T)
      Calculate solvent permittivity using Lichtenecker-Rother mixing rule.
      ln(ε_mix) = Σ(φ_i × ln(ε_i))
      

      This is equivalent to: ε_mix = ε₁^φ₁ × ε₂^φ₂ × ...

      Parameters:
      T - temperature [K]
      Returns:
      solvent permittivity [-]
    • calcSolventPermittivitydT

      public double calcSolventPermittivitydT(double T)
      Calculate temperature derivative of solvent permittivity using the selected mixing rule.
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTMolarAvg

      private double calcSolventPermittivitydTMolarAvg(double T)
      Molar average temperature derivative.
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTVolumeAvg

      private double calcSolventPermittivitydTVolumeAvg(double T)
      Volume average temperature derivative.
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTLooyenga

      private double calcSolventPermittivitydTLooyenga(double T)
      Looyenga temperature derivative.
      d(ε^(1/3))/dT = Σ(φ_i × (1/3) × ε_i^(-2/3) × dε_i/dT)
      dε/dT = 3 × ε^(2/3) × d(ε^(1/3))/dT
      
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTLichtenecker

      private double calcSolventPermittivitydTLichtenecker(double T)
      Lichtenecker temperature derivative.
      d(ln ε)/dT = Σ(φ_i × (1/ε_i) × dε_i/dT)
      dε/dT = ε × d(ln ε)/dT
      
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTOster

      private double calcSolventPermittivitydTOster(double T)
      Oster temperature derivative.
      Parameters:
      T - temperature [K]
      Returns:
      dε/dT [1/K]
    • calcSolventPermittivitydTdT

      public double calcSolventPermittivitydTdT(double T)
      Calculate second temperature derivative of solvent permittivity.

      Currently only implemented for molar average. Other mixing rules use an approximate second derivative based on the molar average.

      Parameters:
      T - temperature [K]
      Returns:
      d²ε/dT² [1/K²]
    • calcSolventPermittivitydn

      public double calcSolventPermittivitydn(int k, double T)
      Calculate composition derivative of solvent permittivity for component k.

      For molar average: dε/dn_k = (ε_k - ε_mix) / n_solvent

      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      dε/dn_k [-]
    • calcSolventPermittivitydnMolarAvg

      private double calcSolventPermittivitydnMolarAvg(int k, double T)
      Molar average composition derivative.
      ε = Σ(n_i × ε_i) / Σ(n_i)
      dε/dn_k = (ε_k - ε) / n_solvent
      
      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      dε/dn_k [-]
    • calcSolventPermittivitydnVolumeAvg

      private double calcSolventPermittivitydnVolumeAvg(int k, double T)
      Volume average composition derivative.
      ε = Σ(φ_i × ε_i) where φ_i = (n_i × V_i) / Σ(n_j × V_j)
      dε/dn_k = V_k × (ε_k - ε) / Σ(n_j × V_j)
      
      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      dε/dn_k [-]
    • calcSolventPermittivitydnLooyenga

      private double calcSolventPermittivitydnLooyenga(int k, double T)
      Looyenga composition derivative.
      ε^(1/3) = Σ(φ_i × ε_i^(1/3))
      d(ε^(1/3))/dn_k = V_k × (ε_k^(1/3) - ε^(1/3)) / Σ(n_j × V_j)
      dε/dn_k = 3 × ε^(2/3) × d(ε^(1/3))/dn_k
      
      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      dε/dn_k [-]
    • calcSolventPermittivitydnLichtenecker

      private double calcSolventPermittivitydnLichtenecker(int k, double T)
      Lichtenecker composition derivative.
      ln(ε) = Σ(φ_i × ln(ε_i))
      d(ln ε)/dn_k = V_k × (ln(ε_k) - ln(ε)) / Σ(n_j × V_j)
      dε/dn_k = ε × d(ln ε)/dn_k
      
      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      dε/dn_k [-]
    • calcSolventPermittivitydndT

      public double calcSolventPermittivitydndT(int k, double T)
      Calculate mixed composition-temperature derivative of solvent permittivity.
      Parameters:
      k - component index
      T - temperature [K]
      Returns:
      d²ε/(dn_k dT) [-/K]
    • calcSolventPermittivitydndn

      public double calcSolventPermittivitydndn(int k, int l, double T)
      Calculate second composition derivative of solvent permittivity.
      For molar average: ε = Σ(n_i × ε_i) / Σ(n_i)
      dε/dn_k = (ε_k - ε) / n_solvent
      d²ε/(dn_k dn_l) = -dε/dn_k / n_solvent - dε/dn_l / n_solvent
                      = -(ε_k - ε)/n_solvent² - (ε_l - ε)/n_solvent²
                      = (2ε - ε_k - ε_l) / n_solvent²
      
      Parameters:
      k - first component index
      l - second component index
      T - temperature [K]
      Returns:
      d²ε/(dn_k dn_l) [-]
    • calcMixturePermittivity

      public double calcMixturePermittivity()
      Calculate mixture permittivity including ion effects. Following Maribo-Mogensen Eq. 5.27: ε = ε_solvent * (1 - η_ion) / (1 + η_ion/2) where η_ion is the ionic packing fraction.
      Returns:
      mixture permittivity [-]
    • calcIonicPackingFraction

      public double calcIonicPackingFraction()
      Calculate ionic packing fraction. η_ion = (π/6V) * Σ(n_i * σ_i³) for ionic species only
      Returns:
      ionic packing fraction [-]
    • calcKappa

      public double calcKappa()
      Calculate Debye screening parameter kappa. From Maribo-Mogensen Eq. 4.14:
      κ² = (e² / (ε₀ε_r k_B T)) * Σ(ρ_i * z_i²)
      
      where ρ_i is number density [1/m³].
      Returns:
      kappa [1/m]
    • calcKappadT

      public double calcKappadT()
      Calculate temperature derivative of kappa.
      Returns:
      dκ/dT [1/(m·K)]
    • calcBornX

      public double calcBornX()
      Calculate Born solvation factor X_Born. X_Born = Σ(n_i * z_i² / R_Born,i)

      Born radius correlations from Maribo-Mogensen Table 6.6:

      • Cations: R_Born = 0.5*σ + 0.1 Å
      • Anions: R_Born = 0.5*σ + 0.85 Å
      Returns:
      X_Born [mol/m]
    • calcBornRadius

      public double calcBornRadius(double sigma, int charge)
      Calculate Born radius using Maribo-Mogensen empirical correlations.
      Parameters:
      sigma - Lennard-Jones diameter [m]
      charge - ionic charge (sign indicates cation/anion)
      Returns:
      Born radius [m]
    • calcPackingFraction

      public double calcPackingFraction()
      Calculate packing fraction (reduced density) η.

      The packing fraction is defined as:

      η = (π N_A / 6V) × Σ n_i σ_i³
      

      This is used in the short-range term following the Furst electrolyte model.

      Returns:
      packing fraction η [-]
    • calcPackingFractiondV

      public double calcPackingFractiondV()
      Calculate volume derivative of packing fraction.

      dη/dV = -η/V

      Returns:
      dη/dV [1/m³]
    • getPackingFraction

      public double getPackingFraction()
      Get the packing fraction.
      Returns:
      packing fraction η [-]
    • getPackingFractiondV

      public double getPackingFractiondV()
      Get the volume derivative of packing fraction.
      Returns:
      dη/dV [1/m³]
    • calcIonSolventW

      public double calcIonSolventW()
      Calculate ion-solvent short-range interaction parameter W using wij mixing rule parameters.

      Following the standard mixing rule convention, W is calculated as:

      W = Σ_i Σ_j (n_i × n_j × wij(T) × R / V)
      

      where wij contains the ΔU_iw values from Maribo-Mogensen thesis Table 6.11, populated via initMixingRuleWij().

      Returns:
      W [K·mol²]
    • calcIonSolventWdT

      public double calcIonSolventWdT()
      Calculate temperature derivative of W using wijT mixing rule parameters.
      Returns:
      dW/dT [mol²] (wijT is dimensionless: d(wij)/dT in K/K = 1)
    • FDebyeHuckel

      public double FDebyeHuckel()
      Debye-Hückel contribution to Helmholtz free energy (extensive). From Maribo-Mogensen Eq. 4.19, the intensive form is:
      A^DH / (V k_B T) = -κ³ / (12π)
      
      Converting to extensive F = A/(RT):
      A^DH = -κ³ V k_B T / (12π)
      F^DH = A^DH / (RT) = -κ³ V / (12π N_A)
      

      The extended DH applies a screening factor τ(χ) = 3ψ(χ)/χ³ where χ = κd and d is the mean ion closest-approach diameter. This accounts for finite ion size and reduces the DH contribution by ~40% at 1 molal NaCl. F^DH_ext = F^DH_DHLL × τ(κd).

      Returns:
      F^DH contribution to Helmholtz energy [-]
    • dFDebyeHuckeldT

      public double dFDebyeHuckeldT()
      Temperature derivative of Debye-Hückel term at constant V. dF^DH/dT

      Uses extended DH: dF/dT = -V/(12π N_A) × dκ³/dT × τ(κd). The dτ/dT correction is neglected since it is O(0.03%/K) relative to the main term.

      Returns:
      dF^DH/dT [-/K]
    • dFDebyeHuckeldV

      public double dFDebyeHuckeldV()
      Volume derivative of Debye-Hückel term. dF^DH/dV = ∂F^DH/∂V at constant κ.

      Uses extended DH: dF/dV = -κ³/(12π N_A) × τ(κd). Since τ is independent of V at constant κ, and F is linear in V, this is exact.

      Returns:
      dF^DH/dV [1/m³]
    • FBorn

      public double FBorn()
      Born solvation contribution to Helmholtz free energy (extensive). From Maribo-Mogensen Eq. 4.23, converted to extensive form:
      F^Born = (N_A * e² / (8π * ε₀ * R * T)) * (1/ε_r - 1) * X_Born
      
      where X_Born = Σ(n_i * z_i² / R_Born,i). This is extensive (scales with n).
      Returns:
      F^Born contribution [-]
    • dFBorndT

      public double dFBorndT()
      Temperature derivative of Born term. dF^Born/dT
      Returns:
      dF^Born/dT [-/K]
    • dFBorndV

      public double dFBorndV()
      Volume derivative of Born term. dF^Born/dV (Born term is almost independent of V)
      Returns:
      dF^Born/dV [1/m³]
    • FBornD

      public double FBornD()
      Derivative of Born term with respect to dielectric constant.
      F^Born = prefactor * (1/ε - 1) * X_Born
      ∂F^Born/∂ε = prefactor * (-1/ε²) * X_Born
      
      This is used for the composition derivative via chain rule: dF^Born/dn_i includes (∂F^Born/∂ε) * (∂ε/∂n_i).
      Returns:
      ∂F^Born/∂ε [-]
    • FBornX

      public double FBornX()
      Derivative of Born term with respect to X_Born.
      F^Born = prefactor * (1/ε - 1) * X_Born
      ∂F^Born/∂X_Born = prefactor * (1/ε - 1)
      
      Returns:
      ∂F^Born/∂X_Born [-]
    • FBornDD

      public double FBornDD()
      Second derivative of Born term with respect to dielectric constant.
      ∂²F^Born/∂ε² = prefactor * 2 * X_Born / ε³
      
      Returns:
      ∂²F^Born/∂ε² [-]
    • FBornDX

      public double FBornDX()
      Mixed derivative of Born term with respect to dielectric constant and X_Born.
      ∂²F^Born/(∂ε ∂X_Born) = prefactor * (-1/ε²)
      
      Returns:
      ∂²F^Born/(∂ε ∂X_Born) [-]
    • FShortRange

      public double FShortRange()
      Short-range ion-solvent contribution to Helmholtz free energy using the Maribo-Mogensen formula.
      F_SR = W / (n_T × T × (1 - η))
      

      where W = Σ_i Σ_j n_i n_j ΔU_ij(T) is computed from the ion-solvent wij parameters (from IonParametersMM ΔU_iw values), n_T is total moles, T is temperature, and η is the packing fraction.

      Returns:
      F_SR contribution [-]
    • FSR2W

      public double FSR2W()
      Partial derivative of FSR2 with respect to W.
      Returns:
      dFSR2/dW [1/(m³)]
    • FSR2V

      public double FSR2V()
      Partial derivative of FSR2 with respect to V.
      Returns:
      dFSR2/dV [1/m³]
    • FSR2eps

      public double FSR2eps()
      Partial derivative of FSR2 with respect to eps (packing fraction).
      Returns:
      dFSR2/deps [-]
    • FSR2VV

      public double FSR2VV()
      d²FSR2/(dV²).
      Returns:
      d²FSR2/dV² [1/m⁶]
    • FSR2epsV

      public double FSR2epsV()
      d²FSR2/(dV deps).
      Returns:
      d²FSR2/(dV deps)
    • FSR2epseps

      public double FSR2epseps()
      d²FSR2/(deps²).
      Returns:
      d²FSR2/deps²
    • FSR2VW

      public double FSR2VW()
      d²FSR2/(dV dW).
      Returns:
      d²FSR2/(dV dW)
    • FSR2epsW

      public double FSR2epsW()
      d²FSR2/(deps dW).
      Returns:
      d²FSR2/(deps dW)
    • FSR2VVV

      public double FSR2VVV()
      d³FSR2/(dV³).
      Returns:
      d³FSR2/dV³
    • FSR2epsepsV

      public double FSR2epsepsV()
      d³FSR2/(deps² dV).
      Returns:
      d³FSR2/(deps² dV)
    • FSR2VVeps

      public double FSR2VVeps()
      d³FSR2/(dV² deps).
      Returns:
      d³FSR2/(dV² deps)
    • FSR2epsepseps

      public double FSR2epsepseps()
      d³FSR2/(deps³).
      Returns:
      d³FSR2/deps³
    • dFShortRangedT

      public double dFShortRangedT()
      Temperature derivative of MM short-range term: dF_SR/dT.

      F = W/(nT × T × (1-η)), so dF/dT = Wt/(nT × T × (1-η)) - W/(nT × T² × (1-η))

      Returns:
      dFSR/dT [-/K]
    • dFShortRangedTdT

      public double dFShortRangedTdT()
      Second temperature derivative of MM short-range term.
      Returns:
      d²FSR/dT² [-/K²]
    • dFShortRangedV

      public double dFShortRangedV()
      Volume derivative of MM short-range term. The MM SR is F = W/(nT × T × (1-η)). Since η depends on V, dF/dV = F × (dη/dV) / (1-η). For dilute solutions η is small and this is negligible.
      Returns:
      dFSR/dV [1/m³]
    • dFShortRangedVdV

      public double dFShortRangedVdV()
      Second volume derivative of MM short-range term. Approximately zero for dilute solutions.
      Returns:
      d²FSR/dV²
    • dFShortRangedVdVdV

      public double dFShortRangedVdVdV()
      Third volume derivative of MM short-range term.
      Returns:
      d³FSR/dV³
    • dFShortRangedTdV

      public double dFShortRangedTdV()
      Mixed temperature-volume derivative of MM short-range term.
      Returns:
      d²FSR/dTdV [1/(m³·K)]
    • calcWi

      public double calcWi(int compNumb, double temperature, double pressure, int numbcomp)
      Compute Wi = dW/dn_i for component i from the electrolyte mixing rule.
      Parameters:
      compNumb - component index
      temperature - temperature in K
      pressure - pressure in bar
      numbcomp - number of components
      Returns:
      Wi value
    • calcMMWi

      public double calcMMWi(int compNumb)
      Compute ∂W/∂n_i for the MM short-range from the wij matrix. W = Σ_a Σ_b n_a n_b wij[a][b], so ∂W/∂n_i = 2 Σ_b n_b wij[i][b].
      Parameters:
      compNumb - component index
      Returns:
      dW/dn_i [K·mol]
    • calcMMWiT

      public double calcMMWiT(int compNumb)
      Compute d²W/(dn_i dT) for the MM short-range from the wijT matrix. Since W_T = Σ_a Σ_b n_a n_b wijT[a][b], the derivative is ∂W_T/∂n_i = 2 Σ_b n_b wijT[i][b].
      Parameters:
      compNumb - component index
      Returns:
      d²W/(dn_i dT) [mol] (wijT is dimensionless)
    • calcWiT

      public double calcWiT(int compNumb, double temperature, double pressure, int numbcomp)
      Compute WiT = d²W/(dn_i dT) for component i (Furst mixing rule version).
      Parameters:
      compNumb - component index
      temperature - temperature in K
      pressure - pressure in bar
      numbcomp - number of components
      Returns:
      WiT value
    • calcWij

      public double calcWij(int i, int j, double temperature, double pressure, int numbcomp)
      Compute Wij = d²W/(dn_i dn_j) from the electrolyte mixing rule.
      Parameters:
      i - first component index
      j - second component index
      temperature - temperature in K
      pressure - pressure in bar
      numbcomp - number of components
      Returns:
      Wij value
    • getElectrolyteMixingRule

      public ElectrolyteMixingRulesInterface getElectrolyteMixingRule()
      Get the electrolyte mixing rule.
      Returns:
      the electrolyte mixing rule interface
    • getSrW

      public double getSrW()
      Get the W parameter from short-range mixing rule.
      Returns:
      W value
    • getSrWT

      public double getSrWT()
      Get the WT parameter (temperature derivative of W).
      Returns:
      WT value
    • dFDebyeHuckeldVdV

      public double dFDebyeHuckeldVdV()
      Second volume derivative of Debye-Hückel term. d²F^DH/dV²
      Returns:
      d²F^DH/dV² [1/m⁶]
    • dFDebyeHuckeldVdVdV

      public double dFDebyeHuckeldVdVdV()
      Third volume derivative of Debye-Hückel term. d³F^DH/dV³
      Returns:
      d³F^DH/dV³ [1/m⁹]
    • dFDebyeHuckeldTdT

      public double dFDebyeHuckeldTdT()
      Second temperature derivative of Debye-Hückel term. d²F^DH/dT²

      Uses extended DH screening factor τ.

      Returns:
      d²F^DH/dT² [-/K²]
    • dFDebyeHuckeldTdV

      public double dFDebyeHuckeldTdV()
      Mixed temperature-volume derivative of Debye-Hückel term. d²F^DH/dTdV

      Uses extended DH screening factor τ.

      Returns:
      d²F^DH/dTdV [1/(m³·K)]
    • dFBorndTdT

      public double dFBorndTdT()
      Second temperature derivative of Born term. d²F^Born/dT²
      Returns:
      d²F^Born/dT² [-/K²]
    • getF

      public double getF()

      getF.

      Overrides:
      getF in class PhaseSrkCPA
      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 PhaseSrkCPA
      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 PhaseSrkCPA
      Returns:
      a double
    • dFdVdV

      public double dFdVdV()

      dFdVdV.

      Specified by:
      dFdVdV in interface PhaseInterface
      Overrides:
      dFdVdV in class PhaseSrkCPA
      Returns:
      a double
    • dFdVdVdV

      public double dFdVdVdV()

      dFdVdVdV.

      Overrides:
      dFdVdVdV in class PhaseSrkCPA
      Returns:
      a double
    • dFdTdT

      public double dFdTdT()

      dFdTdT.

      Specified by:
      dFdTdT in interface PhaseInterface
      Overrides:
      dFdTdT in class PhaseSrkCPA
      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 PhaseSrkCPA
      Returns:
      a double
    • getSolventPermittivity

      public double getSolventPermittivity()
      Get solvent permittivity.
      Returns:
      solvent permittivity [-]
    • getSolventPermittivitydT

      public double getSolventPermittivitydT()
      Get temperature derivative of solvent permittivity.
      Returns:
      d(solvent permittivity)/dT [1/K]
    • getMixturePermittivity

      public double getMixturePermittivity()
      Get mixture permittivity.
      Returns:
      mixture permittivity [-]
    • getKappa

      public double getKappa()
      Get Debye screening parameter.
      Specified by:
      getKappa in interface PhaseInterface
      Overrides:
      getKappa in class PhaseEos
      Returns:
      kappa [1/m]
    • getDebyeLength

      public double getDebyeLength()
      Get Debye screening length.
      Returns:
      1/kappa [m]
    • getBornX

      public double getBornX()
      Get the Born X parameter (sum of z_i²/R_Born,i for all ions).
      Returns:
      Born X parameter [1/m]
    • getIonSolventW

      public double getIonSolventW()
      Get the ion-solvent interaction parameter W.
      Returns:
      W parameter [J/mol]
    • FDebyeHuckelX

      public double FDebyeHuckelX()
      Partial derivative of F^DH with respect to X_DH = Σ(n_i z_i²). F^DH = -κ³V/(12πRTn_T), where κ² ∝ X_DH
      Returns:
      ∂F^DH/∂X_DH
    • calcIonicStrengthSum

      public double calcIonicStrengthSum()
      Calculate ionic strength sum X_DH = Σ(n_i z_i²).
      Returns:
      ionic strength sum [mol]
    • FDebyeHuckelN

      public double FDebyeHuckelN()
      Partial derivative of F^DH with respect to n_T (at constant κ).
      Returns:
      ∂F^DH/∂n_T at constant κ
    • FBornN

      public double FBornN()
      Partial derivative of F^Born with respect to n_T (at constant X_Born).
      Returns:
      ∂F^Born/∂n_T at constant X_Born
    • FBornXT

      public double FBornXT()
      Temperature derivative of FBornX.
      Returns:
      ∂²F^Born/(∂X_Born ∂T)
    • setDebyeHuckelOn

      public void setDebyeHuckelOn(boolean on)
      Enable or disable Debye-Hückel term.
      Parameters:
      on - true to enable
    • setBornOn

      public void setBornOn(boolean on)
      Enable or disable Born term.
      Parameters:
      on - true to enable
    • setShortRangeOn

      public void setShortRangeOn(boolean on)
      Enable or disable short-range ion-solvent term.
      Parameters:
      on - true to enable
    • isShortRangeOn

      public boolean isShortRangeOn()
      Check if short-range ion-solvent term is enabled.
      Returns:
      true if short-range term is enabled
    • isDebyeHuckelOn

      public boolean isDebyeHuckelOn()
      Check if Debye-Hückel term is enabled.
      Returns:
      true if DH term is enabled
    • isBornOn

      public boolean isBornOn()
      Check if Born solvation term is enabled.
      Returns:
      true if Born term is enabled