1078 implement gamma distribution model for plus fraction characteriz…

…ation (#1086) * started inital implementation * feat: V1 Gamma model * feat: last update * some work for verification of code * added get and set gamma model parameters --------- Co-authored-by: Even Solbraa <[email protected]>
equinor · Sep 4, 2024 · 3b7bb38 · 3b7bb38
1 parent b08e894
commit 3b7bb38
Show file tree

Hide file tree

Showing 6 changed files with 262 additions and 49 deletions.
diff --git a/src/main/java/neqsim/thermo/characterization/Characterise.java b/src/main/java/neqsim/thermo/characterization/Characterise.java
@@ -156,7 +156,7 @@ public void characterisePlusFraction() {
     if (plusFractionModel.hasPlusFraction()) {
       if (plusFractionModel.getMPlus() > plusFractionModel.getMaxPlusMolarMass()) {
         logger.error("plus fraction molar mass too heavy for " + plusFractionModel.getName());
-        plusFractionModel = plusFractionModelSelector.getModel("heavyOil");
+        plusFractionModel = plusFractionModelSelector.getModel("Pedersen Heavy Oil");
         logger.info("changing to " + plusFractionModel.getName());
       }
       plusFractionModel.characterizePlusFraction(TBPfractionModel);

diff --git a/src/main/java/neqsim/thermo/characterization/LumpingModel.java b/src/main/java/neqsim/thermo/characterization/LumpingModel.java
@@ -306,6 +306,7 @@ public void generateLumpedComposition(Characterise charac) {
         system.removeComponent(system.getPhase(0)
             .getComponent(charac.getPlusFractionModel().getPlusComponentNumber()).getName());
       }
+      system.init(0);
     }
 
     @Override
@@ -398,18 +399,18 @@ public void generateLumpedComposition(Characterise charac) {
         pseudoNumber++;
         String addName = "C" + Integer.toString(starti) + "-" + Integer.toString(i);
         getLumpedComponentNames()[k] = addName;
-        // System.out.println("adding " + addName);
+        //System.out.println("adding " + addName);
         fractionOfHeavyEnd[k] = zPlus[k] / molFracTot;
-
         system.addTBPfraction(addName, totalNumberOfMoles * zPlus[k], Maverage / zPlus[k],
             denstemp1 / denstemp2);
         k++;
-        starti = i + 1;
+        starti = i;
       }
       if (charac.getPlusFractionModel().hasPlusFraction()) {
         system.removeComponent(system.getPhase(0)
             .getComponent(charac.getPlusFractionModel().getPlusComponentNumber()).getName());
       }
+      system.init(0);
     }
 
     @Override

diff --git a/src/main/java/neqsim/thermo/characterization/PlusFractionModel.java b/src/main/java/neqsim/thermo/characterization/PlusFractionModel.java
@@ -227,9 +227,8 @@ public void characterizePlusFraction(TBPModelInterface TBPModel) {
         z[i] = Math.exp(getCoef(0) + getCoef(1) * i);
         M[i] = PVTsimMolarMass[i - 6] / 1000.0;
         dens[i] = getCoef(2) + getCoef(3) * Math.log(i);
-
-        // System.out.println("z,m,dens " + z[i] + " " + M[i] + " " + dens[i]);
       }
+      // System.out.println("z,m,dens " + z[i] + " " + M[i] + " " + dens[i]);
     }
 
     @Override
@@ -287,6 +286,11 @@ public void setCoefs(double[] coefs) {
     public void setCoefs(double coef, int i) {
       this.coefs[i] = coef;
     }
+
+    @Override
+    public void setLastPlusFractionNumber(int fract) {
+      lastPlusFractionNumber = fract + 1;
+    }
   }
 
   private class PedersenHeavyOilPlusModel extends PedersenPlusModel {
@@ -299,6 +303,179 @@ public PedersenHeavyOilPlusModel() {
     }
   }
 
+  class WhitsonGammaModel extends PedersenPlusModel {
+
+    public double[] zValues;
+    public double[] molarMasses;
+    public double[] densities;
+    public double eta = 90; // minimum molecular weight in C7+
+    public String model = "Whitson";
+    public double alfa = 1.0;
+    public double betta = Double.NaN;
+
+    public WhitsonGammaModel() {
+      name = "Whitson Gamma";
+    }
+
+    public void setCalculationModel(String model) {
+      this.model = model;
+    }
+
+    public void characterizePlusFractionWhitsonGamma() {
+
+    }
+
+    public double gamma(double X) {
+      double[] dataB = {-0.577191652, 0.988205891, -0.897056937, 0.918206857, -0.756704078,
+          0.482199394, -0.193527818, 0.035868343};
+      double const_ = 1.0;
+      double XX = X;
+      if (X < 1.0) {
+        XX = X + 1.0;
+      }
+      while (XX >= 2.0) {
+        XX -= 1.0;
+      }
+      const_ = XX * const_;
+      XX -= 1.0;
+      double Y = 1.0;
+      for (int i = 1; i <= 8; i++) {
+        Y += dataB[i - 1] * XX * i;
+      }
+      double GAMMA = const_ * Y;
+      if (X < 1.0) {
+        GAMMA /= X;
+      }
+      return GAMMA;
+    }
+
+    public double[] P0P1(double MWB) {
+      double P0 = 0.0;
+      double P1 = 0.0;
+      if (MWB == eta) {
+        return new double[] {P0, P1};
+      }
+      double Y = (MWB - eta) / betta;
+      double Q = Math.exp(-Y) * Math.pow(Y, alfa) / gamma(alfa);
+      double TERM = 1.0 / alfa;
+      double S = TERM;
+      for (int j = 1; j <= 10000; j++) {
+        TERM *= Y / (alfa + j);
+        S += TERM;
+        if (Math.abs(TERM) <= 1e-8) {
+          P0 = Q * S;
+          P1 = Q * (S - 1.0 / alfa);
+          break;
+        }
+
+      }
+      return new double[] {P0, P1};
+    }
+
+    public void densityUOP() {
+      // Calculates density of the C+ function with Watson or Universal Oil Products
+      // characterization factor
+      // Experience has shown, that the method is not very accurate for C20+
+      double Kw = 4.5579 * Math.pow(MPlus * 1000, 0.15178) * Math.pow(densPlus, -1.18241);
+      for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
+        densities[i - 1] = 6.0108 * Math.pow(molarMasses[i - 1], 0.17947) * Math.pow(Kw, -1.18241);
+      }
+    }
+
+
+
+    @Override
+    public void characterizePlusFraction(TBPModelInterface TBPModel) {
+      system.init(0);
+      double MWBU = Double.NaN;
+      double MWBL = Double.NaN;
+      double sumZ = 0.0;
+
+      betta = (MPlus * 1000 - eta) / alfa;
+      // Implement the Gamma distribution for the plus fraction
+      zValues = new double[lastPlusFractionNumber];
+      molarMasses = new double[lastPlusFractionNumber];
+      densities = new double[lastPlusFractionNumber];
+
+      if (model.equals("Whitson")) {
+        for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
+          if (i == 1) {
+            MWBU = eta;
+          }
+          MWBL = MWBU;
+          MWBU = MWBL + 14;
+          if (i == lastPlusFractionNumber) {
+            MWBU = 10000.0;
+          }
+          double[] P0LP1L = P0P1(MWBL);
+          double P0L = P0LP1L[0];
+          double P1L = P0LP1L[1];
+
+          double[] P0UP1U = P0P1(MWBU);
+          double P0U = P0UP1U[0];
+          double P1U = P0UP1U[1];
+
+          double Z = P0U - P0L;
+          if (Z < 1E-15) {
+            Z = 1E-15;
+          }
+          zValues[i] = Z * zPlus;
+          double MWAV = eta + alfa * betta * (P1U - P1L) / (P0U - P0L);
+          molarMasses[i] = MWAV / 1000;
+          sumZ = sumZ + zValues[i];
+        }
+        densityUOP();
+
+
+      }
+
+      // Normalize z values to ensure sumZ equals zPlus
+      for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
+        zValues[i] *= zPlus / sumZ;
+      }
+    }
+
+    public void setGammaParameters(double shape, double minMW) {
+      this.alfa = shape;
+      this.eta = minMW;
+    }
+
+    public double[] getGammaParameters() {
+      return new double[] {this.alfa, this.eta};
+    }
+
+    @Override
+    public double[] getCoefs() {
+      return new double[] {alfa, eta};
+    }
+
+    @Override
+    public double getCoef(int i) {
+      if (i == 0)
+        return alfa;
+      if (i == 1)
+        return eta;
+      return 0;
+    }
+
+    @Override
+    public double[] getZ() {
+      return zValues;
+    }
+
+    @Override
+    public double[] getM() {
+      return molarMasses;
+    }
+
+    @Override
+    public double[] getDens() {
+      return densities;
+    }
+
+
+  }
+
   /**
    * <p>
    * getModel.
@@ -312,6 +489,8 @@ public PlusFractionModelInterface getModel(String name) {
       return new PedersenPlusModel();
     } else if (name.equals("Pedersen Heavy Oil")) {
       return new PedersenHeavyOilPlusModel();
+    } else if (name.equals("Whitson Gamma Model")) {
+      return new WhitsonGammaModel();
     } else {
       return new PedersenPlusModel();
     }

diff --git a/src/main/java/neqsim/thermo/characterization/PlusFractionModelInterface.java b/src/main/java/neqsim/thermo/characterization/PlusFractionModelInterface.java
@@ -72,6 +72,8 @@ public interface PlusFractionModelInterface extends java.io.Serializable {
    */
   public int getLastPlusFractionNumber();
 
+  public void setLastPlusFractionNumber(int fract);
+
   /**
    * <p>
    * getPlusComponentNumber.

diff --git a/src/test/java/neqsim/thermo/characterization/CharacteriseTest.java b/src/test/java/neqsim/thermo/characterization/CharacteriseTest.java
@@ -4,6 +4,7 @@
 import org.junit.jupiter.api.Test;
 import neqsim.thermo.system.SystemInterface;
 import neqsim.thermo.system.SystemSrkEos;
+import neqsim.thermodynamicOperations.ThermodynamicOperations;
 
 public class CharacteriseTest extends neqsim.NeqSimTest {
   static SystemInterface thermoSystem = null;
@@ -33,9 +34,37 @@ void testCharacterisePlusFraction() {
     thermoSystem.addTBPfraction("C8", 1.0, 120.0 / 1000.0, 0.76);
     thermoSystem.addTBPfraction("C9", 1.0, 140.0 / 1000.0, 0.79);
     thermoSystem.addPlusFraction("C10", 11.0, 290.0 / 1000.0, 0.82);
-    thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12);
+    //thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(6);
+    //thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel");
+    thermoSystem.getCharacterization().setLumpingModel("no lumping");
+    thermoSystem.getCharacterization().characterisePlusFraction();
+   // assertEquals(15, thermoSystem.getNumberOfComponents());
+    thermoSystem.prettyPrint();
+  }
+
+
+  @Test
+  void testCharacterisePlusFractionGAMMA() {
+    thermoSystem = new SystemSrkEos(298.0, 10.0);
+    thermoSystem.addComponent("methane", 51.0);
+    thermoSystem.addPlusFraction("C10", 11.0, 290.0 / 1000.0, 0.82);
+    thermoSystem.getCharacterization().setPlusFractionModel("Whitson Gamma Model");
     thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel");
+    thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(15);
     thermoSystem.getCharacterization().characterisePlusFraction();
-    assertEquals(15, thermoSystem.getNumberOfComponents());
+    // logger.info("number of components " + thermoSystem.getNumberOfComponents());
+    //assertEquals(86, thermoSystem.getNumberOfComponents());
+    //System.out.println(thermoSystem.getComponent("C1-2_PC").getz());
+    thermoSystem.prettyPrint();
+
+    thermoSystem.setPressure(1, "bara");
+
+
+    ThermodynamicOperations thermoOps = new ThermodynamicOperations(thermoSystem); 
+    thermoOps.TPflash();
+
+    thermoSystem.initProperties();
+    thermoSystem.prettyPrint();
   }
+
 }