From 3b7bb38f799d6599638c1d971522adf63259f299 Mon Sep 17 00:00:00 2001 From: Sviatoslav Eroshkin <109044598+Sviatose@users.noreply.github.com> Date: Wed, 4 Sep 2024 17:57:49 +0200 Subject: [PATCH] 1078 implement gamma distribution model for plus fraction characterization (#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 <41290109+EvenSol@users.noreply.github.com> --- .../thermo/characterization/Characterise.java | 2 +- .../thermo/characterization/LumpingModel.java | 7 +- .../characterization/PlusFractionModel.java | 183 +++++++++++++++++- .../PlusFractionModelInterface.java | 2 + .../characterization/CharacteriseTest.java | 33 +++- .../PlusFractionModelTest.java | 84 ++++---- 6 files changed, 262 insertions(+), 49 deletions(-) diff --git a/src/main/java/neqsim/thermo/characterization/Characterise.java b/src/main/java/neqsim/thermo/characterization/Characterise.java index 34cd86f10a..f46f425a90 100644 --- 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 index b10cdae9e9..56143069bc 100644 --- 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 index 6c6e6eb259..2446dea324 100644 --- 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; + } + + + } + /** *

* 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 index df08e6fa83..bfe437d1ce 100644 --- 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); + /** *

* getPlusComponentNumber. diff --git a/src/test/java/neqsim/thermo/characterization/CharacteriseTest.java b/src/test/java/neqsim/thermo/characterization/CharacteriseTest.java index a7c5f9bf0b..e17c53b8a0 100644 --- 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(); } + } diff --git a/src/test/java/neqsim/thermo/characterization/PlusFractionModelTest.java b/src/test/java/neqsim/thermo/characterization/PlusFractionModelTest.java index 9aeb541864..6cdbe39d68 100644 --- a/src/test/java/neqsim/thermo/characterization/PlusFractionModelTest.java +++ b/src/test/java/neqsim/thermo/characterization/PlusFractionModelTest.java @@ -2,6 +2,7 @@ import static org.junit.jupiter.api.Assertions.assertEquals; import org.junit.jupiter.api.Test; +import neqsim.thermo.characterization.PlusFractionModel.WhitsonGammaModel; import neqsim.thermo.system.SystemInterface; import neqsim.thermo.system.SystemSrkEos; import neqsim.thermodynamicOperations.ThermodynamicOperations; @@ -32,29 +33,27 @@ void testPedersenPlusModelCharacterization() { * component */ thermoSystem.getCharacterization().setPlusFractionModel("Pedersen"); - thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel"); // this is default // lumping model in // neqsim. Needs to be // set before calling // characterisePlusFraction() - thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12); // specify - // numer - // of - // lumped - // components - // (C6-C80 - // components) + thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(9); // specif + // numer + // of + // lumped + // components + // (C6-C80 + // components) thermoSystem.getCharacterization().characterisePlusFraction(); - assertEquals(16, thermoSystem.getNumberOfComponents()); + assertEquals(17, thermoSystem.getNumberOfComponents()); ThermodynamicOperations ops = new ThermodynamicOperations(thermoSystem); ops.TPflash(); - + // thermoSystem.prettyPrint(); assertEquals(0.76652495787, thermoSystem.getBeta(), 1e-4); - } @Test @@ -87,23 +86,22 @@ void testPedersenHeavyOilPlusModelCharacterization() { // neqsim. Needs to be // set before calling // characterisePlusFraction() - - thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12); // specify - // numer - // of - // lumped - // components - // (C6-C80 - // components) + thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(3); + thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(8); + // specify + // numer + // of + // lumped + // components + // (C6-C80 + // components) thermoSystem.getCharacterization().characterisePlusFraction(); - assertEquals(16, thermoSystem.getNumberOfComponents()); + assertEquals(12, thermoSystem.getNumberOfComponents()); ThermodynamicOperations ops = new ThermodynamicOperations(thermoSystem); ops.TPflash(); - - assertEquals(0.767085187, thermoSystem.getBeta(), 1e-4); - - + // thermoSystem.prettyPrint(); + assertEquals(0.767272255056255, thermoSystem.getBeta(), 1e-4); } @Test @@ -115,9 +113,7 @@ void testGammaModelCharacterization() { thermoSystem.addComponent("methane", 51.0); thermoSystem.addComponent("ethane", 1.0); thermoSystem.addComponent("propane", 1.0); - - thermoSystem.getCharacterization().setTBPModel("PedersenSRK"); // this need to be set before - // adding oil components + thermoSystem.getCharacterization().setTBPModel("PedersenSRK"); thermoSystem.addTBPfraction("C6", 1.0, 90.0 / 1000.0, 0.7); thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73); @@ -125,24 +121,30 @@ void testGammaModelCharacterization() { thermoSystem.addTBPfraction("C9", 1.0, 140.0 / 1000.0, 0.79); thermoSystem.addPlusFraction("C10", 11.0, 290.0 / 1000.0, 0.82); - thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel"); // this is default - // lumping model in - // neqsim. Needs to be - // set before calling - // characterisePlusFraction() + thermoSystem.getCharacterization().setPlusFractionModel("Whitson Gamma Model"); + + // Add how to set parameters in the gamma model here + + thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12); - thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12); // specify - // numer - // of - // lumped - // components - // (C6-C80 - // components) thermoSystem.getCharacterization().characterisePlusFraction(); + thermoSystem.setMixingRule("classic"); assertEquals(16, thermoSystem.getNumberOfComponents()); - } - + ThermodynamicOperations ops = new ThermodynamicOperations(thermoSystem); + ops.TPflash(); + // thermoSystem.prettyPrint(); + assertEquals(0.746485111, thermoSystem.getBeta(), 1e-4); + + // illustration of how to set parameters for the gamma model + ((WhitsonGammaModel) thermoSystem.getCharacterization().getPlusFractionModel()) + .setGammaParameters(1.0, 90); + double shape = ((WhitsonGammaModel) thermoSystem.getCharacterization().getPlusFractionModel()) + .getGammaParameters()[0]; + double minMW = ((WhitsonGammaModel) thermoSystem.getCharacterization().getPlusFractionModel()) + .getGammaParameters()[1]; + assertEquals(90.0, minMW, 1e-4); + } }