Apply the SSHIFT volume shift to the density and the viscosity

CMakeLists_files.cmake

@@ -717,6 +717,7 @@ list(APPEND DUNE_TEST_SOURCE_FILES
   tests/material/test_ptflash_ssi_newton_fallback.cpp
   tests/material/test_tabulation.cpp
   tests/material/test_threecomponents_ptflash.cpp
+  tests/material/test_volume_shift.cpp
 )
 
 if(dune-common_FOUND)

opm/material/checkFluidSystem.hpp

@@ -291,14 +291,18 @@ ParameterCache initParamCache()
 /*!
  * \brief Checks whether a fluid system adheres to the specification.
  */
-template <class Scalar, class FluidSystem, class RhsEval, class LhsEval>
-void checkFluidSystem()
+template <class Scalar, class FluidSystem, class RhsEval, class LhsEval, class Initializer>
+void checkFluidSystem(Initializer initializeFluidSystem)
 {
     std::cout << "Testing fluid system '"
               << Opm::getDemangledType<FluidSystem>()
               << ", RhsEval = " << Opm::getDemangledType<RhsEval>()
               << ", LhsEval = " << Opm::getDemangledType<LhsEval>() << "'\n";
 
+    // Initialize before creating the fluid state or its cache. Configurable
+    // systems need to register their components after init() clears them.
+    initializeFluidSystem();
+
     // make sure the fluid system provides the number of phases and
     // the number of components
     static constexpr int numPhases = FluidSystem::numPhases;
@@ -354,7 +358,6 @@ void checkFluidSystem()
     val = 2*val; // get rid of GCC warning (only occurs with paranoid warning flags)
 
     // actually check the fluid system API
-    try { FluidSystem::init(); } catch (...) {};
     for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
         fs.restrictToPhase(static_cast<int>(phaseIdx));
         fs.allowPressure(FluidSystem::isCompressible(phaseIdx));
@@ -405,4 +408,12 @@ void checkFluidSystem()
     }
 }
 
+template <class Scalar, class FluidSystem, class RhsEval, class LhsEval>
+void checkFluidSystem()
+{
+    checkFluidSystem<Scalar, FluidSystem, RhsEval, LhsEval>([] {
+        try { FluidSystem::init(); } catch (...) {};
+    });
+}
+
 #endif

opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp

@@ -93,15 +93,18 @@ namespace Opm {
             Scalar critic_pres; // unit: parscal
             Scalar critic_vol; // unit: m^3/kmol
             Scalar acentric_factor; // unit: dimension less
+            Scalar volume_shift; // dimensionless SSHIFT coefficient
 
             ComponentParam(const std::string_view name_, const Scalar molar_mass_, const Scalar critic_temp_,
-                           const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_)
+                           const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_,
+                           const Scalar volume_shift_ = 0.0)
                     : name(name_),
                       molar_mass(molar_mass_),
                       critic_temp(critic_temp_),
                       critic_pres(critic_pres_),
                       critic_vol(critic_vol_),
-                      acentric_factor(acentric_factor_)
+                      acentric_factor(acentric_factor_),
+                      volume_shift(volume_shift_)
             {}
         };
 
@@ -177,7 +180,10 @@ namespace Opm {
                                                    static_cast<Scalar>(eos_props.critical_temperature[c]),
                                                    static_cast<Scalar>(eos_props.critical_pressure[c]),
                                                    static_cast<Scalar>(eos_props.critical_volume[c] * 1.e3),
-                                                   static_cast<Scalar>(eos_props.acentric_factors[c])});
+                                                   static_cast<Scalar>(eos_props.acentric_factors[c]),
+                                                   c < eos_props.volume_shifts.size()
+                                                       ? static_cast<Scalar>(eos_props.volume_shifts[c])
+                                                       : Scalar{0}});
             }
 
             const auto& bic = eos_props.binary_interaction_coefficient;
@@ -199,7 +205,12 @@ namespace Opm {
         static void init()
         {
             waterPvt_ = std::make_shared<WaterPvt>();
+            // Discard the previous configuration so subsequent component
+            // registrations replace it.
+            component_param_.clear();
             component_param_.reserve(numComponents);
+            interaction_coefficients_.clear();
+            lbc_coefficients_ = ViscosityModel::defaultLBCCoefficients();
         }
 
         /*!
@@ -220,6 +231,23 @@ namespace Opm {
 
             return component_param_[compIdx].acentric_factor;
         }
+
+        /*!
+         * \brief Dimensionless volume-translation coefficient (SSHIFT).
+         *
+         * The component volume correction is s_c b_c, where b_c is the
+         * equation-of-state covolume.
+         *
+         * \copydetails Doxygen::compIdxParam
+         */
+        static Scalar volumeShift(unsigned compIdx)
+        {
+            assert(isConsistent());
+            assert(compIdx < numComponents);
+
+            return component_param_[compIdx].volume_shift;
+        }
+
         /*!
          * \brief Critical temperature of a component [K].
          *
@@ -325,7 +353,10 @@ namespace Opm {
             assert(phaseIdx < numPhases);
 
             if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
-                return decay<LhsEval>(fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx));
+                // Density uses the translated volume; fugacity coefficients
+                // use the unshifted equation-of-state root.
+                const auto Vm = paramCache.correctedMolarVolume(phaseIdx);
+                return decay<LhsEval>(fluidState.averageMolarMass(phaseIdx) / Vm);
             }
             else {
                 const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
@@ -348,8 +379,12 @@ namespace Opm {
             assert(phaseIdx < numPhases);
 
             if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
-                // Use LBC method to calculate viscosity
-                return decay<LhsEval>(ViscosityModel::LBC(fluidState, paramCache, phaseIdx));
+                // LBC is a reduced-density correlation, so use the physical
+                // molar density after applying SSHIFT.
+                const auto rho = density(fluidState, paramCache, phaseIdx);
+                const auto molarDensity = rho / fluidState.averageMolarMass(phaseIdx);
+                return decay<LhsEval>(
+                    ViscosityModel::LBCWithMolarDensity(fluidState, molarDensity, phaseIdx));
             }
             else {
                 const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
@@ -358,7 +393,15 @@ namespace Opm {
             }
         }
 
-        //! \copydoc BaseFluidSystem::fugacityCoefficient
+        /*!
+         * \copydoc BaseFluidSystem::fugacityCoefficient
+         *
+         * Returns the coefficient of the unshifted cubic EOS. The PT flash
+         * uses these coefficients at equal phase pressure and temperature,
+         * where the SSHIFT translation factors cancel from their ratios.
+         * An absolute coefficient of the translated EOS additionally requires
+         * the factor exp(-p s_c b_c / (R T)), which is not included here.
+         */
         template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
         static LhsEval fugacityCoefficient(const FluidState& fluidState,
                                            const ParameterCache<ParamCacheEval>& paramCache,

opm/material/fluidsystems/PTFlashParameterCache.hpp

@@ -30,13 +30,17 @@
 #ifndef OPM_PTFlash_PARAMETER_CACHE_HPP
 #define OPM_PTFlash_PARAMETER_CACHE_HPP
 
+#include <opm/common/Exceptions.hpp>
+
 #include <opm/material/common/Valgrind.hpp>
 #include <opm/material/fluidsystems/ParameterCacheBase.hpp>
 #include <opm/material/eos/CubicEOS.hpp>
 #include <opm/material/eos/CubicEOSParams.hpp>
 
 #include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
 
+#include <fmt/format.h>
+
 #include <cassert>
 
 namespace Opm {
@@ -80,8 +84,10 @@ class PTFlashParameterCache
     {
         VmUpToDate_[oilPhaseIdx] = false;
         Valgrind::SetUndefined(Vm_[oilPhaseIdx]);
+        Valgrind::SetUndefined(volumeShift_[oilPhaseIdx]);
         VmUpToDate_[gasPhaseIdx] = false;
         Valgrind::SetUndefined(Vm_[gasPhaseIdx]);
+        Valgrind::SetUndefined(volumeShift_[gasPhaseIdx]);
 
         oilPhaseParams_.setEOSType(eos_type);
         gasPhaseParams_.setEOSType(eos_type);
@@ -286,6 +292,44 @@ class PTFlashParameterCache
         return Vm_[phaseIdx];
     }
 
+    /*!
+     * \brief Phase volume translation, sum_c x_c s_c b_c [m^3/mol].
+     *
+     * Call updatePhase() before reading this cached value. The translation
+     * is stored separately from the equation-of-state root in molarVolume().
+     * correctedMolarVolume() subtracts it to obtain the physical molar volume.
+     *
+     * \param phaseIdx The fluid phase of interest
+     */
+    Scalar phaseVolumeShift(unsigned phaseIdx) const
+    {
+        return volumeShift_[phaseIdx];
+    }
+
+    /*!
+     * \brief Molar volume after applying SSHIFT [m^3/mol].
+     *
+     * Density, viscosity and phase saturations use this physical volume.
+     * Fugacity coefficients and equation-of-state compressibility factors
+     * use the unshifted molarVolume(). At equal phase pressure and temperature,
+     * the volume-translation factors cancel from the equilibrium ratios.
+     *
+     * \param phaseIdx The fluid phase of interest
+     */
+    Scalar correctedMolarVolume(unsigned phaseIdx) const
+    {
+        const Scalar Vm = molarVolume(phaseIdx) - phaseVolumeShift(phaseIdx);
+
+        // Reject non-positive or NaN volumes before computing fluid properties.
+        if (!(scalarValue(Vm) > 0)) {
+            throw NumericalProblem(
+                fmt::format("The SSHIFT volume shift of phase {} leaves a corrected "
+                            "molar volume of {}, which is not positive.",
+                            phaseIdx, scalarValue(Vm)));
+        }
+        return Vm;
+    }
+
 
     /*!
      * \brief Returns the Peng-Robinson mixture parameters for the oil
@@ -381,6 +425,7 @@ class PTFlashParameterCache
                             unsigned phaseIdx)
     {
         VmUpToDate_[phaseIdx] = true;
+        volumeShift_[phaseIdx] = computeVolumeShift_(fluidState, phaseIdx);
 
         // calculate molar volume of the phase (we will need this for the
         // fugacity coefficients and the density anyway)
@@ -415,8 +460,30 @@ class PTFlashParameterCache
         };
     }
 
+    //! \brief Compute the phase volume translation, sum_c x_c s_c b_c [m^3/mol].
+    template <class FluidState>
+    Scalar computeVolumeShift_(const FluidState& fluidState, unsigned phaseIdx) const
+    {
+        Scalar shift = 0;
+        // Fluid systems without SSHIFT support use a zero translation.
+        if constexpr (requires { FluidSystem::volumeShift(0u); }) {
+            // b_c from the dimensionless B_c = b_c p / (R T).
+            const Scalar T = decay<Scalar>(fluidState.temperature(phaseIdx));
+            const Scalar p = decay<Scalar>(fluidState.pressure(phaseIdx));
+            const Scalar RT_p = Constants<Scalar>::R * T / p;
+
+            for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
+                const Scalar b = decay<Scalar>(Bi(phaseIdx, compIdx)) * RT_p;
+                shift += decay<Scalar>(fluidState.moleFraction(phaseIdx, compIdx))
+                       * FluidSystem::volumeShift(compIdx) * b;
+            }
+        }
+        return shift;
+    }
+
     bool VmUpToDate_[numPhases];
     Scalar Vm_[numPhases];
+    Scalar volumeShift_[numPhases];
 
     OilPhaseParams oilPhaseParams_;
     GasPhaseParams gasPhaseParams_;

opm/material/viscositymodels/ViscosityModels.hpp

@@ -62,11 +62,24 @@ class ViscosityModels
                       const Params& /*paramCache*/,
                       unsigned phaseIdx)
     {
-        const Scalar MPa_atm = 0.101325;
         const Scalar R = Opm::Constants<Scalar>::R;
         const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
         const auto& P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
         const auto& Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
+        const LhsEval molarDensity = P / (R * T * Z);
+
+        return LBCWithMolarDensity<FluidState, LhsEval, LhsEval>(fluidState, molarDensity, phaseIdx);
+    }
+
+    // LBC correlation at the supplied physical molar density [mol/m^3].
+    template <class FluidState, class MolarDensity,
+              class LhsEval = typename FluidState::ValueType>
+    static LhsEval LBCWithMolarDensity(const FluidState& fluidState,
+                                       const MolarDensity& molarDensity,
+                                       unsigned phaseIdx)
+    {
+        const Scalar MPa_atm = 0.101325;
+        const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
 
         LhsEval sumVolume = 0.0;
         for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
@@ -76,9 +89,7 @@ class ViscosityModels
         }
 
         LhsEval rho_pc = 1.0 / sumVolume;
-        LhsEval V = (R * T * Z)/P;
-        LhsEval rho = 1.0 / V;
-        LhsEval rho_r = rho / rho_pc;
+        const LhsEval rho_r = Opm::decay<LhsEval>(molarDensity) / rho_pc;
 
         LhsEval xsum_T_c = 0.0; // mixture pseudocritical temperature
         LhsEval xsum_Mm = 0.0; // mixture molar mass

tests/material/test_fluidsystems.cpp

@@ -393,31 +393,32 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GenericFluidSystem, Scalar, ScalarTypes)
     using Evaluation = Opm::DenseAd::Evaluation<Scalar, 4>;
     using FluidSystem = Opm::GenericOilGasWaterFluidSystem<Scalar, 4, true>;
 
-    registerComponent<FluidSystem, Opm::SimpleCO2<Scalar>>();
-    registerComponent<FluidSystem, Opm::C1<Scalar>>();
-    registerComponent<FluidSystem, Opm::C10<Scalar>>();
-    registerComponent<FluidSystem, Opm::N2<Scalar>>();
-
-    // initialize water pvt
-    using WaterPvt = typename FluidSystem::WaterPvt;
-    std::shared_ptr<WaterPvt> waterPvt;
-    FluidSystem::setWaterPvt(waterPvt);
-
-    checkFluidSystem<Scalar, FluidSystem, Scalar, Scalar>();
-    checkFluidSystem<Scalar, FluidSystem, Evaluation, Scalar>();
-    checkFluidSystem<Scalar, FluidSystem, Evaluation, Evaluation>();
+    const auto initializeFluidSystem = [] {
+        FluidSystem::init();
+        registerComponent<FluidSystem, Opm::SimpleCO2<Scalar>>();
+        registerComponent<FluidSystem, Opm::C1<Scalar>>();
+        registerComponent<FluidSystem, Opm::C10<Scalar>>();
+        registerComponent<FluidSystem, Opm::N2<Scalar>>();
+    };
+
+    checkFluidSystem<Scalar, FluidSystem, Scalar, Scalar>(initializeFluidSystem);
+    checkFluidSystem<Scalar, FluidSystem, Evaluation, Scalar>(initializeFluidSystem);
+    checkFluidSystem<Scalar, FluidSystem, Evaluation, Evaluation>(initializeFluidSystem);
 }
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(GenericFluidSystemNoWater, Scalar, ScalarTypes)
 {
     using Evaluation = Opm::DenseAd::Evaluation<Scalar, 3>;
     using FluidSystem = Opm::GenericOilGasWaterFluidSystem<Scalar, 3, false>;
 
-    registerComponent<FluidSystem, Opm::SimpleCO2<Scalar>>();
-    registerComponent<FluidSystem, Opm::C1<Scalar>>();
-    registerComponent<FluidSystem, Opm::C10<Scalar>>();
+    const auto initializeFluidSystem = [] {
+        FluidSystem::init();
+        registerComponent<FluidSystem, Opm::SimpleCO2<Scalar>>();
+        registerComponent<FluidSystem, Opm::C1<Scalar>>();
+        registerComponent<FluidSystem, Opm::C10<Scalar>>();
+    };
 
-    checkFluidSystem<Scalar, FluidSystem, Scalar, Scalar>();
-    checkFluidSystem<Scalar, FluidSystem, Evaluation, Scalar>();
-    checkFluidSystem<Scalar, FluidSystem, Evaluation, Evaluation>();
+    checkFluidSystem<Scalar, FluidSystem, Scalar, Scalar>(initializeFluidSystem);
+    checkFluidSystem<Scalar, FluidSystem, Evaluation, Scalar>(initializeFluidSystem);
+    checkFluidSystem<Scalar, FluidSystem, Evaluation, Evaluation>(initializeFluidSystem);
 }