Skip to content

Commit

Permalink
initial test rachford rice (#995)
Browse files Browse the repository at this point in the history 
* initial test rachford rice

* update rachford rice - moved to separate class

* modified rachford rice implementation

* added Rachford rice tests
  • Loading branch information
@EvenSol
EvenSol authored May 1, 2024
1 parent c2e40f7 commit caf19b3
Show file tree
Hide file tree
Showing 6 changed files with 218 additions and 195 deletions.
33 changes: 18 additions & 15 deletions src/main/java/neqsim/thermo/system/SystemInterface.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -385,18 +385,6 @@ public void addTBPfraction(String componentName, double numberOfMoles, double mo
*/
public void calc_x_y_nonorm();

/**
* <p>
* calcBeta. For simple gas liquid systems.
* </p>
*
* @return Beta Mole fraction contained in the heaviest phase, i.e., liquid phase.
* @throws neqsim.util.exception.IsNaNException if any.
* @throws neqsim.util.exception.TooManyIterationsException if any.
*/
public double calcBeta()
throws neqsim.util.exception.IsNaNException, neqsim.util.exception.TooManyIterationsException;

/**
* <p>
* calcHenrysConstant.
Expand Down Expand Up @@ -2431,14 +2419,14 @@ public void setImplementedTemperatureDeriativesofFugacity(
*/
public void setPhaseType(int phaseToChange, PhaseType pt);

/**
/**
* Change the phase type of a given phase.
*
* @param phaseToChange the phase number of the phase to set phase type
* @param phaseName String to set
* @param phaseName String to set
*/
public void setPhaseType(int phaseToChange, String phaseName);

/**
* Set the physical property model type for each phase of the System.
*
Expand Down Expand Up @@ -2583,4 +2571,19 @@ public void setImplementedTemperatureDeriativesofFugacity(
* @param newfile a boolean
*/
public void write(String name, String filename, boolean newfile);

/**
* <p>
* getKvector - return vector of equilibrium constants
* </p>
*/
public double[] getKvector();


/**
* <p>
* getzvector - return vector of total molar composition
* </p>
*/
public double[] getzvector();
}
195 changes: 20 additions & 175 deletions src/main/java/neqsim/thermo/system/SystemThermo.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1241,181 +1241,6 @@ public final void calc_x_y_nonorm() {
}
}

/** {@inheritDoc} */
@Override
public final double calcBeta() throws neqsim.util.exception.IsNaNException,
neqsim.util.exception.TooManyIterationsException {
ComponentInterface[] compArray = getPhase(0).getComponents();

int i;
double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
double deriv = 0.0;
double gbeta = 0.0;
double betal = 0;
double nybeta = 0;

double midler = 0;
double minBeta = tolerance;
double maxBeta = 1.0 - tolerance;
double g0 = -1.0;
double g1 = 1.0;
nybeta = this.beta[0];
betal = 1.0 - nybeta;

for (i = 0; i < numberOfComponents; i++) {
midler = (compArray[i].getK() * compArray[i].getz() - 1.0) / (compArray[i].getK() - 1.0);
if ((midler > minBeta) && (compArray[i].getK() > 1.0)) {
minBeta = midler;
}
midler = (1.0 - compArray[i].getz()) / (1.0 - compArray[i].getK());
if ((midler < maxBeta) && (compArray[i].getK() < 1.0)) {
maxBeta = midler;
}
g0 += compArray[i].getz() * compArray[i].getK();
g1 += -compArray[i].getz() / compArray[i].getK();
}

if (g0 < 0) {
this.beta[1] = 1.0 - tolerance;
this.beta[0] = tolerance;
return this.beta[0];
}
if (g1 > 0) {
this.beta[1] = tolerance;
this.beta[0] = 1.0 - tolerance;
return this.beta[0];
}

nybeta = (minBeta + maxBeta) / 2.0;
// System.out.println("guessed beta: " + nybeta + " maxbeta: " +maxBeta + "
// minbeta: " +minBeta );
betal = 1.0 - nybeta;

// ' *l = 1.0-nybeta;
double gtest = 0.0;
for (i = 0; i < numberOfComponents; i++) {
gtest += compArray[i].getz() * (compArray[i].getK() - 1.0)
/ (1.0 - nybeta + nybeta * compArray[i].getK()); // beta
// =
// nybeta
}

if (gtest >= 0) {
minBeta = nybeta;
} else {
maxBeta = nybeta;
}

if (gtest < 0) {
double minold = minBeta;
minBeta = 1.0 - maxBeta;
maxBeta = 1.0 - minold;
}

int iterations = 0;
int maxIterations = 300;
// System.out.println("gtest: " + gtest);
double step = 1.0;
do {
iterations++;
if (gtest >= 0) {
// oldbeta = nybeta;
deriv = 0.0;
gbeta = 0.0;

for (i = 0; i < numberOfComponents; i++) {
double temp1 = (compArray[i].getK() - 1.0);
double temp2 = 1.0 + temp1 * nybeta;
deriv += -(compArray[i].getz() * temp1 * temp1) / (temp2 * temp2);
gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
/ (1.0 + (compArray[i].getK() - 1.0) * nybeta);
}

if (gbeta >= 0) {
minBeta = nybeta;
} else {
maxBeta = nybeta;
}
nybeta -= (gbeta / deriv);

// System.out.println("beta: " + maxBeta);
if (nybeta > maxBeta) {
nybeta = maxBeta;
}
if (nybeta < minBeta) {
nybeta = minBeta;
}

/*
* if ((nybeta > maxBeta) || (nybeta < minBeta)) { // nybeta = 0.5 * (maxBeta + minBeta);
* gbeta = 1.0; }
*/
} else {
// oldbeta = betal;
deriv = 0.0;
gbeta = 0.0;

for (i = 0; i < numberOfComponents; i++) {
deriv -= (compArray[i].getz() * (compArray[i].getK() - 1.0) * (1.0 - compArray[i].getK()))
/ Math.pow((betal + (1 - betal) * compArray[i].getK()), 2);
gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
/ (betal + (-betal + 1.0) * compArray[i].getK());
}

if (gbeta < 0) {
minBeta = betal;
} else {
maxBeta = betal;
}

betal -= (gbeta / deriv);

if (betal > maxBeta) {
betal = maxBeta;
}
if (betal < minBeta) {
betal = minBeta;
}

/*
* if ((betal > maxBeta) || (betal < minBeta)) { gbeta = 1.0; { betal = 0.5 * (maxBeta +
* minBeta); } }
*/
nybeta = 1.0 - betal;
}
step = gbeta / deriv;
// System.out.println("step : " + step);
} while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // &&
// (Math.abs(nybeta)-Math.abs(maxBeta))>0.1);

// System.out.println("beta: " + nybeta + " iterations: " + iterations);
if (nybeta <= tolerance) {
// this.phase = 1;
nybeta = tolerance;
} else if (nybeta >= 1.0 - tolerance) {
// this.phase = 0;
nybeta = 1.0 - tolerance;
// superheated vapour
} else {
// this.phase = 2;
} // two-phase liquid-gas

this.beta[0] = nybeta;
this.beta[1] = 1.0 - nybeta;

if (iterations >= maxIterations) {
throw new neqsim.util.exception.TooManyIterationsException(this, "calcBeta", maxIterations);
}
if (Double.isNaN(beta[1])) {
/*
* for (i = 0; i < numberOfComponents; i++) { System.out.println("K " + compArray[i].getK());
* System.out.println("z " + compArray[i].getz()); }
*/
throw new neqsim.util.exception.IsNaNException(this, "calcBeta", "beta");
}
return this.beta[0];
}

/** {@inheritDoc} */
@Override
public double calcHenrysConstant(String component) {
Expand Down Expand Up @@ -5196,4 +5021,24 @@ public void write(String name, String filename, boolean newfile) {
file.setValues(table);
file.createFile();
}

/** {@inheritDoc} */
@Override
public double[] getKvector() {
double[] K = new double[this.getNumberOfComponents()];
for (int i = 0; i < this.getNumberOfComponents(); i++) {
K[i] = this.getPhase(0).getComponent(i).getK();
}
return K;
}

/** {@inheritDoc} */
@Override
public double[] getzvector() {
double[] z = new double[this.getNumberOfComponents()];
for (int i = 0; i < this.getNumberOfComponents(); i++) {
z[i] = this.getPhase(0).getComponent(i).getz();
}
return z;
}
}
2 changes: 1 addition & 1 deletion src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -353,7 +353,7 @@ public boolean stabilityCheck() {
}
} else {
try {
system.calcBeta();
system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
} catch (Exception ex) {
logger.error(ex.getMessage(), ex);
}
Expand Down
Loading

0 comments on commit caf19b3

Please sign in to comment.