From bf283f86849211f5e858d8e1e6bc312a269dea18 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 8 Mar 2016 14:52:09 +0100 Subject: [PATCH] threshold pressures: consider the saturated case and fix a typo the typo was caused the surface density of the oil phase to be used instead of the one of gas. This caused the density to be off by a factor of typically about 900. using saturated FVFs does not change much, but it does not hurt because it is also done that way in the simulator. This makes the defaults for the threshold pressures reasonable again, but for some reason they are not exactly the same as in the old implementation. (although the differences are very tolerable.) On the question why only "Model 2" is affected by this: the other decks don't use threshold pressures (SPE-X) or do not default any values (Norne). --- opm/simulators/thresholdPressures.hpp | 37 +++++++++++++++++++-------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/opm/simulators/thresholdPressures.hpp b/opm/simulators/thresholdPressures.hpp index 31a916580..30c215351 100644 --- a/opm/simulators/thresholdPressures.hpp +++ b/opm/simulators/thresholdPressures.hpp @@ -198,18 +198,25 @@ void computeMaxDp(std::map, double>& maxDp, double T = initialState.temperature()[cellIdx]; double p = initialState.pressure()[cellIdx]; double Rs = initialState.gasoilratio()[cellIdx]; - double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p); - rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; - if (pu.phase_used[BlackoilPhases::Vapour]) { - int gpos = pu.phase_pos[BlackoilPhases::Vapour]; - rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b; + if (Rs >= RsSat) { + double b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); + rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; + } + else { + double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b; + } } } } if (pu.phase_used[BlackoilPhases::Vapour]) { - const int gpos = pu.phase_pos[BlackoilPhases::Liquid]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; const auto& pvtg = props.gasPvt(); for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { int pvtRegionIdx = pvtRegion[cellIdx]; @@ -217,12 +224,20 @@ void computeMaxDp(std::map, double>& maxDp, double T = initialState.temperature()[cellIdx]; double p = initialState.pressure()[cellIdx]; double Rv = initialState.rv()[cellIdx]; - double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); + double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p); - rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; - if (pu.phase_used[BlackoilPhases::Liquid]) { - int opos = pu.phase_pos[BlackoilPhases::Liquid]; - rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b; + if (Rv >= RvSat) { + double b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); + rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; + } + else { + double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); + rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b; + } } } }