Skip to content

Commit

Permalink
feat: more robust setting of beta (#1253)
Browse files Browse the repository at this point in the history 
* feat: more robust setting of beta

using neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit

* refact: use phaseFractionMinimumLimit directly

* ref: reorder math
  • Loading branch information
@asmfstatoil
asmfstatoil authored Jan 11, 2025
1 parent 0cf54fb commit 0fe37d1
Show file tree
Hide file tree
Showing 8 changed files with 96 additions and 86 deletions.
5 changes: 3 additions & 2 deletions src/main/java/neqsim/thermo/phase/Phase.java
View file
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
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
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
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
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

0 comments on commit 0fe37d1

Please sign in to comment.