Skip to content

Commit

Permalink
improve flash
Browse files Browse the repository at this point in the history
  • Loading branch information
@EvenSol
EvenSol committed Nov 26, 2023
1 parent 3910d39 commit af10b76
Show file tree
Hide file tree
Showing 2 changed files with 351 additions and 20 deletions.
348 changes: 347 additions & 1 deletion src/main/java/neqsim/thermodynamicOperations/flashOps/TPmultiflash.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -573,6 +573,351 @@ public void stabilityAnalysis() {
// system.display();
}

public void stabilityAnalysis3() {
double[] logWi = new double[system.getPhase(0).getNumberOfComponents()];
double[][] Wi = new double[system.getPhase(0).getNumberOfComponents()][system.getPhase(0)
.getNumberOfComponents()];

double[] deltalogWi = new double[system.getPhases()[0].getNumberOfComponents()];
double[] oldDeltalogWi = new double[system.getPhases()[0].getNumberOfComponents()];
double[] oldoldDeltalogWi = new double[system.getPhases()[0].getNumberOfComponents()];
double[] sumw = new double[system.getPhase(0).getNumberOfComponents()];
double err = 0;
double[] oldlogw = new double[system.getPhase(0).getNumberOfComponents()];
double[] oldoldlogw = new double[system.getPhases()[0].getNumberOfComponents()];
double[] oldoldoldlogw = new double[system.getPhases()[0].getNumberOfComponents()];
double[] d = new double[system.getPhase(0).getNumberOfComponents()];
double[][] x = new double[system.getPhase(0).getNumberOfComponents()][system.getPhase(0)
.getNumberOfComponents()];
tm = new double[system.getPhase(0).getNumberOfComponents()];

double[] alpha = null;
// SystemInterface minimumGibbsEnergySystem;
ArrayList<SystemInterface> clonedSystem = new ArrayList<SystemInterface>(1);
// if (minimumGibbsEnergySystem == null) {
// minimumGibbsEnergySystem = system.clone();
// }
minimumGibbsEnergySystem = system;
clonedSystem.add(system.clone());
/*
* for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { if
* (system.getPhase(0).getComponent(i).getx() < 1e-100) { clonedSystem.add(null); continue; }
* double numb = 0; clonedSystem.add(system.clone());
*
* // (clonedSystem.get(i)).init(0); commented out sept 2005, Even S. for (int j = 0; j <
* system.getPhase(0).getNumberOfComponents(); j++) { numb = i == j ? 1.0 : 1.0e-12; // set to 0
* by Even Solbraa 23.01.2013 - chaged back to 1.0e-12 27.04.13 if
* (system.getPhase(0).getComponent(j).getz() < 1e-100) { numb = 0; } (
* clonedSystem.get(i)).getPhase(1).getComponents()[j].setx(numb); } if
* (system.getPhase(0).getComponent(i).getIonicCharge() == 0) { ( clonedSystem.get(i)).init(1);
* } }
*/

lowestGibbsEnergyPhase = 0;
/*
* // logger.info("low gibbs phase " + lowestGibbsEnergyPhase); for (int k = 0; k <
* minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); k++) { for (int i = 0; i <
* minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); i++) { if (!((
* clonedSystem.get(k)) == null)) { sumw[k] += (
* clonedSystem.get(k)).getPhase(1).getComponents()[i].getx(); } } }
*
* for (int k = 0; k < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); k++) { for
* (int i = 0; i < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); i++) { if (!((
* clonedSystem.get(k)) == null) && system.getPhase(0).getComponent(k).getx() > 1e-100) { (
* clonedSystem.get(k)).getPhase(1).getComponents()[i].setx((
* clonedSystem.get(k)).getPhase(1).getComponents()[i].getx() / sumw[0]); } logger.info("x: " +
* ( clonedSystem.get(k)).getPhase(0).getComponents()[i].getx()); } if
* (system.getPhase(0).getComponent(k).getx() > 1e-100) { d[k] =
* Math.log(system.getPhase(0).getComponents()[k].getx()) +
* system.getPhase(0).getComponents()[k].getLogFugacityCoefficient();
* if(minimumGibbsEnergySystem.getPhases()[lowestGibbsEnergyPhase].getComponents
* ()[k].getIonicCharge()!=0) d[k]=0; } //logger.info("dk: " + d[k]); }
*/
for (int k = 0; k < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); k++) {
if (system.getPhase(0).getComponent(k).getx() > 1e-100) {
d[k] = Math.log(system.getPhase(0).getComponents()[k].getx())
+ system.getPhase(0).getComponents()[k].getLogFugacityCoefficient();
// if(minimumGibbsEnergySystem.getPhases()[lowestGibbsEnergyPhase].getComponents()[k].getIonicCharge()!=0)
// d[k]=0;
}
}

for (int j = 0; j < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); j++) {
if (system.getPhase(0).getComponent(j).getz() > 1e-100) {
logWi[j] = 1.0;
} else {
logWi[j] = -10000.0;
}
}

int hydrocarbonTestCompNumb = 0;
int lightTestCompNumb = 0;
double Mmax = 0;
double Mmin = 1e10;
for (int i = 0; i < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); i++) {
if (minimumGibbsEnergySystem.getPhase(0).getComponent(i).isHydrocarbon()) {
if ((minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass()) > Mmax) {
Mmax = minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass();
}
if ((minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass()) < Mmin) {
Mmin = minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass();
}
}
}

for (int i = 0; i < minimumGibbsEnergySystem.getPhase(0).getNumberOfComponents(); i++) {
if (minimumGibbsEnergySystem.getPhase(0).getComponent(i).isHydrocarbon()
&& minimumGibbsEnergySystem.getPhase(0).getComponent(i).getz() > 1e-50) {
if (Math.abs(
(minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass()) - Mmax) < 1e-5) {
hydrocarbonTestCompNumb = i;
// logger.info("CHECKING heavy component " + hydrocarbonTestCompNumb);
}
}

if (minimumGibbsEnergySystem.getPhase(0).getComponent(i).isHydrocarbon()
&& minimumGibbsEnergySystem.getPhase(0).getComponent(i).getz() > 1e-50) {
if (Math.abs(
(minimumGibbsEnergySystem.getPhase(0).getComponent(i).getMolarMass()) - Mmin) < 1e-5) {
lightTestCompNumb = i;
// logger.info("CHECKING light component " + lightTestCompNumb);
}
}
}

for (int j = 0; j < system.getNumberOfComponents(); j++) {
if (minimumGibbsEnergySystem.getPhase(0).getComponent(j).getx() < 1e-100
|| (minimumGibbsEnergySystem.getPhase(0).getComponent(j).getIonicCharge() != 0)
|| (minimumGibbsEnergySystem.getPhase(0).getComponent(j).isHydrocarbon()
&& j != hydrocarbonTestCompNumb && j != lightTestCompNumb)) {
continue;
}

double nomb = 0.0;
for (int cc = 0; cc < system.getPhase(0).getNumberOfComponents(); cc++) {
nomb = cc == j ? 1.0 : 1.0e-12;
if (system.getPhase(0).getComponent(cc).getz() < 1e-100) {
nomb = 0.0;
}

if (clonedSystem.get(0).isPhase(1)) {
try {
clonedSystem.get(0).getPhase(1).getComponents()[cc].setx(nomb);
} catch (Exception ex) {
logger.warn(ex.getMessage());
}
}
}
// if(minimumGibbsEnergySystem.getPhase(0).getComponent(j).getName().equals("water")
// && minimumGibbsEnergySystem.isChemicalSystem()) continue;
// logger.info("STAB CHECK COMP " +
// system.getPhase(0).getComponent(j).getComponentName());
// if(minimumGibbsEnergySystem.getPhase(0).getComponent(j).isInert()) break;
int iter = 0;
double errOld = 1.0e100;
boolean useaccsubst = true;
do {
errOld = err;
iter++;
err = 0;

if (iter <= 150 || !system.isImplementedCompositionDeriativesofFugacity()) {
if (iter % 7 == 0 && useaccsubst) {
double vec1 = 0.0;

double vec2 = 0.0;
double prod1 = 0.0;
double prod2 = 0.0;
for (i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
vec1 = oldDeltalogWi[i] * oldoldDeltalogWi[i];
vec2 = Math.pow(oldoldDeltalogWi[i], 2.0);
prod1 += vec1 * vec2;
prod2 += vec2 * vec2;
}

double lambda = prod1 / prod2;
// logger.info("lambda " + lambda);
for (i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
logWi[i] += lambda / (1.0 - lambda) * deltalogWi[i];
err += Math.abs((logWi[i] - oldlogw[i]) / oldlogw[i]);
Wi[j][i] = Math.exp(logWi[i]);
}
} else {
for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
oldoldoldlogw[i] = oldoldlogw[i];
oldoldlogw[i] = oldlogw[i];
oldlogw[i] = logWi[i];
oldoldDeltalogWi[i] = oldoldlogw[i] - oldoldoldlogw[i];
oldDeltalogWi[i] = oldlogw[i] - oldoldlogw[i];
}
clonedSystem.get(0).init(1, 1);
for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
// oldlogw[i] = logWi[i];
if (!Double.isInfinite(
clonedSystem.get(0).getPhase(1).getComponents()[i].getLogFugacityCoefficient())
&& system.getPhase(0).getComponent(i).getx() > 1e-100) {
logWi[i] = d[i] - clonedSystem.get(0).getPhase(1).getComponents()[i]
.getLogFugacityCoefficient();
if (clonedSystem.get(0).getPhase(1).getComponents()[i].getIonicCharge() != 0) {
logWi[i] = -1000.0;
}
}
deltalogWi[i] = logWi[i] - oldlogw[i];
err += Math.abs(logWi[i] - oldlogw[i]);
Wi[j][i] = Math.exp(logWi[i]);
useaccsubst = true;
}
if (iter > 2 && err > errOld) {
useaccsubst = false;
}
}

} else {
SimpleMatrix f = new SimpleMatrix(system.getPhases()[0].getNumberOfComponents(), 1);
SimpleMatrix df = null;
SimpleMatrix identitytimesConst = null;
// if (!secondOrderStabilityAnalysis) {
for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
oldoldoldlogw[i] = oldoldlogw[i];
oldoldlogw[i] = oldlogw[i];
oldlogw[i] = logWi[i];
oldoldDeltalogWi[i] = oldoldlogw[i] - oldoldoldlogw[i];
oldDeltalogWi[i] = oldlogw[i] - oldoldlogw[i];
}
clonedSystem.get(0).init(3, 1);
alpha = new double[clonedSystem.get(0).getPhases()[0].getNumberOfComponents()];
df = new SimpleMatrix(system.getPhases()[0].getNumberOfComponents(),
system.getPhases()[0].getNumberOfComponents());
identitytimesConst = SimpleMatrix.identity(system.getPhases()[0].getNumberOfComponents());
// ,
// system.getPhases()[0].getNumberOfComponents());
// secondOrderStabilityAnalysis = true;
// }

for (int i = 0; i < clonedSystem.get(0).getPhases()[0].getNumberOfComponents(); i++) {
alpha[i] = 2.0 * Math.sqrt(Wi[j][i]);
}

for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
if (system.getPhase(0).getComponent(i).getz() > 1e-100) {
f.set(i, 0,
Math.sqrt(Wi[j][i])
* (Math.log(Wi[j][i]) + clonedSystem.get(0).getPhases()[1].getComponents()[i]
.getLogFugacityCoefficient() - d[i]));
}
for (int k = 0; k < clonedSystem.get(0).getPhases()[0].getNumberOfComponents(); k++) {
double kronDelt = (i == k) ? 1.0 : 0.0;
if (system.getPhase(0).getComponent(i).getz() > 1e-100) {
df.set(i, k, kronDelt + Math.sqrt(Wi[j][k] * Wi[j][i])
* clonedSystem.get(0).getPhases()[1].getComponents()[i].getdfugdn(k));
// * clonedSystem.getPhases()[j].getNumberOfMolesInPhase());
} else {
df.set(i, k, 0);
// * clonedSystem.getPhases()[j].getNumberOfMolesInPhase());
}
}
}

// f.print(10, 10);
// df.print(10, 10);
SimpleMatrix dx = null;
try {
dx = df.plus(identitytimesConst).solve(f).negative();
} catch (Exception e) {
dx = df.plus(identitytimesConst.scale(0.5)).solve(f).negative();
}

// dx.print(10, 10);

for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
double alphaNew = alpha[i] + dx.get(i, 0);
Wi[j][i] = Math.pow(alphaNew / 2.0, 2.0);
if (system.getPhase(0).getComponent(i).getz() > 1e-100) {
logWi[i] = Math.log(Wi[j][i]);
}
if (system.getPhase(0).getComponent(i).getIonicCharge() != 0) {
logWi[i] = -1000.0;
}
err += Math.abs((logWi[i] - oldlogw[i]) / oldlogw[i]);
}

// logger.info("err newton " + err);
}
// logger.info("err: " + err);
sumw[j] = 0;

for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
sumw[j] += Math.exp(logWi[i]);
}

for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
if (system.getPhase(0).getComponent(i).getx() > 1e-100) {
clonedSystem.get(0).getPhase(1).getComponents()[i].setx(Math.exp(logWi[i]) / sumw[j]);
}
if (system.getPhase(0).getComponent(i).getIonicCharge() != 0) {
clonedSystem.get(0).getPhase(1).getComponents()[i].setx(1e-50);
}
}
} while ((Math.abs(err) > 1e-9 || err > errOld) && iter < 600);
// logger.info("err: " + err + " ITER " + iter);
double xTrivialCheck0 = 0.0;
double xTrivialCheck1 = 0.0;

tm[j] = 1.0;

for (int i = 0; i < system.getPhase(1).getNumberOfComponents(); i++) {
if (system.getPhase(0).getComponent(i).getx() > 1e-100) {
tm[j] -= Math.exp(logWi[i]);
}
x[j][i] = clonedSystem.get(0).getPhase(1).getComponents()[i].getx();
// logger.info("txji: " + x[j][i]);

xTrivialCheck0 += Math.abs(x[j][i] - system.getPhase(0).getComponent(i).getx());
xTrivialCheck1 += Math.abs(x[j][i] - system.getPhase(1).getComponent(i).getx());
}
if (iter >= 599) {
// logger.info("iter > maxiter multiphase stability ");
// logger.info("error " + Math.abs(err));
// logger.info("tm: " + tm[j]);
}

if (Math.abs(xTrivialCheck0) < 1e-4 || Math.abs(xTrivialCheck1) < 1e-4) {
tm[j] = 10.0;
}

if (tm[j] < -1e-8) {
break;
}
}

int unstabcomp = 0;
for (int k = system.getPhase(0).getNumberOfComponents() - 1; k >= 0; k--) {
if (tm[k] < -1e-8 && !(Double.isNaN(tm[k]))) {
system.addPhase();
unstabcomp = k;
for (int i = 0; i < system.getPhase(1).getNumberOfComponents(); i++) {
system.getPhase(system.getNumberOfPhases() - 1).getComponents()[i].setx(x[k][i]);
}
system.getPhases()[system.getNumberOfPhases() - 1].normalize();
multiPhaseTest = true;
system.setBeta(system.getNumberOfPhases() - 1,
system.getPhase(0).getComponent(unstabcomp).getz());
system.init(1);
system.normalizeBeta();

// logger.info("STABILITY ANALYSIS: ");
// logger.info("tm1: " + k + " "+ tm[k]);
// system.display();
return;
}
}

system.normalizeBeta();
// logger.info("STABILITY ANALYSIS: ");
// logger.info("tm1: " + tm[0] + " tm2: " + tm[1]);
// system.display();
}

/**
* <p>
* stabilityAnalysis2.
Expand Down Expand Up @@ -715,7 +1060,7 @@ public void stabilityAnalysis2() {
iter++;
err = 0;

if (iter <= 20 || !system.isImplementedCompositionDeriativesofFugacity()) {
if (iter <= 150 || !system.isImplementedCompositionDeriativesofFugacity()) {
if (iter % 7 == 0) {
double vec1 = 0.0;

Expand Down Expand Up @@ -1028,6 +1373,7 @@ public void run() {
*/
if (hasRemovedPhase && !secondTime) {
secondTime = true;
stabilityAnalysis3();
run();
}
/*
Expand Down
Loading

0 comments on commit af10b76

Please sign in to comment.