Skip to content

Commit 0fe37d1

Browse files
asmfstatoilasmfstatoil
authored
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
1 parent 0cf54fb commit 0fe37d1

File tree

8 files changed

+96
-86
lines changed

8 files changed

+96
-86
lines changed

‎src/main/java/neqsim/thermo/phase/Phase.java

+3-2
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66

77
package neqsim.thermo.phase;
88

9+
import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
910
import java.util.ArrayList;
1011
import org.apache.logging.log4j.LogManager;
1112
import org.apache.logging.log4j.Logger;
@@ -2020,10 +2021,10 @@ public final double getBeta() {
20202021
@Override
20212022
public final void setBeta(double b) {
20222023
if (b < 0) {
2023-
b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
2024+
b = phaseFractionMinimumLimit;
20242025
}
20252026
if (b > 1) {
2026-
b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
2027+
b = 1.0 - phaseFractionMinimumLimit;
20272028
}
20282029
this.beta = b;
20292030
}

‎src/main/java/neqsim/thermo/system/SystemThermo.java

+7-7
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
package neqsim.thermo.system;
22

3+
import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
34
import java.io.ByteArrayInputStream;
45
import java.io.ByteArrayOutputStream;
56
import java.io.FileInputStream;
@@ -3335,9 +3336,8 @@ public final void initBeta() {
33353336
for (int i = 0; i < numberOfPhases; i++) {
33363337
this.beta[phaseIndex[i]] = getPhase(i).getNumberOfMolesInPhase() / getTotalNumberOfMoles();
33373338
}
3338-
if (!isInitialized
3339-
&& this.getSumBeta() < 1.0 - ThermodynamicModelSettings.phaseFractionMinimumLimit
3340-
|| this.getSumBeta() > 1.0 + ThermodynamicModelSettings.phaseFractionMinimumLimit) {
3339+
if (!isInitialized && this.getSumBeta() < 1.0 - phaseFractionMinimumLimit
3340+
|| this.getSumBeta() > 1.0 + phaseFractionMinimumLimit) {
33413341
logger.warn("SystemThermo:initBeta - Sum of beta does not equal 1.0. " + beta);
33423342
}
33433343
}
@@ -4282,11 +4282,11 @@ public final void setBeta(double b) {
42824282
// TODO: if number of phases > 2, should fail
42834283
if (b < 0) {
42844284
logger.warn("setBeta - Tried to set beta < 0: " + beta);
4285-
b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
4285+
b = phaseFractionMinimumLimit;
42864286
}
42874287
if (b > 1) {
42884288
logger.warn("setBeta - Tried to set beta > 1: " + beta);
4289-
b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
4289+
b = 1.0 - phaseFractionMinimumLimit;
42904290
}
42914291
beta[0] = b;
42924292
beta[1] = 1.0 - b;
@@ -4297,11 +4297,11 @@ public final void setBeta(double b) {
42974297
public final void setBeta(int phaseNum, double b) {
42984298
if (b < 0) {
42994299
logger.warn("setBeta - Tried to set beta < 0: " + beta);
4300-
b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
4300+
b = phaseFractionMinimumLimit;
43014301
}
43024302
if (b > 1) {
43034303
logger.warn("setBeta - Tried to set beta > 1: " + beta);
4304-
b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
4304+
b = 1.0 - phaseFractionMinimumLimit;
43054305
}
43064306
beta[phaseIndex[phaseNum]] = b;
43074307
}

‎src/main/java/neqsim/thermodynamicoperations/flashops/RachfordRice.java

+25-28
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@
33
*
44
* Created on 1.april 2024
55
*/
6-
76
package neqsim.thermodynamicoperations.flashops;
87

8+
import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
99
import java.io.Serializable;
1010
import java.util.Arrays;
1111
import org.apache.logging.log4j.LogManager;
@@ -85,10 +85,9 @@ public double calcBetaMichelsen2001(double[] K, double[] z)
8585
throws neqsim.util.exception.IsNaNException,
8686
neqsim.util.exception.TooManyIterationsException {
8787
int i;
88-
double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
8988
double midler = 0;
90-
double minBeta = tolerance;
91-
double maxBeta = 1.0 - tolerance;
89+
double minBeta = phaseFractionMinimumLimit;
90+
double maxBeta = 1.0 - phaseFractionMinimumLimit;
9291
double g0 = -1.0;
9392
double g1 = 1.0;
9493

@@ -108,10 +107,10 @@ public double calcBetaMichelsen2001(double[] K, double[] z)
108107
// logger.debug("Max beta " + maxBeta + " min beta " + minBeta);
109108

110109
if (g0 < 0) {
111-
return tolerance;
110+
return phaseFractionMinimumLimit;
112111
}
113112
if (g1 > 0) {
114-
return 1.0 - tolerance;
113+
return 1.0 - phaseFractionMinimumLimit;
115114
}
116115

117116
double nybeta = (minBeta + maxBeta) / 2.0;
@@ -191,10 +190,10 @@ public double calcBetaMichelsen2001(double[] K, double[] z)
191190
}
192191
step = gbeta / deriv;
193192
} while (Math.abs(step) >= 1.0e-11 && iterations < maxIterations);
194-
if (nybeta <= tolerance) {
195-
nybeta = tolerance;
196-
} else if (nybeta >= 1.0 - tolerance) {
197-
nybeta = 1.0 - tolerance;
193+
if (nybeta <= phaseFractionMinimumLimit) {
194+
nybeta = phaseFractionMinimumLimit;
195+
} else if (nybeta >= 1.0 - phaseFractionMinimumLimit) {
196+
nybeta = 1.0 - phaseFractionMinimumLimit;
198197
}
199198
beta[0] = nybeta;
200199
beta[1] = 1.0 - nybeta;
@@ -229,7 +228,6 @@ public double calcBetaMichelsen2001(double[] K, double[] z)
229228
public double calcBetaNielsen2023(double[] K, double[] z)
230229
throws neqsim.util.exception.IsNaNException,
231230
neqsim.util.exception.TooManyIterationsException {
232-
double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
233231
double g0 = -1.0;
234232
double g1 = 1.0;
235233

@@ -239,10 +237,10 @@ public double calcBetaNielsen2023(double[] K, double[] z)
239237
}
240238

241239
if (g0 < 0) {
242-
return tolerance;
240+
return phaseFractionMinimumLimit;
243241
}
244242
if (g1 > 0) {
245-
return 1.0 - tolerance;
243+
return 1.0 - phaseFractionMinimumLimit;
246244
}
247245

248246
double V = 0.5;
@@ -329,10 +327,10 @@ public double calcBetaNielsen2023(double[] K, double[] z)
329327
V = 1 - V;
330328
}
331329

332-
if (V <= tolerance) {
333-
V = tolerance;
334-
} else if (V >= 1.0 - tolerance) {
335-
V = 1.0 - tolerance;
330+
if (V <= phaseFractionMinimumLimit) {
331+
V = phaseFractionMinimumLimit;
332+
} else if (V >= 1.0 - phaseFractionMinimumLimit) {
333+
V = 1.0 - phaseFractionMinimumLimit;
336334
}
337335

338336
beta[0] = V;
@@ -380,10 +378,9 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti
380378
ComponentInterface[] compArray = system.getPhase(0).getComponents();
381379

382380
int i;
383-
double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
384381
double midler = 0;
385-
double minBeta = tolerance;
386-
double maxBeta = 1.0 - tolerance;
382+
double minBeta = phaseFractionMinimumLimit;
383+
double maxBeta = 1.0 - phaseFractionMinimumLimit;
387384
double g0 = -1.0;
388385
double g1 = 1.0;
389386

@@ -401,13 +398,13 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti
401398
}
402399

403400
if (g0 < 0) {
404-
this.beta[1] = 1.0 - tolerance;
405-
this.beta[0] = tolerance;
401+
this.beta[1] = 1.0 - phaseFractionMinimumLimit;
402+
this.beta[0] = phaseFractionMinimumLimit;
406403
return this.beta[0];
407404
}
408405
if (g1 > 0) {
409-
this.beta[1] = tolerance;
410-
this.beta[0] = 1.0 - tolerance;
406+
this.beta[1] = phaseFractionMinimumLimit;
407+
this.beta[0] = 1.0 - phaseFractionMinimumLimit;
411408
return this.beta[0];
412409
}
413410

@@ -496,12 +493,12 @@ public final double calcBetaS(SystemInterface system) throws neqsim.util.excepti
496493
step = gbeta / deriv;
497494
} while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // &&
498495

499-
if (nybeta <= tolerance) {
496+
if (nybeta <= phaseFractionMinimumLimit) {
500497
// this.phase = 1;
501-
nybeta = tolerance;
502-
} else if (nybeta >= 1.0 - tolerance) {
498+
nybeta = phaseFractionMinimumLimit;
499+
} else if (nybeta >= 1.0 - phaseFractionMinimumLimit) {
503500
// this.phase = 0;
504-
nybeta = 1.0 - tolerance;
501+
nybeta = 1.0 - phaseFractionMinimumLimit;
505502
// superheated vapour
506503
} else {
507504
// this.phase = 2;

‎src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash.java

+8-9
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
package neqsim.thermodynamicoperations.flashops;
22

3+
import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
34
import org.apache.logging.log4j.LogManager;
45
import org.apache.logging.log4j.Logger;
56
import Jama.Matrix;
@@ -236,16 +237,14 @@ public void solveBeta(boolean ideal) {
236237
}
237238

238239
betaMatrix.minusEquals(ans.times((iter + 1.0) / (10.0 + iter)));
239-
// betaMatrix.print(10, 2);
240-
// betaMatrix.print(10,2);
241-
242240
for (int k = 0; k < system.getNumberOfPhases() - 1; k++) {
243-
system.setBeta(k, betaMatrix.get(k, 0));
244-
if (betaMatrix.get(k, 0) < 0) {
245-
system.setBeta(k, 1e-9);
246-
}
247-
if (betaMatrix.get(k, 0) > 1) {
248-
system.setBeta(k, 1 - 1e-9);
241+
double currBeta = betaMatrix.get(k, 0);
242+
if (currBeta < phaseFractionMinimumLimit) {
243+
system.setBeta(k, phaseFractionMinimumLimit);
244+
} else if (currBeta > 1) {
245+
system.setBeta(k, 1 - phaseFractionMinimumLimit);
246+
} else {
247+
system.setBeta(k, currBeta);
249248
}
250249
}
251250

‎src/main/java/neqsim/thermodynamicoperations/flashops/SolidFlash1.java

+15-10
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
package neqsim.thermodynamicoperations.flashops;
22

3+
import static neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
34
import org.apache.logging.log4j.LogManager;
45
import org.apache.logging.log4j.Logger;
56
import Jama.Matrix;
@@ -272,10 +273,10 @@ public void solveBeta() {
272273
double minBetaTem = 1000000;
273274
int minBetaTemIndex = 0;
274275

275-
for (int i = 0; i < system.getNumberOfPhases() - solidsNumber; i++) {
276-
if (betaMatrixTemp.get(i, 0) * FluidPhaseActiveDescriptors[i] < minBetaTem) {
277-
minBetaTem = betaMatrixTemp.get(i, 0);
278-
minBetaTemIndex = i;
276+
for (int k = 0; k < system.getNumberOfPhases() - solidsNumber; k++) {
277+
if (betaMatrixTemp.get(k, 0) * FluidPhaseActiveDescriptors[k] < minBetaTem) {
278+
minBetaTem = betaMatrixTemp.get(k, 0);
279+
minBetaTemIndex = k;
279280
}
280281
}
281282

@@ -294,20 +295,24 @@ public void solveBeta() {
294295
betaMatrixOld = betaMatrix.copy();
295296
betaMatrix = betaMatrixTemp.copy();
296297
boolean deactivatedPhase = false;
297-
for (int i = 0; i < system.getNumberOfPhases() - solidsNumber; i++) {
298-
system.setBeta(i, betaMatrix.get(i, 0));
299-
// logger.info("Fluid Phase fraction" + system.getBeta(i));
300-
if (Math.abs(system.getBeta(i)) < 1.0e-10) {
301-
FluidPhaseActiveDescriptors[i] = 0;
298+
for (int k = 0; k < system.getNumberOfPhases() - solidsNumber; k++) {
299+
double currBeta = betaMatrix.get(k, 0);
300+
if (currBeta < phaseFractionMinimumLimit) {
301+
system.setBeta(k, phaseFractionMinimumLimit);
302+
// This used to be verified with abs(betamatrix.get(i,0)) < 1.0e-10
303+
FluidPhaseActiveDescriptors[k] = 0;
302304
deactivatedPhase = true;
305+
} else if (currBeta > (1.0 - phaseFractionMinimumLimit)) {
306+
system.setBeta(k, 1.0 - phaseFractionMinimumLimit);
307+
} else {
308+
system.setBeta(k, currBeta);
303309
}
304310
}
305311

306312
Qnew = calcQ();
307313

308314
// logger.info("Qold = " + Qold);
309315
// logger.info("Qnew = " + Qnew);
310-
311316
if (Qnew > Qold + 1.0e-10 && !deactivatedPhase) {
312317
// logger.info("Qnew > Qold...............................");
313318
int iter2 = 0;

0 commit comments

Comments
 (0)