diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp index 6512c41b3..8559e6be0 100644 --- a/opm/material/eos/PengRobinson.hpp +++ b/opm/material/eos/PengRobinson.hpp @@ -191,10 +191,13 @@ public: // the EOS has three intersections with the pressure, // i.e. the molar volume of gas is the largest one and the // molar volume of liquid is the smallest one +#warning HACK, should investigate why if (isGasPhase) - Vm = Z[2]*RT/p; + // Vm = Z[2]*RT/p; + Vm = max(1e-7, Z[2]*RT/p); else - Vm = Z[0]*RT/p; + // Vm = Z[2]*RT/p; + Vm = max(1e-7, Z[0]*RT/p); } else if (numSol == 1) { // the EOS only has one intersection with the pressure, @@ -203,27 +206,28 @@ public: Evaluation VmCubic = Z[0]*RT/p; Vm = VmCubic; +#warning, should investigate why here // find the extrema (if they are present) - Evaluation Vmin, Vmax, pmin, pmax; - if (findExtrema_(Vmin, Vmax, - pmin, pmax, - a, b, T)) - { - if (isGasPhase) - Vm = std::max(Vmax, VmCubic); - else { - if (Vmin > 0) - Vm = std::min(Vmin, VmCubic); - else - Vm = VmCubic; - } - } - else { - // the EOS does not exhibit any physically meaningful - // extrema, and the fluid is critical... - Vm = VmCubic; - handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase); - } + // Evaluation Vmin, Vmax, pmin, pmax; + // if (findExtrema_(Vmin, Vmax, + // pmin, pmax, + // a, b, T)) + // { + // if (isGasPhase) + // Vm = std::max(Vmax, VmCubic); + // else { + // if (Vmin > 0) + // Vm = std::min(Vmin, VmCubic); + // else + // Vm = VmCubic; + // } + // } + // else { + // // the EOS does not exhibit any physically meaningful + // // extrema, and the fluid is critical... + // Vm = VmCubic; + // handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase); + // } } Valgrind::CheckDefined(Vm); diff --git a/opm/material/eos/PengRobinsonMixture.hpp b/opm/material/eos/PengRobinsonMixture.hpp index ce987177f..d93453999 100644 --- a/opm/material/eos/PengRobinsonMixture.hpp +++ b/opm/material/eos/PengRobinsonMixture.hpp @@ -87,12 +87,28 @@ public: * R. Reid, et al.: The Properties of Gases and Liquids, * 4th edition, McGraw-Hill, 1987, pp. 42-44, 143-145 */ +#warning should check why this function is changed template - static LhsEval computeFugacityCoefficient(const FluidState& fs, - const Params& params, + static LhsEval computeFugacityCoefficient(const FluidState& fs2, + const Params& params2, unsigned phaseIdx, unsigned compIdx) { +#warning also a HACK, should investigate why + auto fs = fs2; + double sumx = 0.0; + for (int i = 0; i < FluidState::numComponents; ++i) + sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i)); + if (sumx < 0.95) { + double alpha = 0.95/sumx; + std::cerr << "normalize: " << sumx + << " alpha: " << alpha << "\n"; + for (int i = 0; i < FluidState::numComponents; ++i) + fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i)); + } + + auto params = params2; + params.updatePhase(fs, phaseIdx); // note that we normalize the component mole fractions, so // that their sum is 100%. This increases numerical stability // considerably if the fluid state is not physical. diff --git a/tests/test_chiflash.cpp b/tests/test_chiflash.cpp index 4d7823727..f81f2e0a3 100644 --- a/tests/test_chiflash.cpp +++ b/tests/test_chiflash.cpp @@ -78,7 +78,18 @@ void testChiFlash() fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx)); ComponentVector zInit(0.); // TODO; zInit needs to be normalized. - const double flash_tolerance = 1.e-8; + { + Scalar sumMoles = 0.0; + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + Scalar tmp = Opm::getValue(fs.molarity(phaseIdx, compIdx) * fs.saturation(phaseIdx)); + zInit[compIdx] += Opm::max(tmp, 1e-8); + sumMoles += tmp; + } + } + zInit /= sumMoles; + } + const double flash_tolerance = -1.; // just to test the setup in co2-compositional const int flash_verbosity = 1; const std::string flash_twophase_method = "ssi";