start adding more fluid property models (#1075)

* start adding more fluid property models * add java test to sdcode * further implementation * futher updates on Twu model * updated test * added PC, TC and VC
equinor · Aug 12, 2024 · 19251bb · 19251bb
1 parent e9eacdb
commit 19251bb
Show file tree

Hide file tree

Showing 4 changed files with 237 additions and 6 deletions.
diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
@@ -24,7 +24,8 @@
 				"ms-vscode.test-adapter-converter",
 				"shengchen.vscode-checkstyle",
 				"mechatroner.rainbow-csv",
-				"redhat.vscode-xml"
+				"redhat.vscode-xml",
+				"vscjava.vscode-java-test"
 			],
 			"settings": {
 				"java.configuration.updateBuildConfiguration": "interactive",

diff --git a/src/main/java/neqsim/thermo/characterization/TBPfractionModel.java b/src/main/java/neqsim/thermo/characterization/TBPfractionModel.java
@@ -104,6 +104,16 @@ public double calcCriticalViscosity(double molarMass, double density) {
           * 1e-7;
     }
 
+    @Override
+    public double calcRacketZ(SystemInterface thermoSystem, double molarMass, double density) {
+      throw new RuntimeException("calcm() method not defined");
+    }
+
+    @Override
+    public double calcm(double molarMass, double density) {
+      throw new RuntimeException("calcm() method not defined");
+    }
+
     @Override
     public boolean isCalcm() {
       return calcm;
@@ -154,13 +164,11 @@ public double calcPC(double molarMass, double density) {
     public double calcm(double molarMass, double density) {
       if (molarMass < 1120) {
         TBPfractionCoefs = TBPfractionCoefOil;
-        return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * molarMass
-            + TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.pow(molarMass, 2.0);
       } else {
         TBPfractionCoefs = TBPfractionCoefsHeavyOil;
-        return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * Math.log(molarMass)
-            + TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.sqrt(molarMass);
       }
+      return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * molarMass
+          + TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.pow(molarMass, 2.0);
     }
 
     @Override
@@ -312,6 +320,167 @@ public double calcAcentricFactor(double molarMass, double density) {
     }
   }
 
+  /**
+   * Lee-Kesler property estimation method
+   */
+  public class LeeKesler extends TBPBaseModel {
+    private static final long serialVersionUID = 1000;
+
+    public LeeKesler() {
+      calcm = false;
+    }
+
+    @Override
+    public double calcTC(double molarMass, double density) {
+      double sg = density;
+      double TB = calcTB(molarMass, density);
+      double TC =
+          189.8 + 450.6 * sg + (0.4244 + 0.1174 * sg) * TB + (0.1441 - 1.0069 * sg) * 1e5 / TB;
+      return TC;
+    }
+
+    @Override
+    public double calcPC(double molarMass, double density) {
+      double sg = density;
+      double TB = calcTB(molarMass, density);
+      double logpc =
+          3.3864 - 0.0566 / sg - ((0.43639 + 4.1216 / sg + 0.21343 / sg / sg) * 1e-3 * TB)
+              + ((0.47579 + 1.182 / sg + 0.15302 / sg / sg) * 1e-6 * TB * TB)
+              - ((2.4505 + 9.9099 / sg / sg) * 1e-10 * TB * TB * TB);
+      double PC = Math.exp(logpc) * 10;
+      return PC;
+    }
+
+    public double calcAcentricFactor(double molarMass, double density) {
+      return super.calcAcentricFactorKeslerLee(molarMass, density);
+    }
+  }
+
+  /**
+   * Two property estimation method
+   */
+  public class TwuModel extends TBPBaseModel {
+    private static final long serialVersionUID = 1000;
+
+    public TwuModel() {
+      calcm = false;
+    }
+
+    @Override
+    public double calcTC(double molarMass, double density) {
+      double sg = density;
+      double TB = calcTB(molarMass, density);
+      double MW = solveMW(TB);
+      double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
+          - 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
+      double phi = 1.0 - TB / Tcnalkane;
+      double SGalkane =
+          0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
+      double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
+          + 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
+      double VCnalkane = Math.pow(
+          (0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
+      double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
+      double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
+          + (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
+      double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
+      return TC;
+    }
+
+    public double calculateTfunc(double MW_alkane, double TB) {
+      double phi = Math.log(MW_alkane);
+      return Math
+          .exp(5.1264 + 2.71579 * phi - 0.28659 * phi * phi - 39.8544 / phi - 0.122488 / phi / phi)
+          - 13.7512 * phi + 19.6197 * phi * phi - TB;
+    }
+
+    public double computeGradient(double MW_alkane, double TB) {
+      double delta = 1;
+      double TfuncPlus = calculateTfunc(MW_alkane + delta, TB);
+      double TfuncMinus = calculateTfunc(MW_alkane - delta, TB);
+      return (TfuncPlus - TfuncMinus) / (2 * delta);
+    }
+
+    public double solveMW(double TB) {
+      double MW_alkane = TB / (5.8 - 0.0052 * TB);
+      double tolerance = 1e-6;
+      double prevMW_alkane;
+      double error = 1.0;
+      int iter = 0;
+
+      do {
+        iter++;
+        prevMW_alkane = MW_alkane;
+        double gradient = computeGradient(MW_alkane, TB);
+        MW_alkane -= 0.5 * calculateTfunc(MW_alkane, TB) / gradient;
+        error = Math.abs(MW_alkane - prevMW_alkane);
+      } while (Math.abs(error) > tolerance && iter < 1000 || iter < 3);
+
+      return MW_alkane;
+    }
+
+    @Override
+    public double calcPC(double molarMass, double density) {
+      double sg = density;
+      double TB = calcTB(molarMass, density);
+      double MW = solveMW(TB);
+      double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
+          - 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
+      double phi = 1.0 - TB / Tcnalkane;
+      double SGalkane =
+          0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
+      double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
+          + 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
+      double VCnalkane = Math.pow(
+          (0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
+      double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
+      double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
+          + (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
+      double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
+      double deltaSP = Math.exp(0.5 * (SGalkane - sg)) - 1.0;
+      double deltaSV = Math.exp(4.0 * (SGalkane * SGalkane - sg * sg)) - 1.0;
+      double fV = deltaSV
+          * (0.347776 * Math.pow(TB, -0.5) + (-0.182421 + 2.24890 * Math.pow(TB, -0.5)) * deltaSV);
+      double VC = VCnalkane * Math.pow(((1 + 2 * fV) / (1 - 2 * fV)), 2);
+      double fP = deltaSP * ((2.53262 - 34.4321 * Math.pow(TB, -0.5) - 0.00230193 * TB)
+          + (-11.4277 + 187.934 * Math.pow(TB, -0.5) + 0.00414963 * TB) * deltaSP);
+      double PC = PCnalkane * (TC / Tcnalkane) * (VCnalkane / VC)
+          * Math.pow(((1 + 2 * fP) / (1 - 2 * fP)), 2);
+      return PC * 10.0; // * 10 due to conversion MPa to bar
+    }
+
+    @Override
+    public double calcCriticalVolume(double molarMass, double density) {
+      double sg = density;
+      double TB = calcTB(molarMass, density);
+      double MW = solveMW(TB);
+      double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
+          - 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
+      double phi = 1.0 - TB / Tcnalkane;
+      double SGalkane =
+          0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
+      double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
+          + 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
+      double VCnalkane = Math.pow(
+          (0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
+      double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
+      double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
+          + (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
+      double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
+      double deltaSP = Math.exp(0.5 * (SGalkane - sg)) - 1.0;
+      double deltaSV = Math.exp(4.0 * (SGalkane * SGalkane - sg * sg)) - 1.0;
+      double fV = deltaSV
+          * (0.347776 * Math.pow(TB, -0.5) + (-0.182421 + 2.24890 * Math.pow(TB, -0.5)) * deltaSV);
+      double VC = VCnalkane * Math.pow(((1 + 2 * fV) / (1 - 2 * fV)), 2);
+      double fP = deltaSP * ((2.53262 - 34.4321 * Math.pow(TB, -0.5) - 0.00230193 * TB)
+          + (-11.4277 + 187.934 * Math.pow(TB, -0.5) + 0.00414963 * TB) * deltaSP);
+      double PC = PCnalkane * (TC / Tcnalkane) * (VCnalkane / VC)
+          * Math.pow(((1 + 2 * fP) / (1 - 2 * fP)), 2);
+      return VC * 1e3; // m3/mol
+    }
+
+  }
+
   /**
    * <p>
    * getModel.
@@ -334,6 +503,10 @@ public TBPModelInterface getModel(String name) {
       return new PedersenTBPModelPRHeavyOil();
     } else if (name.equals("RiaziDaubert")) {
       return new RiaziDaubert();
+    } else if (name.equals("Lee-Kesler")) {
+      return new LeeKesler();
+    } else if (name.equals("Twu")) {
+      return new TwuModel();
     } else {
       // System.out.println("not a valid TBPModelName.................");
       return new PedersenTBPModelSRK();

diff --git a/src/main/java/neqsim/thermo/system/SystemThermo.java b/src/main/java/neqsim/thermo/system/SystemThermo.java
@@ -869,7 +869,9 @@ public void addTBPfraction(String componentName, double numberOfMoles, double mo
       molarMass = 1000 * molarMass;
       TC = characterization.getTBPModel().calcTC(molarMass, density);
       PC = characterization.getTBPModel().calcPC(molarMass, density);
-      m = characterization.getTBPModel().calcm(molarMass, density);
+      if (characterization.getTBPModel().isCalcm()) {
+        m = characterization.getTBPModel().calcm(molarMass, density);
+      }
       acs = characterization.getTBPModel().calcAcentricFactor(molarMass, density);
       // TBPfractionCoefs[2][0]+TBPfractionCoefs[2][1]*molarMass+TBPfractionCoefs[2][2]*density+TBPfractionCoefs[2][3]*Math.pow(molarMass,2.0);
       TB = characterization.getTBPModel().calcTB(molarMass, density);

diff --git a/src/test/java/neqsim/thermo/characterization/TBPfractionModelTest.java b/src/test/java/neqsim/thermo/characterization/TBPfractionModelTest.java
@@ -0,0 +1,55 @@
+package neqsim.thermo.characterization;
+
+import static org.junit.jupiter.api.Assertions.assertEquals;
+import org.junit.jupiter.api.Test;
+import neqsim.thermo.system.SystemInterface;
+import neqsim.thermo.system.SystemPrEos;
+import neqsim.thermo.system.SystemSrkEos;
+
+public class TBPfractionModelTest {
+  @Test
+  void testTwuModel() {
+    SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
+    thermoSystem.getCharacterization().setTBPModel("Twu");
+    thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
+    assertEquals(536.173400, thermoSystem.getComponent(0).getTC(), 1e-3);
+    assertEquals(26.52357312690, thermoSystem.getComponent(0).getPC(), 1e-3);
+    assertEquals(0.56001213933, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
+    assertEquals(437.335493, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
+  }
+
+  @Test
+  void testLeeKeslerModel() {
+    SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
+    thermoSystem.getCharacterization().setTBPModel("Lee-Kesler");
+    thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
+    assertEquals(562.4229803010662, thermoSystem.getComponent(0).getTC(), 1e-3);
+    assertEquals(28.322987349048354, thermoSystem.getComponent(0).getPC(), 1e-3);
+    assertEquals(0.3509842412742902, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
+    assertEquals(427.99744457199, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
+  }
+
+  @Test
+  void testPedersenSRKModel() {
+    SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
+    thermoSystem.getCharacterization().setTBPModel("PedersenSRK");
+    thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
+    assertEquals(554.3185637098962, thermoSystem.getComponent(0).getTC(), 1e-3);
+    assertEquals(26.007549082822628, thermoSystem.getComponent(0).getPC(), 1e-3);
+    assertEquals(0.508241, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
+    assertEquals(384.6714299777243, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
+  }
+
+  @Test
+  void testPedersenPRModel() {
+    SystemInterface thermoSystem = new SystemPrEos(298.0, 10.0);
+    thermoSystem.getCharacterization().setTBPModel("PedersenPR");
+    thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
+    assertEquals(560.546, thermoSystem.getComponent(0).getTC(), 1e-3);
+    assertEquals(25.838137535018557, thermoSystem.getComponent(0).getPC(), 1e-3);
+    assertEquals(0.3838836222383, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
+    assertEquals(405.0890245138075, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
+  }
+
+
+}