adding a better water model

Author: RealAntEngineer

src/main/java/com/rae/formicapi/thermal_utilities/PengRobinsonEOS.java

@@ -4,7 +4,6 @@
 import com.rae.formicapi.FormicAPI;
 import com.rae.formicapi.math.Solvers;
 import net.createmod.catnip.data.Couple;
-import net.minecraft.util.Mth;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.ArrayList;
@@ -37,6 +36,15 @@ public double alpha(double T) {
         return Math.pow(1 + kappa * (1 - Math.sqrt(Tr)), 2);
     }
 
+    private double dAlpha_dT(double T) {
+        double sqrtTr = Math.sqrt(T / Tc);
+        return -kappa * (1 + kappa * (1 - sqrtTr)) / (Tc * sqrtTr);
+    }
+
+    private double da_dT(double T) {
+        return a0 * dAlpha_dT(T); // ac is the "a" at critical temperature
+    }
+
     public double a(double T) {
         return a0 * alpha(T);
     }
@@ -80,12 +88,14 @@ public double volumeMolar(double T, double P, double vaporFraction) {
 
     /**
      *  What does this represent exactly ?
-     * @param T the temperature in Kelvin
-     * @param P the pressure in Pascals
-     * @return the possible compressibility coef for the given isotherm and Pressure : minimal is the value at the saturation curve for liquid
-     * the middle on is the normal on that is non physical, and the max is for the vapor
+     * @param T the temperature in Kelvin, needs to be strictly positive
+     * @param P the pressure in Pascals, needs to be strictly positive
+     * @return the possible compressibility coefficient for the given isotherm and Pressure : minimal is the value at the saturation curve for liquid
+     * the middle is non-physical, and the max is for the vapor
      */
     public double @NotNull [] getZFactors(double T, double P) {
+        assert T > 0;
+        assert P > 0;
         double A = a(T) * P / (R * R * T * T);
         double B = b * P / (R * T);
 
@@ -101,6 +111,7 @@ public double volumeMolar(double T, double P, double vaporFraction) {
                 .sorted()
                 .toArray();
     }
+    //not robust at low pressure
     public double fugacityCoefficient(double T, double P, double Z) {
         double R = PengRobinsonEOS.R;
         double a = a(T); // temperature-dependent 'a'
@@ -114,6 +125,8 @@ public double fugacityCoefficient(double T, double P, double Z) {
         double term3 = term2 * Math.log(logArgument);
 
         double lnPhi = term1 - term3;
+
+        System.out.println("fug : "+Math.exp(lnPhi));
         return Math.exp(lnPhi);
     }
 
@@ -141,7 +154,7 @@ public List<Double> findSpinodalPoints(double T, double vMin, double vMax) {
         return spinodals;
     }
 
-    public double computeCoexistencePressure(double T) {
+    public double computeSaturationPressure(double T) {
         if (T > Tc) {
             throw new IllegalArgumentException("T must be lower than the critical temperature to compute a coexistence pressure");
         }
@@ -157,67 +170,92 @@ public double computeCoexistencePressure(double T) {
         assert spinodals.size() == 2;
 
         double pMin = Math.max(1e-3, pressure(T, spinodals.get(0)));  // lower pressure bound
-        double pMax = pressure(T, spinodals.get(1));                 // upper pressure bound
+        double pMax = pressure(T, spinodals.get(1));               // upper pressure bound
 
         // Fugacity difference function
-        Function<Float, Float> fugacityDifference = (Float PGuess) -> {
-            double P = PGuess;
-            if (P <= 0) return Float.NaN;
-
-            double[] ZRoots = getZFactors(T, P);
-            if (ZRoots.length < 2) return Float.NaN; // Need at least two roots for coexistence
-
-            double Zl = Arrays.stream(ZRoots).min().getAsDouble();
-            double Zv = Arrays.stream(ZRoots).max().getAsDouble();
-
-            double phiL = fugacityCoefficient(T, P, Zl);
-            double phiV = fugacityCoefficient(T, P, Zv);
-            if (phiL <= 0 || phiV <= 0 || Double.isNaN(phiL) || Double.isNaN(phiV)) return Float.NaN;
+        Function<Float, Float> areaDifferenceFunction = (Float Guess) -> {
+            double P = Guess;
+            double[] roots = getZFactors(T, P);
+            if (roots.length < 2) return Float.MAX_VALUE;
+
+            double Zl = Arrays.stream(roots).min().getAsDouble();
+            double Zv = Arrays.stream(roots).max().getAsDouble();
+            double Vl = Zl * PengRobinsonEOS.R * T / P;
+            double Vv = Zv * PengRobinsonEOS.R * T / P;
+            if (Math.abs(Vv - Vl) < 1e-9) return Float.MAX_VALUE;
+
+            int n = 1000;
+            double dV = (Vv - Vl) / n;
+            double area = 0;
+            for (int i = 0; i < n; i++) {
+                double V1 = Vl + i * dV;
+                double V2 = V1 + dV;
+                double p1 = pressure(T, V1);
+                double p2 = pressure(T, V2);
+                area += 0.5 * (p1 + p2) * dV;
+            }
 
-            double diff = Math.abs(Math.log(phiL) - Math.log(phiV));
-            return (float) diff;
+            double rect = P * (Vv - Vl);
+            return (float)Math.abs(area - rect);
+        };
+        Function<Float, Float> fugacityDifference = (Float P) -> {
+            try {
+                double[] Zroots = getZFactors(T, P);
+                if (Zroots.length < 2) return Float.MAX_VALUE;
+
+                double Zl = Arrays.stream(Zroots).min().getAsDouble();
+                double Zv = Arrays.stream(Zroots).max().getAsDouble();
+
+                double phiL = fugacityCoefficient(Zl, T, P);
+                double phiV = fugacityCoefficient(Zv, T, P);
+                return (float) Math.abs(Math.log(phiV) - Math.log(phiL));
+            } catch (Exception e) {
+                return Float.MAX_VALUE;
+            }
         };
 
-        float PCoexistence = Solvers.dichotomy(fugacityDifference, (float) pMin, (float) pMax, 1e-3f);
-
-        if (Float.isNaN(PCoexistence)) {
-            throw new RuntimeException("Dichotomy failed: no valid fugacity-based root found.");
-        }
-
-        return PCoexistence;
+        return Solvers.gradientDecent(areaDifferenceFunction, (float)Math.max(1e-3,(pMin+pMax)/2), 10,1e-3f,1e-4f);
     }
 
     /**
-     *
+     * Waarning only use this if you know that the fuild is in the LV phase.
      * @param T temperature
      * @return return the couple Vl, Vv or the intermediate V if it's super critical.
      */
+    @Deprecated
     public Couple<Double> getSaturationVolumes(double T) {
         try {
-            double Pcoex = computeCoexistencePressure(T); // try normal coexistence
-            double Vl = volumeMolar(T, Pcoex, 0.0); // saturated liquid
-            double Vv = volumeMolar(T, Pcoex, 1.0); // saturated vapor
-            return Couple.create(Vl, Vv);
+            if (T < Tc ) {
+
+                double PSat = computeSaturationPressure(T); // try normal coexistence
+                double Vl = volumeMolar(T, PSat, 0.0); // saturated liquid
+                double Vv = volumeMolar(T, PSat, 1.0); // saturated vapor
+                return Couple.create(Vl, Vv);
+            } else {
+                double[] critRoot = getZFactors(Tc, Pc);
+
+                // Start guess: approximate ideal gas volume
+                double Z = Arrays.stream(critRoot).max().getAsDouble();
+                //we need to find a better interpolation function.
+                double VmInflection = Z * R * Tc / Pc;
+                return Couple.create(VmInflection, VmInflection);
+            }
         } catch (RuntimeException e) {
-            FormicAPI.LOGGER.debug("found outside LV phase {}", T);
+            FormicAPI.LOGGER.debug("Unexpectedly found outside LV phase {}", T);
             FormicAPI.LOGGER.debug(e);
             // Fall back to inflection point
-            double[] critRoot = getZFactors(Tc, Pc);
-
-            // Start guess: approximate ideal gas volume
-            double Z = Arrays.stream(critRoot).max().getAsDouble();
-            //we need to find a better interpolation function.
-            double VmInflection = Z * R * Tc / Pc;
-            return Couple.create(VmInflection, VmInflection);
+            return Couple.create(Double.NaN, Double.NaN);
         }
     }
 
-    public double idealGasEntropy(double T) {
-        // Ideal gas: s = ∫(Cv/T)dT + R ln(P_ref/P)
-        // Approximate Cv/R = constant for now, e.g. 3.5 for a monatomic, ~2.0 for H2O vapor
-        double Cv = 2.0 * R;
-        return Cv * Math.log(T); // Ignore pressure-dependent term (constant reference)
+    public double idealGasEntropy(double T, double P) {
+        double Cp = 3.5 * R;
+        double Pref = 101325.0;
+        double Tref = 298.15;
+
+        return Cp * Math.log(T / Tref) - R * Math.log(P / Pref);
     }
+
     public double entropyDeparturePR(double T, double Vm) {
         double a = a(T);  // temperature-dependent 'a' from PR EOS
         double sqrt2 = Math.sqrt(2);
@@ -230,8 +268,7 @@ public double entropyDeparturePR(double T, double Vm) {
     }
     public double totalEntropy(double T, double Vm) {
         // Ideal gas entropy (reference + integration)
-        double sIdeal = idealGasEntropy(T) - idealGasEntropy(298.15);  // reference s = 0 at 298.15 K
-
+        double sIdeal = idealGasEntropy(T, pressure(T, Vm));
         // Entropy departure from Peng-Robinson EOS
         double sDeparture = entropyDeparturePR(T, Vm);
 
@@ -246,30 +283,36 @@ public double totalEntropy(double T, double Vm) {
      */
     private double idealGasEnthalpy(double T) {
         double Cp = 75.0; // J/mol·K
-        return Cp * T;
+        return Cp * (T - 273.15);
     }
 
     /**
      *
      * @param T temperature
      * @param Vm specific volume
-     * @return residual Enthalpy
+     * @return the residual enthalpy (H^res) at given temperature and molar volume.
+     * Units: J/mol
      */
     private double residualEnthalpy(double T, double Vm) {
-        double a = a(T);
+        double a = a(T); // attraction parameter at T
+        double da_dT = da_dT(T); // temperature derivative of a
         double P = pressure(T, Vm);
         double Z = P * Vm / (R * T);
 
-        double logTerm = Math.log((Vm + (1 + Math.sqrt(2)) * b) / (Vm + (1 - Math.sqrt(2)) * b));
-        double hRes = T * (Z - 1) * R - a / (b * Math.sqrt(8)) * logTerm;
-        if (Double.isNaN(hRes)) {
-            System.out.println("log therm is NaN for T "+T + " and Vm "+Vm);
-        }
+        // PR-specific constants
+        double sqrt2 = Math.sqrt(2);
+        double B = b;
+        double logTerm = Math.log((Vm + (1 + sqrt2) * B) / (Vm + (1 - sqrt2) * B));
+
+        // Correct residual enthalpy formula
+        double hRes = R * T * (Z - 1) + (T * da_dT - a) / (B * sqrt2) * logTerm;
+
         return hRes;
     }
 
+
     public double totalEnthalpy(double T, double Vm) {
-        return idealGasEnthalpy(T) + residualEnthalpy(T, Vm);
+        return idealGasEnthalpy(Math.max(T,1)) + residualEnthalpy(Math.max(T,1), Math.max(Vm,b));
     }
 
     public double getB() {