Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: more robust setting of beta #1253

Merged
merged 3 commits into from
Jan 11, 2025
Merged

feat: more robust setting of beta #1253

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6aa7a78
feat: more robust setting of beta
asmfstatoil Jan 11, 2025
204108a
refact: use phaseFractionMinimumLimit directly
asmfstatoil Jan 11, 2025
deb6b27
ref: reorder math
asmfstatoil Jan 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/main/java/neqsim/thermo/phase/Phase.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand Down
14 changes: 7 additions & 7 deletions src/main/java/neqsim/thermo/system/SystemThermo.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -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);
}
}
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand Down
53 changes: 25 additions & 28 deletions src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand All @@ -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];
}

Expand Down Expand Up @@ -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;
Expand Down
17 changes: 8 additions & 9 deletions src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -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);
}
}

Expand Down
25 changes: 15 additions & 10 deletions src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -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;
}
}

Expand All @@ -294,20 +295,24 @@ 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);
}
}

Qnew = calcQ();

// logger.info("Qold = " + Qold);
// logger.info("Qnew = " + Qnew);

if (Qnew > Qold + 1.0e-10 && !deactivatedPhase) {
// logger.info("Qnew > Qold...............................");
int iter2 = 0;
Expand Down
Loading
Loading