diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index a8f3d6378f..af7f782aa4 100644 --- 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 index 582de57f78..ad61a4cf2b 100644 --- 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 + } + + } + /** *

* 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 index 0ffda8de69..3177abb004 100644 --- 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 new file mode 100644 index 0000000000..2ae7dcb7bb --- /dev/null +++ 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); + } + + +}