compiled and running now.

This commit is contained in:
Kai Bao 2021-11-12 13:07:01 +01:00 committed by Trine Mykkeltvedt
parent c6064aa3e9
commit e7da7cdad6
3 changed files with 56 additions and 25 deletions

View File

@ -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);

View File

@ -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 <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
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.

View File

@ -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";