From 0fe37d1b13f5d2daa6298600f545db8c8e52667b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=85smund=20V=C3=A5ge=20Fannemel?= <34712686+asmfstatoil@users.noreply.github.com> Date: Sat, 11 Jan 2025 18:35:15 +0100 Subject: [PATCH] feat: more robust setting of beta (#1253) * feat: more robust setting of beta using neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit * refact: use phaseFractionMinimumLimit directly * ref: reorder math --- src/main/java/neqsim/thermo/phase/Phase.java | 5 +- .../neqsim/thermo/system/SystemThermo.java | 14 ++--- .../flashops/RachfordRice.java | 53 +++++++++--------- .../flashops/SolidFlash.java | 17 +++--- .../flashops/SolidFlash1.java | 25 +++++---- .../flashops/TPflash.java | 54 +++++++++++-------- .../flashops/TPmultiflash.java | 11 ++-- .../flashops/dTPflash.java | 3 -- 8 files changed, 96 insertions(+), 86 deletions(-) diff --git a/src/main/java/neqsim/thermo/phase/Phase.java b/src/main/java/neqsim/thermo/phase/Phase.java index 44b69fbbd..c106d3591 100644 --- a/src/main/java/neqsim/thermo/phase/Phase.java +++ b/src/main/java/neqsim/thermo/phase/Phase.java @@ -6,6 +6,7 @@ package neqsim.thermo.phase; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import java.util.ArrayList; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; @@ -2020,10 +2021,10 @@ public final double getBeta() { @Override public final void setBeta(double b) { if (b < 0) { - b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = phaseFractionMinimumLimit; } if (b > 1) { - b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = 1.0 - phaseFractionMinimumLimit; } this.beta = b; } diff --git a/src/main/java/neqsim/thermo/system/SystemThermo.java b/src/main/java/neqsim/thermo/system/SystemThermo.java index bdf730cf0..0333a915c 100644 --- a/src/main/java/neqsim/thermo/system/SystemThermo.java +++ b/src/main/java/neqsim/thermo/system/SystemThermo.java @@ -1,5 +1,6 @@ package neqsim.thermo.system; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import java.io.ByteArrayInputStream; import java.io.ByteArrayOutputStream; import java.io.FileInputStream; @@ -3335,9 +3336,8 @@ public final void initBeta() { for (int i = 0; i < numberOfPhases; i++) { this.beta[phaseIndex[i]] = getPhase(i).getNumberOfMolesInPhase() / getTotalNumberOfMoles(); } - if (!isInitialized - && this.getSumBeta() < 1.0 - ThermodynamicModelSettings.phaseFractionMinimumLimit - || this.getSumBeta() > 1.0 + ThermodynamicModelSettings.phaseFractionMinimumLimit) { + if (!isInitialized && this.getSumBeta() < 1.0 - phaseFractionMinimumLimit + || this.getSumBeta() > 1.0 + phaseFractionMinimumLimit) { logger.warn("SystemThermo:initBeta - Sum of beta does not equal 1.0. " + beta); } } @@ -4282,11 +4282,11 @@ public final void setBeta(double b) { // TODO: if number of phases > 2, should fail if (b < 0) { logger.warn("setBeta - Tried to set beta < 0: " + beta); - b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = phaseFractionMinimumLimit; } if (b > 1) { logger.warn("setBeta - Tried to set beta > 1: " + beta); - b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = 1.0 - phaseFractionMinimumLimit; } beta[0] = b; beta[1] = 1.0 - b; @@ -4297,11 +4297,11 @@ public final void setBeta(double b) { public final void setBeta(int phaseNum, double b) { if (b < 0) { logger.warn("setBeta - Tried to set beta < 0: " + beta); - b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = phaseFractionMinimumLimit; } if (b > 1) { logger.warn("setBeta - Tried to set beta > 1: " + beta); - b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; + b = 1.0 - phaseFractionMinimumLimit; } beta[phaseIndex[phaseNum]] = b; } diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java b/src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java index 65117fd64..cf22889ff 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java @@ -3,9 +3,9 @@ * * Created on 1.april 2024 */ - package neqsim.thermodynamicoperations.flashops; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import java.io.Serializable; import java.util.Arrays; import org.apache.logging.log4j.LogManager; @@ -85,10 +85,9 @@ public double calcBetaMichelsen2001(double[] K, double[] z) throws neqsim.util.exception.IsNaNException, neqsim.util.exception.TooManyIterationsException { int i; - double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; double midler = 0; - double minBeta = tolerance; - double maxBeta = 1.0 - tolerance; + double minBeta = phaseFractionMinimumLimit; + double maxBeta = 1.0 - phaseFractionMinimumLimit; double g0 = -1.0; double g1 = 1.0; @@ -108,10 +107,10 @@ public double calcBetaMichelsen2001(double[] K, double[] z) // logger.debug("Max beta " + maxBeta + " min beta " + minBeta); if (g0 < 0) { - return tolerance; + return phaseFractionMinimumLimit; } if (g1 > 0) { - return 1.0 - tolerance; + return 1.0 - phaseFractionMinimumLimit; } double nybeta = (minBeta + maxBeta) / 2.0; @@ -191,10 +190,10 @@ public double calcBetaMichelsen2001(double[] K, double[] z) } step = gbeta / deriv; } while (Math.abs(step) >= 1.0e-11 && iterations < maxIterations); - if (nybeta <= tolerance) { - nybeta = tolerance; - } else if (nybeta >= 1.0 - tolerance) { - nybeta = 1.0 - tolerance; + if (nybeta <= phaseFractionMinimumLimit) { + nybeta = phaseFractionMinimumLimit; + } else if (nybeta >= 1.0 - phaseFractionMinimumLimit) { + nybeta = 1.0 - phaseFractionMinimumLimit; } beta[0] = nybeta; beta[1] = 1.0 - nybeta; @@ -229,7 +228,6 @@ public double calcBetaMichelsen2001(double[] K, double[] z) public double calcBetaNielsen2023(double[] K, double[] z) throws neqsim.util.exception.IsNaNException, neqsim.util.exception.TooManyIterationsException { - double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; double g0 = -1.0; double g1 = 1.0; @@ -239,10 +237,10 @@ public double calcBetaNielsen2023(double[] K, double[] z) } if (g0 < 0) { - return tolerance; + return phaseFractionMinimumLimit; } if (g1 > 0) { - return 1.0 - tolerance; + return 1.0 - phaseFractionMinimumLimit; } double V = 0.5; @@ -329,10 +327,10 @@ public double calcBetaNielsen2023(double[] K, double[] z) V = 1 - V; } - if (V <= tolerance) { - V = tolerance; - } else if (V >= 1.0 - tolerance) { - V = 1.0 - tolerance; + if (V <= phaseFractionMinimumLimit) { + V = phaseFractionMinimumLimit; + } else if (V >= 1.0 - phaseFractionMinimumLimit) { + V = 1.0 - phaseFractionMinimumLimit; } beta[0] = V; @@ -380,10 +378,9 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti ComponentInterface[] compArray = system.getPhase(0).getComponents(); int i; - double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; double midler = 0; - double minBeta = tolerance; - double maxBeta = 1.0 - tolerance; + double minBeta = phaseFractionMinimumLimit; + double maxBeta = 1.0 - phaseFractionMinimumLimit; double g0 = -1.0; double g1 = 1.0; @@ -401,13 +398,13 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti } if (g0 < 0) { - this.beta[1] = 1.0 - tolerance; - this.beta[0] = tolerance; + this.beta[1] = 1.0 - phaseFractionMinimumLimit; + this.beta[0] = phaseFractionMinimumLimit; return this.beta[0]; } if (g1 > 0) { - this.beta[1] = tolerance; - this.beta[0] = 1.0 - tolerance; + this.beta[1] = phaseFractionMinimumLimit; + this.beta[0] = 1.0 - phaseFractionMinimumLimit; return this.beta[0]; } @@ -496,12 +493,12 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti step = gbeta / deriv; } while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // && - if (nybeta <= tolerance) { + if (nybeta <= phaseFractionMinimumLimit) { // this.phase = 1; - nybeta = tolerance; - } else if (nybeta >= 1.0 - tolerance) { + nybeta = phaseFractionMinimumLimit; + } else if (nybeta >= 1.0 - phaseFractionMinimumLimit) { // this.phase = 0; - nybeta = 1.0 - tolerance; + nybeta = 1.0 - phaseFractionMinimumLimit; // superheated vapour } else { // this.phase = 2; diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java index 2848fb179..1f8a57c64 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java @@ -1,5 +1,6 @@ package neqsim.thermodynamicoperations.flashops; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import Jama.Matrix; @@ -236,16 +237,14 @@ public void solveBeta(boolean ideal) { } betaMatrix.minusEquals(ans.times((iter + 1.0) / (10.0 + iter))); - // betaMatrix.print(10, 2); - // betaMatrix.print(10,2); - for (int k = 0; k < system.getNumberOfPhases() - 1; k++) { - system.setBeta(k, betaMatrix.get(k, 0)); - if (betaMatrix.get(k, 0) < 0) { - system.setBeta(k, 1e-9); - } - if (betaMatrix.get(k, 0) > 1) { - system.setBeta(k, 1 - 1e-9); + double currBeta = betaMatrix.get(k, 0); + if (currBeta < phaseFractionMinimumLimit) { + system.setBeta(k, phaseFractionMinimumLimit); + } else if (currBeta > 1) { + system.setBeta(k, 1 - phaseFractionMinimumLimit); + } else { + system.setBeta(k, currBeta); } } diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java b/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java index a2f1aead0..47f7c902e 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java @@ -1,5 +1,6 @@ package neqsim.thermodynamicoperations.flashops; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import Jama.Matrix; @@ -272,10 +273,10 @@ public void solveBeta() { double minBetaTem = 1000000; int minBetaTemIndex = 0; - for (int i = 0; i < system.getNumberOfPhases() - solidsNumber; i++) { - if (betaMatrixTemp.get(i, 0) * FluidPhaseActiveDescriptors[i] < minBetaTem) { - minBetaTem = betaMatrixTemp.get(i, 0); - minBetaTemIndex = i; + for (int k = 0; k < system.getNumberOfPhases() - solidsNumber; k++) { + if (betaMatrixTemp.get(k, 0) * FluidPhaseActiveDescriptors[k] < minBetaTem) { + minBetaTem = betaMatrixTemp.get(k, 0); + minBetaTemIndex = k; } } @@ -294,12 +295,17 @@ public void solveBeta() { betaMatrixOld = betaMatrix.copy(); betaMatrix = betaMatrixTemp.copy(); boolean deactivatedPhase = false; - for (int i = 0; i < system.getNumberOfPhases() - solidsNumber; i++) { - system.setBeta(i, betaMatrix.get(i, 0)); - // logger.info("Fluid Phase fraction" + system.getBeta(i)); - if (Math.abs(system.getBeta(i)) < 1.0e-10) { - FluidPhaseActiveDescriptors[i] = 0; + for (int k = 0; k < system.getNumberOfPhases() - solidsNumber; k++) { + double currBeta = betaMatrix.get(k, 0); + if (currBeta < phaseFractionMinimumLimit) { + system.setBeta(k, phaseFractionMinimumLimit); + // This used to be verified with abs(betamatrix.get(i,0)) < 1.0e-10 + FluidPhaseActiveDescriptors[k] = 0; deactivatedPhase = true; + } else if (currBeta > (1.0 - phaseFractionMinimumLimit)) { + system.setBeta(k, 1.0 - phaseFractionMinimumLimit); + } else { + system.setBeta(k, currBeta); } } @@ -307,7 +313,6 @@ public void solveBeta() { // logger.info("Qold = " + Qold); // logger.info("Qnew = " + Qnew); - if (Qnew > Qold + 1.0e-10 && !deactivatedPhase) { // logger.info("Qnew > Qold..............................."); int iter2 = 0; diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java index a8e4d6955..98930cb7a 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java @@ -1,5 +1,6 @@ package neqsim.thermodynamicoperations.flashops; +import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import neqsim.thermo.phase.PhaseType; @@ -23,7 +24,6 @@ public class TPflash extends Flash { static Logger logger = LogManager.getLogger(TPflash.class); SystemInterface clonedSystem; - double betaTolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit; double presdiff = 1.0; /** @@ -66,7 +66,7 @@ public TPflash(SystemInterface system, boolean checkForSolids) { /** *

- * sucsSubs. + * sucsSubs. Successive substitutions. *

*/ public void sucsSubs() { @@ -104,11 +104,11 @@ public void sucsSubs() { logger.warn("Not able to calculate beta, calculation is not converging."); system.setBeta(oldBeta); } - if (system.getBeta() > 1.0 - betaTolerance) { - system.setBeta(1.0 - betaTolerance); + if (system.getBeta() > 1.0 - phaseFractionMinimumLimit) { + system.setBeta(1.0 - phaseFractionMinimumLimit); } - if (system.getBeta() < betaTolerance) { - system.setBeta(betaTolerance); + if (system.getBeta() < phaseFractionMinimumLimit) { + system.setBeta(phaseFractionMinimumLimit); } system.calc_x_y(); system.init(1); @@ -142,7 +142,8 @@ public void accselerateSucsSubs() { system.setBeta(rachfordRice.calcBeta(system.getKvector(), system.getzvector())); } catch (Exception ex) { system.setBeta(rachfordRice.getBeta()[0]); - if (system.getBeta() > 1.0 - betaTolerance || system.getBeta() < betaTolerance) { + if (system.getBeta() > 1.0 - phaseFractionMinimumLimit + || system.getBeta() < phaseFractionMinimumLimit) { system.setBeta(oldBeta); } // logger.info("temperature " + system.getTemperature() + " pressure " + @@ -195,7 +196,20 @@ public void resetK() { } } - /** {@inheritDoc} */ + /** + * {@inheritDoc} + * + *

+ * Calculate the following properties: + *

+ * + */ @Override public void run() { if (system.isForcePhaseTypes() && system.getMaxNumberOfPhases() == 1) { @@ -273,16 +287,16 @@ public void run() { // This solves some problems when we have high volumes of water and heavy // hydrocarbons returning only one liquid phase (and this phase desolves all // gas) - if (system.getBeta() > (1.0 - betaTolerance * 1.1) - || system.getBeta() < (betaTolerance * 1.1)) { + if (system.getBeta() > (1.0 - 1.1 * phaseFractionMinimumLimit) + || system.getBeta() < (1.1 * phaseFractionMinimumLimit)) { system.setBeta(0.5); sucsSubs(); } // Performs three iterations of successive substitution for (int k = 0; k < 3; k++) { - if (system.getBeta() < (1.0 - betaTolerance * 1.1) - && system.getBeta() > (betaTolerance * 1.1)) { + if (system.getBeta() < (1.0 - 1.1 * phaseFractionMinimumLimit) + && system.getBeta() > (1.1 * phaseFractionMinimumLimit)) { sucsSubs(); if ((system.getGibbsEnergy() - minimumGibbsEnergy) / Math.abs(minimumGibbsEnergy) < -1e-12) { @@ -292,14 +306,13 @@ public void run() { } // System.out.println("beta " + system.getBeta()); - int totiter = 0; double tpdx = 1.0; double tpdy = 1.0; double dgonRT = 1.0; boolean passedTests = false; - if (system.getBeta() > (1.0 - 1.1 * betaTolerance) - || system.getBeta() < (1.1 * betaTolerance)) { + if (system.getBeta() > (1.0 - 1.1 * phaseFractionMinimumLimit) + || system.getBeta() < (1.1 * phaseFractionMinimumLimit)) { tpdx = 1.0; tpdy = 1.0; dgonRT = 1.0; @@ -322,7 +335,6 @@ public void run() { } dgonRT = system.getPhase(0).getBeta() * tpdy + (1.0 - system.getPhase(0).getBeta()) * tpdx; - if (dgonRT > 0) { if (tpdx < 0) { for (i = 0; i < system.getPhases()[0].getNumberOfComponents(); i++) { @@ -363,8 +375,6 @@ public void run() { system.orderByDensity(); system.init(1); - // commented out by Even Solbraa 6/2-2012k - // system.init(3); return; } } @@ -424,9 +434,9 @@ public void run() { gibbsEnergy = system.getGibbsEnergy(); if (((gibbsEnergy - gibbsEnergyOld) / Math.abs(gibbsEnergyOld) > 1e-8 - || system.getBeta() < betaTolerance * 1.01 - || system.getBeta() > (1 - betaTolerance * 1.01)) && !system.isChemicalSystem() - && timeFromLastGibbsFail > 0) { + || system.getBeta() < phaseFractionMinimumLimit * 1.01 + || system.getBeta() > (1 - phaseFractionMinimumLimit * 1.01)) + && !system.isChemicalSystem() && timeFromLastGibbsFail > 0) { resetK(); timeFromLastGibbsFail = 0; // logger.info("gibbs decrease " + (gibbsEnergy - gibbsEnergyOld) / @@ -505,7 +515,7 @@ public void run() { } for (int i = 0; i < system.getNumberOfPhases(); i++) { - if (system.getBeta(i) < betaTolerance * 1.01) { + if (system.getBeta(i) < phaseFractionMinimumLimit * 1.01) { system.removePhase(i); } } diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPmultiflash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPmultiflash.java index f01b5571c..023ad86c2 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPmultiflash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPmultiflash.java @@ -11,7 +11,6 @@ import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import org.ejml.simple.SimpleMatrix; -import neqsim.thermo.ThermodynamicModelSettings; import neqsim.thermo.phase.PhaseType; import neqsim.thermo.system.SystemInterface; @@ -207,16 +206,18 @@ public double solveBeta() { betaMatrix = betaMatrix.minus(ans.scale(iter / (iter + 3.0))); removePhase = false; for (int k = 0; k < system.getNumberOfPhases(); k++) { - system.setBeta(k, betaMatrix.get(0, k)); - if (betaMatrix.get(0, k) < phaseFractionMinimumLimit) { + double currBeta = betaMatrix.get(0, k); + if (currBeta < phaseFractionMinimumLimit) { system.setBeta(k, phaseFractionMinimumLimit); if (checkOneRemove) { checkOneRemove = false; removePhase = true; } checkOneRemove = true; - } else if (betaMatrix.get(0, k) > (1.0 - phaseFractionMinimumLimit)) { + } else if (currBeta > (1.0 - phaseFractionMinimumLimit)) { system.setBeta(k, 1.0 - phaseFractionMinimumLimit); + } else { + system.setBeta(k, currBeta); } } system.normalizeBeta(); @@ -1347,7 +1348,7 @@ public void run() { boolean hasRemovedPhase = false; for (int i = 0; i < system.getNumberOfPhases(); i++) { - if (system.getBeta(i) < ThermodynamicModelSettings.phaseFractionMinimumLimit * 1.1) { + if (system.getBeta(i) < 1.1 * phaseFractionMinimumLimit) { system.removePhaseKeepTotalComposition(i); doStabilityAnalysis = false; hasRemovedPhase = true; diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/dTPflash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/dTPflash.java index f8405db90..5fcd92f15 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/dTPflash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/dTPflash.java @@ -46,7 +46,6 @@ public void run() { double diff = 0.0; // double fracdiff = 0.0; - // system.setBeta(0.5); do { diff = 0.0; // fracdiff = 0.0; @@ -81,8 +80,6 @@ public void run() { // if(!hasgot) system.getPhase(1).getComponent(i).setx(1e-16); } - // system.setBeta(0.5+fracdiff); - system.getPhase(1).normalize(); logger.info("diff " + diff); } while (diff > 1e-10 && iterations < 1000);