Merge pull request #36 from andlaus/inconsistent_hysteresis_updates

EclDefaultMaterial: use inconsistent saturations for the hysteresis update
This commit is contained in:
Andreas Lauser 2015-08-14 15:32:03 +02:00
commit 8e51155f24
2 changed files with 38 additions and 12 deletions

View File

@ -181,10 +181,7 @@ public:
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result);
return result;
return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
}
/*!
@ -327,18 +324,35 @@ public:
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Swco = params.Swl();
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx));
Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
Sw = std::min(1.0, std::max(Swco, Sw));
Sg = std::min(1.0, std::max(0.0, Sg));
Scalar Sw_ow = Sg + Sw;
Scalar So_go = 1 - Sw_ow + Swco;
if (params.inconsistentHysteresisUpdate()) {
Sg = std::min(1.0, std::max(0.0, Sg));
// NOTE: the saturations which are passed to update the hysteresis curves are
// inconsistent with the ones used to calculate the relative permabilities. We do
// it like this anyway because (a) the saturation functions of opm-core do it
// this way (b) the simulations seem to converge better (which is not too much
// surprising actually, because the time step does not start on a kink in the
// solution) and (c) the Eclipse 100 simulator may do the same.
//
// Though be aware that from a physical perspective this is definitively
// incorrect!
params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
}
else {
Scalar Swco = params.Swl();
Sw = std::min(1.0, std::max(0.0, Sw));
Sg = std::min(1.0, std::max(0.0, Sg));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
Scalar Sw_ow = Sg + std::max(Swco, Sw);
Scalar So_go = 1 + Swco - Sw_ow;
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
}
}
};
} // namespace Opm

View File

@ -124,6 +124,18 @@ public:
Scalar Swl() const
{ assertFinalized_(); return Swl_; }
/*!
* \brief Specify whether inconsistent saturations should be used to update the
* hysteresis parameters.
*
* Returning 'true' is wrong from a physical point of view because the saturations
* which are used to update the hysteresis parameters are calculated differently than
* the ones used to calculate the relperms and capillary pressures. Since Eclipse
* E100 probably uses inconsistent saturations, we return true here anyway.
*/
bool inconsistentHysteresisUpdate() const
{ return true; }
private:
#ifndef NDEBUG
void assertFinalized_() const