diff --git a/src/main/java/neqsim/thermodynamicOperations/ThermodynamicOperations.java b/src/main/java/neqsim/thermodynamicOperations/ThermodynamicOperations.java
index b36c97fd1c..36f5ce8163 100644
--- a/src/main/java/neqsim/thermodynamicOperations/ThermodynamicOperations.java
+++ b/src/main/java/neqsim/thermodynamicOperations/ThermodynamicOperations.java
@@ -59,6 +59,7 @@
import neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.CricondenThermFlash;
import neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.HPTphaseEnvelope;
import neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.pTphaseEnvelope;
+import neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.pTphaseEnvelopeNew2;
import neqsim.thermodynamicOperations.phaseEnvelopeOps.reactiveCurves.pLoadingCurve2;
import neqsim.thermodynamicOperations.propertyGenerator.OLGApropertyTableGeneratorWaterStudents;
import neqsim.thermodynamicOperations.propertyGenerator.OLGApropertyTableGeneratorWaterStudentsPH;
@@ -1539,6 +1540,23 @@ public void calcPTphaseEnvelope() {
getOperation().run();
}
+ // public void dewPointPressureFlash(){
+ // constantDutyFlashInterface operation = new constantDutyPressureFlash(system);
+ // operation.setBeta((1-1e-7));
+ // operation.run();
+ // }
+ /**
+ *
+ * calcPTphaseEnvelopeNew.
+ *
+ */
+ public void calcPTphaseEnvelope2() {
+ operation = new pTphaseEnvelopeNew2(system, fileName, (1.0 - 1e-10), 1.0, false);
+ // thisThread = new Thread(operation);
+ // thisThread.start();
+ getOperation().run();
+ }
+
/**
*
* calcPTphaseEnvelope.
@@ -2074,7 +2092,8 @@ public CalculationResult propertyFlash(List Spec1, List Spec2, i
}
if (hasOnlineFractions) {
- // Assure that fraction is inserted for the correct component (in case of mismatch of
+ // Assure that fraction is inserted for the correct component (in case of
+ // mismatch of
// component input and fluid component list)
double[] fraction = new double[this.system.getNumberOfComponents()];
// For all components defined in system
diff --git a/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/pTphaseEnvelopeNew2.java b/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/pTphaseEnvelopeNew2.java
new file mode 100644
index 0000000000..35712d84ff
--- /dev/null
+++ b/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/pTphaseEnvelopeNew2.java
@@ -0,0 +1,599 @@
+/*
+ * pTphaseEnvelope.java
+ *
+ * Created on 14. oktober 2000, 21:59 Updated on May 2019, by Nefeli
+ */
+
+package neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps;
+
+import java.util.ArrayList;
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+import neqsim.thermo.system.SystemInterface;
+import neqsim.thermodynamicOperations.BaseOperation;
+import neqsim.thermodynamicOperations.ThermodynamicOperations;
+
+/**
+ *
+ * pTphaseEnvelopeNew2 class.
+ *
+ *
+ * @author asmund
+ * @version $Id: $Id
+ */
+public class pTphaseEnvelopeNew2 extends BaseOperation {
+ private static final long serialVersionUID = 1000;
+ static Logger logger = LogManager.getLogger(pTphaseEnvelopeNew2.class);
+
+ double maxPressure = 1000.0;
+ double minPressure = 1.0;
+ SystemInterface system;
+ boolean bubblePointFirst = true;
+ boolean calculatesDewPoint = true;
+ double[] cricondenTherm = new double[3];
+ double[] cricondenBar = new double[3];
+ double[] cricondenThermX = new double[100];
+ double[] cricondenThermY = new double[100];
+ double[] cricondenBarX = new double[100];
+ double[] cricondenBarY = new double[100];
+ double phaseFraction = 1e-10;
+ int i;
+ int j = 0;
+ int nummer = 0;
+ int iterations = 0;
+ int maxNumberOfIterations = 10000;
+ double gibbsEnergy = 0;
+ double gibbsEnergyOld = 0;
+ double Kold;
+ double deviation = 0;
+ double g0 = 0;
+ double g1 = 0;
+ double lowPres = 1.0;
+ double[] lnOldOldK;
+ double[] lnK;
+ boolean outputToFile = false;
+ double[] lnOldK;
+ double[] lnKwil;
+ double[] oldDeltalnK;
+ double[] deltalnK;
+ double[] tm = {1, 1};
+ double beta = 1e-5;
+ int lowestGibbsEnergyPhase = 0;
+ String fileName = "c:/file";
+ double temp = 0;
+ double pres = 0;
+ double startPres = 0;
+ boolean moreLines = false;
+ boolean restart = true;
+ int np = 0;
+ // points[2] = new double[1000];
+ int speceq = 0;
+ double Tcfirst;
+ double Pcfirst;
+ double Tmin = 0.0;
+
+ ArrayList dewPointTemperature = new ArrayList();
+ ArrayList dewPointPressure = new ArrayList();
+ ArrayList bubblePointTemperature = new ArrayList();
+ ArrayList bubblePointPressure = new ArrayList();
+
+ ArrayList bubblePointEnthalpy = new ArrayList();
+ ArrayList dewPointEnthalpy = new ArrayList();
+
+ ArrayList bubblePointEntropy = new ArrayList();
+ ArrayList dewPointEntropy = new ArrayList();
+
+ ArrayList bubblePointVolume = new ArrayList();
+ ArrayList dewPointVolume = new ArrayList();
+
+ double[] cricondenThermfirst = new double[3];
+ double[] cricondenBarfirst = new double[3];
+ double[] cricondenThermXfirst = new double[100];
+ double[] cricondenThermYfirst = new double[100];
+ double[] cricondenBarXfirst = new double[100];
+ double[] cricondenBarYfirst = new double[100];
+
+ double[] dewPointTemperatureArray;
+ double[] dewPointPressureArray;
+ double[] dewPointEnthalpyArray;
+ double[] dewPointVolumeArray;
+ double[] dewPointEntropyArray;
+
+ double[] bubblePointTemperatureArray;
+ double[] bubblePointPressureArray;
+ double[] bubblePointEnthalpyArray;
+ double[] bubblePointVolumeArray;
+ double[] bubblePointEntropyArray;
+
+ /**
+ *
+ * Constructor for pTphaseEnvelope.
+ *
+ */
+ public pTphaseEnvelopeNew2() {}
+
+ /**
+ *
+ * Constructor for pTphaseEnvelope.
+ *
+ *
+ * @param system a {@link neqsim.thermo.system.SystemInterface} object
+ * @param name a {@link java.lang.String} object
+ * @param phaseFraction a double
+ * @param lowPres a double
+ * @param bubfirst a boolean
+ */
+ public pTphaseEnvelopeNew2(SystemInterface system, String name, double phaseFraction,
+ double lowPres, boolean bubfirst) {
+ this.bubblePointFirst = bubfirst;
+ if (name != null) {
+ outputToFile = true;
+ fileName = name;
+ }
+ this.system = system;
+ this.phaseFraction = phaseFraction;
+ lnOldOldK = new double[system.getPhase(0).getNumberOfComponents()];
+ lnOldK = new double[system.getPhase(0).getNumberOfComponents()];
+ lnK = new double[system.getPhase(0).getNumberOfComponents()];
+ this.lowPres = lowPres;
+ oldDeltalnK = new double[system.getPhase(0).getNumberOfComponents()];
+ deltalnK = new double[system.getPhase(0).getNumberOfComponents()];
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public void run() {
+ speceq = 0; // initialization
+ try {
+ system.init(0); // initialization
+
+ // selects the most volatile and least volatile component based on Tc values
+ // afterwards it uses them to define the speceq of the first point
+ // based on the desired first point, dew/bubble
+ for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
+ if (system.getComponent(i).getz() < 1e-10) {
+ continue;
+ }
+ if (system.getPhase(0).getComponent(i).getIonicCharge() == 0) {
+ if (bubblePointFirst && system.getPhase(0).getComponents()[speceq]
+ .getTC() > system.getPhase(0).getComponents()[i].getTC()) {
+ speceq = system.getPhase(0).getComponent(i).getComponentNumber();
+ }
+ if (!bubblePointFirst && system.getPhase(0).getComponents()[speceq]
+ .getTC() < system.getPhase(0).getComponents()[i].getTC()) {
+ speceq = system.getPhase(0).getComponent(i).getComponentNumber();
+ }
+ }
+ }
+
+ // initialized the first step of the phase envelope
+ // pressure is already defined
+ // temperature is the antoine vapor pressure of the selected component
+ // (least or most volatile.
+ pres = lowPres;
+ // temp =
+ // system.getPhase(0).getComponent(speceq).getAntoineVaporTemperature(pres);
+ temp = tempKWilson(phaseFraction, pres);
+
+ if (Double.isNaN(temp)) {
+ temp = system.getPhase(0).getComponent(speceq).getTC() - 20.0;
+ }
+ system.setTemperature(temp);
+ system.setPressure(pres);
+
+ ThermodynamicOperations testOps = new ThermodynamicOperations(system);
+
+ // this part converges the first phase envelope point.
+ // if the plasefraction is more than 0.5 it does a dew point initiallization
+ // else a bubble point initiallization
+
+ for (int i = 0; i < 5; i++) {
+ try {
+ if (phaseFraction < 0.5) {
+ temp += i * 2;
+ system.setTemperature(temp);
+ testOps.bubblePointTemperatureFlash();
+ } else {
+ temp += i * 2;
+ system.setTemperature(temp);
+ testOps.dewPointTemperatureFlash();
+ }
+ } catch (Exception ex) {
+ // ex.toString();
+ }
+ double tempNy = system.getTemperature();
+
+ if (!Double.isNaN(tempNy)) {
+ temp = tempNy;
+ break;
+ }
+ }
+
+ // this part sets the first envelope point into the system
+ system.setBeta(phaseFraction);
+ system.setPressure(pres);
+ system.setTemperature(temp);
+
+ sysNewtonRhapsonPhaseEnvelope nonLinSolver =
+ new sysNewtonRhapsonPhaseEnvelope(system, 2, system.getPhase(0).getNumberOfComponents());
+ startPres = system.getPressure();
+ nonLinSolver.setu();
+
+ for (np = 1; np < 9980; np++) {
+ try {
+ // solves the np point of the envelope
+ nonLinSolver.calcInc(np);
+ nonLinSolver.solve(np);
+
+ // this catches the exceptions
+ // double TT = system.getPhase(0).getTemperature();
+ // double PP = system.getPhase(0).getPressure();
+ } catch (Exception e0) {
+ // the envelope crushed.
+ // this part keeps the old values
+ // restarts the envelope from the other side
+ // and then stops
+
+ if (restart) {
+ calculatesDewPoint = false;
+ restart = !restart;
+ Tcfirst = system.getTC();
+ Pcfirst = system.getPC();
+
+ cricondenBarfirst = cricondenBar;
+ cricondenBarXfirst = cricondenBarX;
+ cricondenBarYfirst = cricondenBarY;
+
+ cricondenThermfirst = cricondenTherm;
+ cricondenThermXfirst = cricondenThermX;
+ cricondenThermYfirst = cricondenThermY;
+
+ // new settings
+ phaseFraction = 1.0 - phaseFraction;
+ bubblePointFirst = !bubblePointFirst;
+ run();
+ /**/
+ break;
+ } else {
+ np = np - 1;
+ break;
+ }
+ }
+
+ // check for critical point
+ double Kvallc = system.getPhase(0).getComponent(nonLinSolver.lc).getx()
+ / system.getPhase(1).getComponent(nonLinSolver.lc).getx();
+ double Kvalhc = system.getPhase(0).getComponent(nonLinSolver.hc).getx()
+ / system.getPhase(1).getComponent(nonLinSolver.hc).getx();
+ // double densV = system.getPhase(0).getDensity();
+ // double densL = system.getPhase(1).getDensity();
+
+ // System.out.println(np + " " + system.getTemperature() + " " +
+ // system.getPressure() + " " + densV + " " + densL );
+
+ if (!nonLinSolver.etterCP) {
+ if (Kvallc < 1.05 && Kvalhc > 0.95) {
+ // close to the critical point
+ // invert phase types and find the CP Temp and Press
+
+ // System.out.println("critical point");
+ nonLinSolver.npCrit = np;
+ system.invertPhaseTypes();
+ nonLinSolver.etterCP = true;
+ calculatesDewPoint = false;
+ // the critical point is found from interpolation polynomials based on K=1 of
+ // the most or least volatile component
+ nonLinSolver.calcCrit();
+ }
+ }
+
+ // stores critondenbar and cricondentherm
+ // HERE the new cricoT and crico P values will be called instead
+ if (system.getTemperature() > cricondenTherm[0]) {
+ cricondenTherm[1] = system.getPressure();
+ cricondenTherm[0] = system.getTemperature();
+ for (int ii = 0; ii < nonLinSolver.numberOfComponents; ii++) {
+ cricondenThermX[ii] = system.getPhase(1).getComponent(ii).getx();
+ cricondenThermY[ii] = system.getPhase(0).getComponent(ii).getx();
+ }
+ } else {
+ nonLinSolver.ettercricoT = true;
+ }
+ if (system.getPressure() > cricondenBar[1]) {
+ cricondenBar[0] = system.getTemperature();
+ cricondenBar[1] = system.getPressure();
+ for (int ii = 0; ii < nonLinSolver.numberOfComponents; ii++) {
+ cricondenBarX[ii] = system.getPhase(1).getComponent(ii).getx();
+ cricondenBarY[ii] = system.getPhase(0).getComponent(ii).getx();
+ }
+ }
+
+ // Exit criteria
+ if ((system.getPressure() < minPressure && nonLinSolver.ettercricoT)) {
+ break;
+ }
+ if (system.getPressure() > maxPressure) {
+ break;
+ }
+ if (system.getTemperature() < Tmin) {
+ break;
+ }
+
+ if (system.getTemperature() > 1e-6 && system.getPressure() > 1e-6
+ && !(Double.isNaN(system.getTemperature()) || Double.isNaN(system.getPressure()))) {
+ if (calculatesDewPoint) {
+ dewPointTemperature.add(system.getTemperature());
+ dewPointPressure.add(system.getPressure());
+ dewPointEnthalpy
+ .add(system.getPhase(1).getEnthalpy() / system.getPhase(1).getNumberOfMolesInPhase()
+ / system.getPhase(1).getMolarMass() / 1e3);
+ dewPointVolume.add(system.getPhase(1).getDensity());
+ dewPointEntropy
+ .add(system.getPhase(1).getEntropy() / system.getPhase(1).getNumberOfMolesInPhase()
+ / system.getPhase(1).getMolarMass() / 1e3);
+ } else {
+ bubblePointTemperature.add(system.getTemperature());
+ bubblePointPressure.add(system.getPressure());
+ bubblePointEnthalpy
+ .add(system.getPhase(1).getEnthalpy() / system.getPhase(1).getNumberOfMolesInPhase()
+ / system.getPhase(1).getMolarMass() / 1e3);
+ bubblePointVolume.add(system.getPhase(1).getDensity());
+ bubblePointEntropy
+ .add(system.getPhase(1).getEntropy() / system.getPhase(1).getNumberOfMolesInPhase()
+ / system.getPhase(1).getMolarMass() / 1e3);
+ }
+ }
+ }
+
+ system.setTemperature(system.getTC());
+ system.setPressure(system.getPC());
+ } catch (Exception ex) {
+ logger.error(ex.getMessage(), ex);
+ throw ex;
+ }
+
+ dewPointTemperatureArray = new double[dewPointTemperature.size()];
+ dewPointPressureArray = new double[dewPointPressure.size()];
+ dewPointEnthalpyArray = new double[dewPointTemperature.size()];
+ dewPointVolumeArray = new double[dewPointPressure.size()];
+ dewPointEntropyArray = new double[dewPointPressure.size()];
+
+ bubblePointTemperatureArray = new double[bubblePointTemperature.size()];
+ bubblePointPressureArray = new double[bubblePointPressure.size()];
+ bubblePointEnthalpyArray = new double[bubblePointPressure.size()];
+ bubblePointVolumeArray = new double[bubblePointPressure.size()];
+ bubblePointEntropyArray = new double[bubblePointPressure.size()];
+
+ for (int i = 0; i < dewPointTemperature.size(); i++) {
+ dewPointTemperatureArray[i] = dewPointTemperature.get(i);
+ dewPointPressureArray[i] = dewPointPressure.get(i);
+ dewPointEnthalpyArray[i] = dewPointEnthalpy.get(i);
+ dewPointVolumeArray[i] = dewPointVolume.get(i);
+ dewPointEntropyArray[i] = dewPointEntropy.get(i);
+ }
+
+ for (int i = 0; i < bubblePointTemperature.size(); i++) {
+ bubblePointTemperatureArray[i] = bubblePointTemperature.get(i);
+ bubblePointPressureArray[i] = bubblePointPressure.get(i);
+ bubblePointEnthalpyArray[i] = bubblePointEnthalpy.get(i);
+ bubblePointEntropyArray[i] = bubblePointEntropy.get(i);
+ bubblePointVolumeArray[i] = bubblePointVolume.get(i);
+ }
+ }
+
+ /**
+ *
+ * calcHydrateLine.
+ *
+ */
+ public void calcHydrateLine() {
+ ThermodynamicOperations opsHyd = new ThermodynamicOperations(system);
+ try {
+ opsHyd.hydrateEquilibriumLine(10.0, 300.0);
+ } catch (Exception ex) {
+ logger.error(ex.getMessage(), ex);
+ }
+
+ // double[][] hydData = opsHyd.getData();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public void printToFile(String name) {}
+
+
+ /** {@inheritDoc} */
+ @Override
+ public double[] get(String name) {
+ if (name.equals("dewT")) {
+ return dewPointTemperatureArray;
+ }
+ if (name.equals("dewP")) {
+ return dewPointPressureArray;
+ }
+ if (name.equals("bubT")) {
+ return bubblePointTemperatureArray;
+ }
+ if (name.equals("bubP")) {
+ return bubblePointPressureArray;
+ }
+
+ if (name.equals("dewH")) {
+ return dewPointEnthalpyArray;
+ }
+ if (name.equals("dewDens")) {
+ return dewPointVolumeArray;
+ }
+ if (name.equals("dewS")) {
+ return dewPointEntropyArray;
+ }
+ if (name.equals("bubH")) {
+ return bubblePointEnthalpyArray;
+ }
+ if (name.equals("bubDens")) {
+ return bubblePointVolumeArray;
+ }
+ if (name.equals("bubS")) {
+ return bubblePointEntropyArray;
+ }
+ if (name.equals("cricondentherm")) {
+ return cricondenTherm;
+ }
+ if (name.equals("cricondenthermX")) {
+ return cricondenThermX;
+ }
+ if (name.equals("cricondenthermY")) {
+ return cricondenThermY;
+ }
+ if (name.equals("cricondenbar")) {
+ return cricondenBar;
+ }
+ if (name.equals("cricondenbarX")) {
+ return cricondenBarX;
+ }
+ if (name.equals("cricondenbarY")) {
+ return cricondenBarY;
+ }
+ if (name.equals("criticalPoint1")) {
+ return new double[] {system.getTC(), system.getPC()};
+ }
+ if (name.equals("criticalPoint2")) {
+ return new double[] {0, 0};
+ } else {
+ return null;
+ }
+ }
+
+ /**
+ * Getter for property bubblePointFirst.
+ *
+ * @return Value of property bubblePointFirst.
+ */
+ public boolean isBubblePointFirst() {
+ return bubblePointFirst;
+ }
+
+ /**
+ * Setter for property bubblePointFirst.
+ *
+ * @param bubblePointFirst New value of property bubblePointFirst.
+ */
+ public void setBubblePointFirst(boolean bubblePointFirst) {
+ this.bubblePointFirst = bubblePointFirst;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String[][] getResultTable() {
+ return null;
+ }
+
+ /**
+ *
+ * tempKWilson.
+ *
+ *
+ * @param beta a double
+ * @param P a double
+ * @return a double
+ */
+ public double tempKWilson(double beta, double P) {
+ // Initiallizes the temperature of a saturation point for given pressure
+ // based on K values of Wilson
+ // see Michelsen book thermodynamics & computational aspects
+
+ double initTc = 0.;
+ double initPc = 0.;
+ double initAc = 0.;
+ double Tstart = 0.;
+ double Tstartold = 0.;
+ double initT = 0;
+ double dinitT = 0;
+ int numberOfComponents = system.getPhase(0).getNumberOfComponents();
+ int lc = 0;
+ int hc = 0;
+
+ double[] Kwil = new double[numberOfComponents];
+
+ double min = 100000.;
+ double max = 0.;
+
+ for (int i = 0; i < numberOfComponents; i++) {
+ if (system.getPhase(0).getComponents()[i].getTC() > max) {
+ max = system.getPhase(0).getComponents()[i].getTC();
+ hc = i;
+ }
+ if (system.getPhase(0).getComponents()[i].getTC() < min) {
+ min = system.getPhase(0).getComponents()[i].getTC();
+ lc = i;
+ }
+ }
+
+ try {
+ if (beta <= 0.5) {
+ // closer to bubble point get the lightest component
+
+ initTc = system.getPhase(0).getComponents()[lc].getTC();
+ initPc = system.getPhase(0).getComponents()[lc].getPC();
+ initAc = system.getPhase(0).getComponents()[lc].getAcentricFactor();
+ } else if (beta > 0.5) {
+ // closer to dew point get the heaviest component
+ initTc = system.getPhase(0).getComponents()[hc].getTC();
+ initPc = system.getPhase(0).getComponents()[hc].getPC();
+ initAc = system.getPhase(0).getComponents()[hc].getAcentricFactor();
+ }
+
+ // initial T based on the lightest/heaviest component
+ Tstart = initTc * 5.373 * (1 + initAc) / (5.373 * (1 + initAc) - Math.log(P / initPc));
+
+ // solve for Tstart with Newton
+ for (int i = 0; i < 1000; i++) {
+ initT = 0.;
+ dinitT = 0.;
+ for (int j = 0; j < numberOfComponents; j++) {
+ Kwil[j] = system.getPhase(0).getComponents()[j].getPC() / P
+ * Math.exp(5.373 * (1. + system.getPhase(0).getComponents()[j].getAcentricFactor())
+ * (1. - system.getPhase(0).getComponents()[j].getTC() / Tstart));
+ // system.getPhases()[0].getComponents()[j].setK(Kwil[j]);
+ }
+
+ for (int j = 0; j < numberOfComponents; j++) {
+ if (beta < 0.5) {
+ initT = initT + system.getPhase(0).getComponents()[j].getz() * Kwil[j];
+ dinitT = dinitT + system.getPhase(0).getComponents()[j].getz() * Kwil[j] * 5.373
+ * (1 + system.getPhase(0).getComponents()[j].getAcentricFactor())
+ * system.getPhase(0).getComponents()[j].getTC() / (Tstart * Tstart);
+ } else {
+ initT = initT + system.getPhase(0).getComponents()[j].getz() / Kwil[j];
+ dinitT = dinitT - system.getPhase(0).getComponents()[j].getz() / Kwil[j] * 5.373
+ * (1 + system.getPhase(0).getComponents()[j].getAcentricFactor())
+ * system.getPhase(0).getComponents()[j].getTC() / (Tstart * Tstart);
+ }
+ }
+
+ initT = initT - 1.;
+ if (Math.abs(initT / dinitT) > 0.1 * Tstart) {
+ Tstart = Tstart - 0.001 * initT / dinitT;
+ } else {
+ Tstart = Tstart - initT / dinitT;
+ }
+ if (Math.abs(Tstart - Tstartold) < 1.e-5) {
+ return Tstart;
+ }
+ Tstartold = Tstart;
+ }
+ } catch (Exception ex) {
+ Tstart = initTc * 5.373 * (1 + initAc) / (5.373 * (1 + initAc) - Math.log(P / initPc));
+ }
+ if (Double.isNaN(Tstart) || Double.isInfinite(Tstart)) {
+ Tstart = initTc * 5.373 * (1 + initAc) / (5.373 * (1 + initAc) - Math.log(P / initPc));
+ }
+ return Tstart;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public void displayResult() {
+
+ }
+}
diff --git a/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/sysNewtonRhapsonPhaseEnvelope.java b/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/sysNewtonRhapsonPhaseEnvelope.java
index ba7060e99f..71c3265a19 100644
--- a/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/sysNewtonRhapsonPhaseEnvelope.java
+++ b/src/main/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/sysNewtonRhapsonPhaseEnvelope.java
@@ -413,27 +413,6 @@ public void calcInc(int np) {
// calculate the dxds of the system
dxds = Jac.solve(fvec);
- // check for critical point
-
- // check density
- // double densV = system.getPhase(0).getDensity();
- // double densL = system.getPhase(1).getDensity();
- // check the proximity to the critical point by adding the lnKs and finding the highest
- double Kvallc =
- system.getPhase(0).getComponent(lc).getx() / system.getPhase(1).getComponent(lc).getx();
- double Kvalhc =
- system.getPhase(0).getComponent(hc).getx() / system.getPhase(1).getComponent(hc).getx();
-
- if (!etterCP) {
- if (Kvallc < 1.05 && Kvalhc > 0.95) {
- calcCP = true;
- etterCP = true;
- npCrit = np;
- system.invertPhaseTypes();
- // System.out.println("critical point");
- }
- }
-
// manipulate stepsize according to the number of iterations of the previous point
if (iter > 6) {
ds *= 0.5;
diff --git a/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTPhaseEnvelopeTest.java b/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTPhaseEnvelopeTest.java
new file mode 100644
index 0000000000..887f064f59
--- /dev/null
+++ b/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTPhaseEnvelopeTest.java
@@ -0,0 +1,178 @@
+package neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps;
+
+import static org.junit.jupiter.api.Assertions.assertArrayEquals;
+import static org.junit.jupiter.api.Assertions.assertThrows;
+import static org.junit.jupiter.api.Assertions.assertTrue;
+import org.junit.jupiter.api.BeforeEach;
+import org.junit.jupiter.api.Test;
+import neqsim.thermodynamicOperations.ThermodynamicOperations;
+
+public class PTPhaseEnvelopeTest {
+ static neqsim.thermo.system.SystemInterface testSystem = null;
+ static ThermodynamicOperations testOps = null;
+
+ @BeforeEach
+ void setUp() {
+ testSystem = new neqsim.thermo.system.SystemSrkEos(298.0, 50.0);
+ }
+
+ /**
+ * Test method for
+ * {@link neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.pTphaseEnvelope}.
+ */
+ @Test
+ void testDewP() {
+ testSystem.addComponent("nitrogen", 0.01);
+ testSystem.addComponent("CO2", 0.01);
+ testSystem.addComponent("methane", 0.98);
+ testSystem.setMixingRule("classic");
+
+ testOps = new ThermodynamicOperations(testSystem);
+ testOps.TPflash();
+ testSystem.initProperties();
+ testOps.calcPTphaseEnvelope();
+ double[] dewPointPressures = testOps.get("dewP");
+ double[] expectedDewPointPressures =
+ new double[] {1.1051709180756477, 1.2214027581601699, 1.3498588075760032,
+ 1.4918246976412703, 1.6652911949458864, 1.8794891289619104, 2.1418131227502055,
+ 2.4690864123141987, 2.881197018974799, 3.404779997613969, 4.075230307874481,
+ 4.938583914869986, 6.051801019586486, 7.477304695462727, 9.260793952051571,
+ 11.364101185282063, 13.480106047577934, 14.53423776629387, 13.607498029406681,
+ 11.181207439509638, 9.189487040488075, 9.612827246459474, 10.706126846063928,
+ 12.501491987760147, 15.075672692089958, 18.51283799420178, 23.330378296334104,
+ 29.71319711031059, 37.25532259549197, 43.660805656603934, 45.75836660678656,
+ 46.42490219574348, 46.83203503669948, 46.869568345957006, 46.903557772489435};
+ // System.out.println(Arrays.toString(dewPointPressures));
+ assertArrayEquals(expectedDewPointPressures, dewPointPressures, 10E-10);
+ }
+
+ @Test
+ void testFailingCaseWithWater() {
+ testSystem.addComponent("nitrogen", 0.04);
+ testSystem.addComponent("CO2", 0.06);
+ testSystem.addComponent("methane", 0.80);
+ testSystem.addComponent("water", 0.00000000001);
+
+ testSystem.setMixingRule("classic");
+
+ testOps = new ThermodynamicOperations(testSystem);
+ testOps.TPflash();
+ testSystem.initProperties();
+
+ Exception exception =
+ assertThrows(ArrayIndexOutOfBoundsException.class, () -> testOps.calcPTphaseEnvelope());
+ }
+
+ @Test
+ void testSimpleCase() {
+ testSystem.addComponent("nitrogen", 0.88);
+ testSystem.addComponent("CO2", 5.7);
+ testSystem.addComponent("methane", 86.89);
+ testSystem.addComponent("ethane", 3.59);
+ testSystem.addComponent("propane", 1.25);
+ testSystem.addComponent("i-butane", 0.19);
+ testSystem.addComponent("n-butane", 0.35);
+ testSystem.addComponent("i-pentane", 0.12);
+ testSystem.addComponent("n-pentane", 0.12);
+ testSystem.setMixingRule("classic");
+ testSystem.setMixingRule("classic");
+
+ testOps = new ThermodynamicOperations(testSystem);
+ testOps.calcPTphaseEnvelope2();
+ double[] dewPointPressures = testOps.get("dewP");
+ double[] dewPointTemperautres = testOps.get("dewT");
+ double[] bubblePointPressures = testOps.get("bubP");
+ double[] bubblePointTemperautres = testOps.get("bubT");
+ double[] cricondenbar = testOps.get("cricondenbar");
+ double[] criticalPoint1 = testOps.get("criticalPoint1");
+
+
+ assertTrue(dewPointTemperautres.length > 20);
+ assertTrue(bubblePointTemperautres.length > 10);
+
+ }
+
+
+ @Test
+ void testFailingCase1() {
+ // testSystem.setTemperature(40, "C");
+ // testSystem.setPressure(50, "bara");
+ testSystem.addComponent("nitrogen", 0.88);
+ testSystem.addComponent("CO2", 5.7);
+ testSystem.addComponent("methane", 86.89);
+ testSystem.addComponent("ethane", 3.59);
+ testSystem.addComponent("propane", 1.25);
+ testSystem.addComponent("i-butane", 0.19);
+ testSystem.addComponent("n-butane", 0.35);
+ testSystem.addComponent("i-pentane", 0.12);
+ testSystem.addComponent("n-pentane", 0.12);
+ testSystem.addTBPfraction("C6", 0.15, 86 / 1000.0, 0.672);
+ testSystem.addTBPfraction("C7", 0.2, 96 / 1000.0, 0.737);
+ testSystem.addTBPfraction("C8", 0.22, 106 / 1000.0, 0.767);
+ testSystem.addTBPfraction("C9", 0.13, 121 / 1000.0, 0.783);
+ testSystem.addPlusFraction("C10+", 0.21, 172 / 1000.0, 0.818);
+ testSystem.setMixingRule("classic");
+ testSystem.setMultiPhaseCheck(true);
+ testSystem.useVolumeCorrection(true);
+ testSystem.initPhysicalProperties();
+ testSystem.setMixingRule("classic");
+
+ testOps = new ThermodynamicOperations(testSystem);
+ testOps.TPflash();
+ testSystem.initProperties();
+
+ testOps.calcPTphaseEnvelope2();
+ double[] dewPointPressures = testOps.get("dewP");
+ double[] dewPointTemperautres = testOps.get("dewT");
+ double[] bubblePointPressures = testOps.get("bubP");
+ double[] bubblePointTemperautres = testOps.get("bubT");
+ double[] cricondenbar = testOps.get("cricondenbar");
+ double[] criticalPoint1 = testOps.get("criticalPoint1");
+
+
+ assertTrue(dewPointTemperautres.length > 20);
+ assertTrue(bubblePointTemperautres.length > 10);
+
+ }
+
+ @Test
+ void testFailingCase2() {
+ // testSystem.setTemperature(40, "C");
+ // testSystem.setPressure(50, "bara");
+ neqsim.thermo.system.SystemInterface fluid0_HC =
+ new neqsim.thermo.system.SystemUMRPRUMCEos(298.0, 50.0);
+ fluid0_HC.addComponent("nitrogen", 2.5);
+ fluid0_HC.addComponent("CO2", 4.5);
+ fluid0_HC.addComponent("methane", 79.45);
+ fluid0_HC.addComponent("ethane", 10);
+ fluid0_HC.addComponent("propane", 2.5);
+ fluid0_HC.addComponent("i-butane", 0.3);
+ fluid0_HC.addComponent("n-butane", 0.5);
+ fluid0_HC.addComponent("22-dim-C3", 0.01);
+ fluid0_HC.addComponent("i-pentane", 0.05);
+ fluid0_HC.addComponent("n-pentane", 0.05);
+ fluid0_HC.addComponent("n-hexane", 0.05);
+ fluid0_HC.addComponent("benzene", 0.02);
+ fluid0_HC.addComponent("c-hexane", 0.02);
+ fluid0_HC.addComponent("n-heptane", 0.02);
+ fluid0_HC.addComponent("toluene", 0.01);
+ fluid0_HC.addComponent("n-octane", 0.01);
+ fluid0_HC.setMixingRule("HV", "UNIFAC_UMRPRU");
+ testOps = new ThermodynamicOperations(fluid0_HC);
+ testOps.calcPTphaseEnvelope2();
+ double[] dewPointPressures = testOps.get("dewP");
+ double[] dewPointTemperautres = testOps.get("dewT");
+ double[] bubblePointPressures = testOps.get("bubP");
+ double[] bubblePointTemperautres = testOps.get("bubT");
+ double[] bubblePointEnthalpies = testOps.get("bubH");
+ double[] bubblePointVolumes = testOps.get("bubDens");
+ double[] cricondenbar = testOps.get("cricondenbar");
+ double[] criticalPoint1 = testOps.get("criticalPoint1");
+
+ assertTrue(dewPointTemperautres.length > 20);
+ assertTrue(bubblePointTemperautres.length > 20);
+ assertTrue(bubblePointEnthalpies.length > 20);
+ assertTrue(bubblePointVolumes.length > 20);
+
+ }
+}
diff --git a/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTphaseEnvelopeTest.java b/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTphaseEnvelopeTest.java
deleted file mode 100644
index 919a6d7340..0000000000
--- a/src/test/java/neqsim/thermodynamicOperations/phaseEnvelopeOps/multicomponentEnvelopeOps/PTphaseEnvelopeTest.java
+++ /dev/null
@@ -1,64 +0,0 @@
-package neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps;
-
-import static org.junit.jupiter.api.Assertions.assertArrayEquals;
-import static org.junit.jupiter.api.Assertions.assertThrows;
-import org.junit.jupiter.api.BeforeEach;
-import org.junit.jupiter.api.Test;
-import neqsim.thermodynamicOperations.ThermodynamicOperations;
-
-class PTPhaseEnvelopeTest {
- static neqsim.thermo.system.SystemInterface testSystem = null;
- static ThermodynamicOperations testOps = null;
-
- @BeforeEach
- void setUp() {
- testSystem = new neqsim.thermo.system.SystemSrkEos(298.0, 50.0);
- }
-
- /**
- * Test method for
- * {@link neqsim.thermodynamicOperations.phaseEnvelopeOps.multicomponentEnvelopeOps.pTphaseEnvelope}.
- */
- @Test
- void testDewP() {
- testSystem.addComponent("nitrogen", 0.01);
- testSystem.addComponent("CO2", 0.01);
- testSystem.addComponent("methane", 0.98);
- testSystem.setMixingRule("classic");
-
- testOps = new ThermodynamicOperations(testSystem);
- testOps.TPflash();
- testSystem.initProperties();
- testOps.calcPTphaseEnvelope();
- double[] dewPointPressures = testOps.get("dewP");
- double[] expectedDewPointPressures =
- new double[] {1.1051709180756477, 1.2214027581601699, 1.3498588075760032,
- 1.4918246976412703, 1.6652911949458864, 1.8794891289619104, 2.1418131227502055,
- 2.4690864123141987, 2.881197018974799, 3.404779997613969, 4.075230307874481,
- 4.938583914869986, 6.051801019586486, 7.477304695462727, 9.260793952051571,
- 11.364101185282063, 13.480106047577934, 14.53423776629387, 13.607498029406681,
- 11.181207439509638, 9.189487040488075, 9.612827246459474, 10.706126846063928,
- 12.501491987760147, 15.075672692089958, 18.51283799420178, 23.330378296334104,
- 29.71319711031059, 37.25532259549197, 43.660805656603934, 45.75836660678656,
- 46.42490219574348, 46.83203503669948, 46.869568345957006, 46.903557772489435};
- // System.out.println(Arrays.toString(dewPointPressures));
- assertArrayEquals(expectedDewPointPressures, dewPointPressures, 10E-10);
- }
-
- @Test
- void testFailingCaseWithWater() {
- testSystem.addComponent("nitrogen", 0.04);
- testSystem.addComponent("CO2", 0.06);
- testSystem.addComponent("methane", 0.80);
- testSystem.addComponent("water", 0.00000000001);
-
- testSystem.setMixingRule("classic");
-
- testOps = new ThermodynamicOperations(testSystem);
- testOps.TPflash();
- testSystem.initProperties();
-
- Exception exception =
- assertThrows(ArrayIndexOutOfBoundsException.class, () -> testOps.calcPTphaseEnvelope());
- }
-}