saft updates

src/main/java/neqsim/thermo/component/ComponentPCSAFT.java

@@ -111,6 +111,55 @@ public void Finit(PhaseInterface phase, double temp, double pres, double totMole
     // System.out.println("fugacity " + getFugacityCoefficient());
   }
 
+  /**
+   * Initializes only the PC-SAFT specific component derivative quantities from the current phase
+   * state. Unlike {@link #Finit}, this method does not call super.Finit() and therefore avoids
+   * triggering dFdNdV/dFdNdN computation, which is essential for preventing infinite recursion when
+   * those methods are implemented numerically.
+   *
+   * @param phase the phase to compute quantities for
+   * @param numberOfComponents number of components in the phase
+   * @param temp temperature in K
+   * @param pres pressure in bara
+   */
+  void initSAFTDerivatives(PhaseInterface phase, int numberOfComponents, double temp, double pres) {
+    setDnSAFTdi(calcdnSAFTdi(phase, numberOfComponents, temp, pres));
+    setDghsSAFTdi(calcdghsSAFTdi(phase, numberOfComponents, temp, pres));
+    setDlogghsSAFTdi(1.0 / ((PhasePCSAFT) phase).getGhsSAFT() * getDghsSAFTdi());
+    setDmSAFTdi(calcdmSAFTdi(phase, numberOfComponents, temp, pres));
+    setdahsSAFTdi(calcdahsSAFTdi(phase, numberOfComponents, temp, pres));
+
+    F1dispVolTermdn = 1.0 * ThermodynamicConstantsInterface.avagadroNumber
+        / ((PhasePCSAFT) phase).getVolumeSAFT();
+    F2dispVolTermdn = F1dispVolTermdn;
+    F1dispSumTermdn = calcF1dispSumTermdn(phase, numberOfComponents, temp, pres);
+    F2dispSumTermdn = calcF2dispSumTermdn(phase, numberOfComponents, temp, pres);
+    F1dispI1dn = ((PhasePCSAFT) phase).calcF1dispI1dN() * getDnSAFTdi()
+        + ((PhasePCSAFT) phase).calcF1dispI1dm() * getDmSAFTdi();
+    F2dispI2dn = ((PhasePCSAFT) phase).calcF2dispI2dN() * getDnSAFTdi()
+        + ((PhasePCSAFT) phase).calcF2dispI2dm() * getDmSAFTdi();
+    F2dispZHCdn = ((PhasePCSAFT) phase).getF2dispZHCdN() * getDnSAFTdi()
+        + ((PhasePCSAFT) phase).getF2dispZHCdm() * getDmSAFTdi();
+  }
+
+  /**
+   * Reinitializes PC-SAFT quantities on a perturbed phase clone. Calls volInit() on the phase and
+   * then initSAFTDerivatives() on each component, without triggering the full Finit chain.
+   *
+   * @param perturbedPhase the cloned and perturbed phase
+   * @param numberOfComponents number of components
+   * @param temperature temperature in K
+   * @param pressure pressure in bara
+   */
+  private static void reinitSAFTOnPhase(PhasePCSAFT perturbedPhase, int numberOfComponents,
+      double temperature, double pressure) {
+    perturbedPhase.volInit();
+    for (int i = 0; i < numberOfComponents; i++) {
+      ((ComponentPCSAFT) perturbedPhase.getComponent(i)).initSAFTDerivatives(perturbedPhase,
+          numberOfComponents, temperature, pressure);
+    }
+  }
+
   /** {@inheritDoc} */
   @Override
   public double dFdN(PhaseInterface phase, int numberOfComponents, double temperature,
@@ -242,11 +291,7 @@ public double dFdNdT(PhaseInterface phase, int numberOfComponents, double temper
       ((ComponentPCSAFT) plus.getComponent(i)).init(plus.getTemperature(), pressure,
           plus.getNumberOfMolesInPhase(), 1.0, 0);
     }
-    plus.volInit();
-    for (int i = 0; i < numberOfComponents; i++) {
-      ((ComponentPCSAFT) plus.getComponent(i)).Finit(plus, plus.getTemperature(), pressure,
-          plus.getNumberOfMolesInPhase(), 1.0, numberOfComponents, 0);
-    }
+    reinitSAFTOnPhase(plus, numberOfComponents, plus.getTemperature(), pressure);
     double dFplus = ((ComponentPCSAFT) plus.getComponent(getComponentNumber())).dFdN(plus,
         numberOfComponents, plus.getTemperature(), pressure);
 
@@ -256,17 +301,86 @@ public double dFdNdT(PhaseInterface phase, int numberOfComponents, double temper
       ((ComponentPCSAFT) minus.getComponent(i)).init(minus.getTemperature(), pressure,
           minus.getNumberOfMolesInPhase(), 1.0, 0);
     }
-    minus.volInit();
-    for (int i = 0; i < numberOfComponents; i++) {
-      ((ComponentPCSAFT) minus.getComponent(i)).Finit(minus, minus.getTemperature(), pressure,
-          minus.getNumberOfMolesInPhase(), 1.0, numberOfComponents, 0);
-    }
+    reinitSAFTOnPhase(minus, numberOfComponents, minus.getTemperature(), pressure);
     double dFminus = ((ComponentPCSAFT) minus.getComponent(getComponentNumber())).dFdN(minus,
         numberOfComponents, minus.getTemperature(), pressure);
 
     return (dFplus - dFminus) / (2.0 * h);
   }
 
+  /**
+   * {@inheritDoc}
+   *
+   * <p>
+   * Computes the cross-derivative of the reduced residual Helmholtz energy with respect to moles of
+   * component i and total volume using numerical central difference. This derivative is essential
+   * for computing partial molar volumes and fugacity coefficient derivatives with respect to
+   * pressure.
+   * </p>
+   */
+  @Override
+  public double dFdNdV(PhaseInterface phase, int numberOfComponents, double temperature,
+      double pressure) {
+    double totalVolume = phase.getMolarVolume() * phase.getNumberOfMolesInPhase();
+    double h = totalVolume * 1.0e-6;
+    if (h < 1.0e-10) {
+      h = 1.0e-10;
+    }
+
+    PhasePCSAFT plus = (PhasePCSAFT) ((PhasePCSAFT) phase).clone();
+    plus.setMolarVolume((totalVolume + h) / phase.getNumberOfMolesInPhase());
+    reinitSAFTOnPhase(plus, numberOfComponents, temperature, pressure);
+    double dFdNplus = ((ComponentPCSAFT) plus.getComponent(getComponentNumber())).dFdN(plus,
+        numberOfComponents, temperature, pressure);
+
+    PhasePCSAFT minus = (PhasePCSAFT) ((PhasePCSAFT) phase).clone();
+    minus.setMolarVolume((totalVolume - h) / phase.getNumberOfMolesInPhase());
+    reinitSAFTOnPhase(minus, numberOfComponents, temperature, pressure);
+    double dFdNminus = ((ComponentPCSAFT) minus.getComponent(getComponentNumber())).dFdN(minus,
+        numberOfComponents, temperature, pressure);
+
+    return (dFdNplus - dFdNminus) / (2.0 * h);
+  }
+
+  /**
+   * {@inheritDoc}
+   *
+   * <p>
+   * Computes the second cross-derivative of the reduced residual Helmholtz energy with respect to
+   * moles of components i and j using numerical central difference at constant volume and
+   * temperature. This derivative is essential for Newton-Raphson flash convergence and stability
+   * analysis.
+   * </p>
+   */
+  @Override
+  public double dFdNdN(int j, PhaseInterface phase, int numberOfComponents, double temperature,
+      double pressure) {
+    double nj = phase.getComponent(j).getNumberOfMolesInPhase();
+    double h = Math.max(Math.abs(nj) * 1.0e-6, 1.0e-10);
+    double totalVolume = phase.getMolarVolume() * phase.getNumberOfMolesInPhase();
+    double totalMoles = phase.getNumberOfMolesInPhase();
+
+    PhasePCSAFT plus = (PhasePCSAFT) ((PhasePCSAFT) phase).clone();
+    plus.getComponent(j).setNumberOfMolesInPhase(nj + h);
+    double newTotalPlus = totalMoles + h;
+    plus.numberOfMolesInPhase = newTotalPlus;
+    plus.setMolarVolume(totalVolume / newTotalPlus);
+    reinitSAFTOnPhase(plus, numberOfComponents, temperature, pressure);
+    double dFdNplus = ((ComponentPCSAFT) plus.getComponent(getComponentNumber())).dFdN(plus,
+        numberOfComponents, temperature, pressure);
+
+    PhasePCSAFT minus = (PhasePCSAFT) ((PhasePCSAFT) phase).clone();
+    minus.getComponent(j).setNumberOfMolesInPhase(nj - h);
+    double newTotalMinus = totalMoles - h;
+    minus.numberOfMolesInPhase = newTotalMinus;
+    minus.setMolarVolume(totalVolume / newTotalMinus);
+    reinitSAFTOnPhase(minus, numberOfComponents, temperature, pressure);
+    double dFdNminus = ((ComponentPCSAFT) minus.getComponent(getComponentNumber())).dFdN(minus,
+        numberOfComponents, temperature, pressure);
+
+    return (dFdNplus - dFdNminus) / (2.0 * h);
+  }
+
   /**
    * <p>
    * calcF1dispSumTermdn.

src/main/java/neqsim/thermo/phase/PhasePCSAFTRahmat.java

@@ -19,6 +19,9 @@ public class PhasePCSAFTRahmat extends PhasePCSAFT {
   /** Logger object for class. */
   static Logger logger = LogManager.getLogger(PhasePCSAFTRahmat.class);
 
+  /** Cached molar volume from last converged solution for faster initial guess. */
+  private transient double cachedMolarVolume = -1.0;
+
   double dnSAFTdVdVdV = 1.0;
 
   double daHSSAFTdNdNdN = 1.0;
@@ -1142,7 +1145,13 @@ public double molarVolume(double pressure, double temperature, double A, double
     if (Btemp <= 0) {
       logger.info("b negative in volume calc");
     }
-    setMolarVolume(1.0 / BonV * Btemp / numberOfMolesInPhase);
+
+    // Use cached molar volume from previous converged solution as initial guess
+    if (cachedMolarVolume > 1.0e-10) {
+      setMolarVolume(cachedMolarVolume);
+    } else {
+      setMolarVolume(1.0 / BonV * Btemp / numberOfMolesInPhase);
+    }
     int iterations = 0;
     double oldMolarVolume = 0.0;
     // System.out.println("volume " + getVolume());
@@ -1203,6 +1212,7 @@ public double molarVolume(double pressure, double temperature, double A, double
      * " "+"  itert: " + iterations +" " + "  phase " + pt+ "  " + h + " " +dh + " B " + Btemp +
      * "  D " + Dtemp + " gv" + gV() + " fv " + fv() + " fvv" + fVV());
      */
+    cachedMolarVolume = getMolarVolume();
     return getMolarVolume();
   }
 }

src/main/java/neqsim/thermo/phase/PhasePCSAFTa.java

@@ -125,7 +125,42 @@ public double dFdT() {
   /** {@inheritDoc} */
   @Override
   public double dFdTdV() {
-    return super.dFdTdV();
+    return super.dFdTdV() + cpaon * dFCPAdTdV();
+  }
+
+  /**
+   * Computes the cross-derivative of the CPA association Helmholtz energy contribution with respect
+   * to temperature and volume using numerical central difference in volume.
+   *
+   * @return the CPA contribution to d2F/dTdV
+   */
+  public double dFCPAdTdV() {
+    double totalV = getTotalVolume();
+    double h = totalV * 1.0e-6;
+    if (h < 1.0e-10) {
+      h = 1.0e-10;
+    }
+    // Perturb volume up
+    PhasePCSAFTa plus = (PhasePCSAFTa) this.clone();
+    double molarVolPlus = (totalV + h) / numberOfMolesInPhase;
+    plus.setMolarVolume(molarVolPlus);
+    plus.volInit();
+    plus.solveX();
+    plus.hcpatot = plus.calc_hCPA();
+    plus.hcpatotdT = plus.calc_hCPAdT();
+    double dFdTplus = plus.dFCPAdT();
+
+    // Perturb volume down
+    PhasePCSAFTa minus = (PhasePCSAFTa) this.clone();
+    double molarVolMinus = (totalV - h) / numberOfMolesInPhase;
+    minus.setMolarVolume(molarVolMinus);
+    minus.volInit();
+    minus.solveX();
+    minus.hcpatot = minus.calc_hCPA();
+    minus.hcpatotdT = minus.calc_hCPAdT();
+    double dFdTminus = minus.dFCPAdT();
+
+    return (dFdTplus - dFdTminus) / (2.0 * h);
   }
 
   /** {@inheritDoc} */

src/main/java/neqsim/thermo/system/SystemPCSAFT.java

@@ -82,24 +82,33 @@ public SystemPCSAFT(double T, double P, boolean checkForSolids) {
   @Override
   public void addTBPfraction(String componentName2, double numberOfMoles, double molarMass,
       double density) {
-    // componentName = (componentName + "_" + getFluidName());
     super.addTBPfraction(componentName2, numberOfMoles, molarMass, density);
-    // addComponent(componentName2, numberOfMoles, 290.0, 30.0, 0.11);
     String componentName =
         getPhase(0).getComponent(getPhase(0).getNumberOfComponents() - 1).getComponentName();
-    for (int i = 0; i < numberOfPhases; i++) {
-      // getPhase(phaseIndex[i]).getComponent(componentName).setMolarMass(molarMass);
-      // getPhase(phaseIndex[i]).getComponent(componentName).setIsTBPfraction(true);
+    double mwG = molarMass * 1e3; // convert kg/mol to g/mol
+
+    // Gross & Sadowski (2001, I&EC Res, 40, 1244-1260) n-alkane correlations
+    double mSaft = 0.0257 * mwG + 0.8444;
+    double epskSaftm = 7.0106 * mwG + 117.10; // m * (eps/k) [K]
+    double msigm3 = 1.7296 * mwG + 15.050; // m * sigma^3 [Angstrom^3]
 
-      double mSaft = 0.0249 * molarMass * 1e3 + 0.9711;
-      double epskSaftm = 6.5446 * molarMass * 1e3 + 177.92;
-      double msigm = 1.6947 * molarMass * 1e3 + 23.27;
+    // Density-based correction for sigma when specific gravity is available.
+    // Higher SG (e.g. aromatics) implies smaller molecular diameter at same MW.
+    if (density > 0.5 && density < 1.5) {
+      // Approximate n-alkane specific gravity from Riazi-Daubert type fit (C5-C20)
+      double sgAlkane = 1.07 - 3.5649 / Math.pow(mwG, 0.25);
+      if (sgAlkane > 0.5 && sgAlkane < 1.0) {
+        msigm3 = msigm3 * (sgAlkane / density);
+      }
+    }
+
+    for (int i = 0; i < numberOfPhases; i++) {
       getPhase(phaseIndex[i]).getComponent(componentName).setmSAFTi(mSaft);
       getPhase(phaseIndex[i]).getComponent(componentName).setEpsikSAFT(epskSaftm / mSaft);
       getPhase(phaseIndex[i]).getComponent(componentName)
-          .setSigmaSAFTi(Math.pow(msigm / mSaft, 1.0 / 3.0) / 1.0e10);
-      logger.info("Saft parameters: m " + mSaft + " epsk " + epskSaftm / mSaft + " sigma "
-          + Math.pow(msigm / mSaft, 1.0 / 3.0));
+          .setSigmaSAFTi(Math.pow(msigm3 / mSaft, 1.0 / 3.0) / 1.0e10);
+      logger.info("PC-SAFT TBP params: m=" + mSaft + " eps/k=" + epskSaftm / mSaft + " sigma="
+          + Math.pow(msigm3 / mSaft, 1.0 / 3.0) + " Angstrom (SG=" + density + ")");
     }
   }
 
@@ -122,8 +131,8 @@ public SystemPCSAFT clone() {
    * </p>
    */
   public void commonInitialization() {
-    setImplementedCompositionDeriativesofFugacity(false);
-    setImplementedPressureDeriativesofFugacity(false);
-    setImplementedTemperatureDeriativesofFugacity(false);
+    setImplementedCompositionDeriativesofFugacity(true);
+    setImplementedPressureDeriativesofFugacity(true);
+    setImplementedTemperatureDeriativesofFugacity(true);
   }
 }

src/main/java/neqsim/thermo/system/SystemPCSAFTa.java

@@ -71,6 +71,17 @@ public SystemPCSAFTa(double T, double P, boolean checkForSolids) {
       phaseArray[numberOfPhases - 1].setRefPhase(phaseArray[1].getRefPhase());
     }
     this.useVolumeCorrection(false);
+    commonInitialization();
+  }
+
+  /**
+   * Common initialization for PC-SAFT with association. Enables fugacity derivative flags for
+   * Newton-Raphson flash convergence.
+   */
+  private void commonInitialization() {
+    setImplementedCompositionDeriativesofFugacity(true);
+    setImplementedPressureDeriativesofFugacity(true);
+    setImplementedTemperatureDeriativesofFugacity(true);
   }
 
   /** {@inheritDoc} */