diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index 9be32dacf5..7beb6d616f 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -339,12 +339,17 @@ typename FlowGenericProblem::Scalar FlowGenericProblem:: rockFraction(unsigned elementIdx, unsigned timeIdx) const { - // the reference porosity is defined as the accumulated pore volume divided by the - // geometric volume of the element. Note that it can - // be larger than 1.0 if porevolume multipliers are used - // to for instance implement larger boundary cells - auto porosity = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(), "PORO", elementIdx); - return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity); + // For the energy equation, we need the volume of the rock. + // The volume of the rock is computed by rockFraction * geometric volume of the element. + // The reference porosity is defined as porosity * ntg * pore-volume-multiplier. + // A common practice in reservoir simulation is to use large pore-volume-multipliers in boundary cells + // to model boundary conditions other than no-flow. This may result in reference porosities that are larger than 1. + // A simple (1-reference porosity) * geometric volume of the element may give unphysical results. + // We therefore instead consider the pore-volume-multiplier as a volume multiplier. The rock fraction is thus given by + // (1 - porosity * ntg) * pore-volume-multiplier = (1 - porosity * ntg) * reference porosity / (porosity * ntg) + const auto ntg = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(), "NTG", elementIdx); + const auto poro_eff = ntg * this->lookUpData_.fieldPropDouble(eclState_.fieldProps(), "PORO", elementIdx); + return (1 - poro_eff) * referencePorosity(elementIdx, timeIdx) / poro_eff; } template