Spe5FluidSystem.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_SPE5_FLUID_SYSTEM_HPP
#define OPM_SPE5_FLUID_SYSTEM_HPP

#include "BaseFluidSystem.hpp"
#include "Spe5ParameterCache.hpp"

#include <opm/material/Constants.hpp>
#include <opm/material/eos/PengRobinsonMixture.hpp>

#include <opm/material/common/Spline.hpp>

#include <string_view>

namespace Opm {

template <class Scalar>
class Spe5FluidSystem
    : public BaseFluidSystem<Scalar, Spe5FluidSystem<Scalar> >
{
    typedef Spe5FluidSystem<Scalar> ThisType;

    typedef typename ::Opm::PengRobinsonMixture<Scalar, ThisType> PengRobinsonMixture;
    typedef typename ::Opm::PengRobinson<Scalar> PengRobinson;

    static const Scalar R;

public:
    template <class Evaluation>
    struct ParameterCache : public Spe5ParameterCache<Evaluation, ThisType>
    {};

    /****************************************
     * Fluid phase parameters
     ****************************************/

    static const int numPhases = 3;

    static const int gasPhaseIdx = 0;
    static const int waterPhaseIdx = 1;
    static const int oilPhaseIdx = 2;

    typedef ::Opm::H2O<Scalar> H2O;

    static std::string_view phaseName(unsigned phaseIdx)
    {
        static const std::string_view name[] = {
            "gas",
            "water",
            "oil",
        };

        assert(phaseIdx < numPhases);
        return name[phaseIdx];
    }

    static bool isLiquid(unsigned phaseIdx)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx != gasPhaseIdx;
    }

    static bool isCompressible(unsigned /*phaseIdx*/)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }

    static bool isIdealGas(unsigned /*phaseIdx*/)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        return false; // gas is not ideal here!
    }

    static bool isIdealMixture(unsigned phaseIdx)
    {
        // always use the reference oil for the fugacity coefficents,
        // so they cannot be dependent on composition and they the
        // phases thus always an ideal mixture
        return phaseIdx == waterPhaseIdx;
    }

    /****************************************
     * Component related parameters
     ****************************************/

    static const int numComponents = 7;

    static const int H2OIdx = 0; 
    static const int C1Idx = 1; 
    static const int C3Idx = 2; 
    static const int C6Idx = 3; 
    static const int C10Idx = 4; 
    static const int C15Idx = 5; 
    static const int C20Idx = 6; 

    static std::string_view componentName(unsigned compIdx)
    {
        static const std::string_view name[] = {
            H2O::name(),
            "C1",
            "C3",
            "C6",
            "C10",
            "C15",
            "C20"
        };

        assert(compIdx < numComponents);
        return name[compIdx];
    }

    static Scalar molarMass(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::molarMass()
            : (compIdx == C1Idx)
            ? 16.04e-3
            : (compIdx == C3Idx)
            ? 44.10e-3
            : (compIdx == C6Idx)
            ? 86.18e-3
            : (compIdx == C10Idx)
            ? 142.29e-3
            : (compIdx == C15Idx)
            ? 206.00e-3
            : (compIdx == C20Idx)
            ? 282.00e-3
            : 1e30;
    }

    static Scalar criticalTemperature(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::criticalTemperature()
            : (compIdx == C1Idx)
            ? 343.0*5/9
            : (compIdx == C3Idx)
            ? 665.7*5/9
            : (compIdx == C6Idx)
            ? 913.4*5/9
            : (compIdx == C10Idx)
            ? 1111.8*5/9
            : (compIdx == C15Idx)
            ? 1270.0*5/9
            : (compIdx == C20Idx)
            ? 1380.0*5/9
            : 1e30;
    }

    static Scalar criticalPressure(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::criticalPressure()
            : (compIdx == C1Idx)
            ? 667.8 * 6894.7573
            : (compIdx == C3Idx)
            ? 616.3 * 6894.7573
            : (compIdx == C6Idx)
            ? 436.9 * 6894.7573
            : (compIdx == C10Idx)
            ? 304.0 * 6894.7573
            : (compIdx == C15Idx)
            ? 200.0 * 6894.7573
            : (compIdx == C20Idx)
            ? 162.0 * 6894.7573
            : 1e30;
    }

    static Scalar criticalMolarVolume(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::criticalMolarVolume()
            : (compIdx == C1Idx)
            ? 0.290*R*criticalTemperature(C1Idx)/criticalPressure(C1Idx)
            : (compIdx == C3Idx)
            ? 0.277*R*criticalTemperature(C3Idx)/criticalPressure(C3Idx)
            : (compIdx == C6Idx)
            ? 0.264*R*criticalTemperature(C6Idx)/criticalPressure(C6Idx)
            : (compIdx == C10Idx)
            ? 0.257*R*criticalTemperature(C10Idx)/criticalPressure(C10Idx)
            : (compIdx == C15Idx)
            ? 0.245*R*criticalTemperature(C15Idx)/criticalPressure(C15Idx)
            : (compIdx == C20Idx)
            ? 0.235*R*criticalTemperature(C20Idx)/criticalPressure(C20Idx)
            : 1e30;
    }

    static Scalar acentricFactor(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::acentricFactor()
            : (compIdx == C1Idx)
            ? 0.0130
            : (compIdx == C3Idx)
            ? 0.1524
            : (compIdx == C6Idx)
            ? 0.3007
            : (compIdx == C10Idx)
            ? 0.4885
            : (compIdx == C15Idx)
            ? 0.6500
            : (compIdx == C20Idx)
            ? 0.8500
            : 1e30;
    }

    static Scalar interactionCoefficient(unsigned comp1Idx, unsigned comp2Idx)
    {
        unsigned i = std::min(comp1Idx, comp2Idx);
        unsigned j = std::max(comp1Idx, comp2Idx);
        if (i == C1Idx && (j == C15Idx || j == C20Idx))
            return 0.05;
        else if (i == C3Idx && (j == C15Idx || j == C20Idx))
            return 0.005;
        return 0;
    }

    /****************************************
     * Methods which compute stuff
     ****************************************/

    static void init(Scalar minT = 273.15,
                     Scalar maxT = 373.15,
                     Scalar minP = 1e4,
                     Scalar maxP = 100e6)
    {
        PengRobinsonParamsMixture<Scalar, ThisType, gasPhaseIdx, /*useSpe5=*/true> prParams;

        // find envelopes of the 'a' and 'b' parameters for the range
        // minT <= T <= maxT and minP <= p <= maxP. For
        // this we take advantage of the fact that 'a' and 'b' for
        // mixtures is just a convex combination of the attractive and
        // repulsive parameters of the pure components

        Scalar minA = 1e30, maxA = -1e30;
        Scalar minB = 1e30, maxB = -1e30;

        prParams.updatePure(minT, minP);
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = std::min(prParams.pureParams(compIdx).a(), minA);
            maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
            minB = std::min(prParams.pureParams(compIdx).b(), minB);
            maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
        };

        prParams.updatePure(maxT, minP);
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = std::min(prParams.pureParams(compIdx).a(), minA);
            maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
            minB = std::min(prParams.pureParams(compIdx).b(), minB);
            maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
        };

        prParams.updatePure(minT, maxP);
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = std::min(prParams.pureParams(compIdx).a(), minA);
            maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
            minB = std::min(prParams.pureParams(compIdx).b(), minB);
            maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
        };

        prParams.updatePure(maxT, maxP);
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = std::min(prParams.pureParams(compIdx).a(), minA);
            maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
            minB = std::min(prParams.pureParams(compIdx).b(), minB);
            maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
        };

        PengRobinson::init(/*aMin=*/minA, /*aMax=*/maxA, /*na=*/100,
                           /*bMin=*/minB, /*bMax=*/maxB, /*nb=*/200);
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval density(const FluidState& fluidState,
                           const ParameterCache<ParamCacheEval>& paramCache,
                           unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval viscosity(const FluidState& /*fluidState*/,
                             const ParameterCache<ParamCacheEval>& /*paramCache*/,
                             unsigned phaseIdx)
    {
        assert(phaseIdx <= numPhases);

        if (phaseIdx == gasPhaseIdx) {
            // given by SPE-5 in table on page 64. we use a constant
            // viscosity, though...
            return 0.0170e-2 * 0.1;
        }
        else if (phaseIdx == waterPhaseIdx)
            // given by SPE-5: 0.7 centi-Poise  = 0.0007 Pa s
            return 0.7e-2 * 0.1;
        else {
            assert(phaseIdx == oilPhaseIdx);
            // given by SPE-5 in table on page 64. we use a constant
            // viscosity, though...
            return 0.208e-2 * 0.1;
        }
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval fugacityCoefficient(const FluidState& fluidState,
                                       const ParameterCache<ParamCacheEval>& paramCache,
                                       unsigned phaseIdx,
                                       unsigned compIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(compIdx <= numComponents);

        if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx)
            return PengRobinsonMixture::computeFugacityCoefficient(fluidState,
                                                                   paramCache,
                                                                   phaseIdx,
                                                                   compIdx);
        else {
            assert(phaseIdx == waterPhaseIdx);
            return
                henryCoeffWater_(compIdx, fluidState.temperature(waterPhaseIdx))
                / fluidState.pressure(waterPhaseIdx);
        }
    }

protected:
    template <class LhsEval>
    static LhsEval henryCoeffWater_(unsigned compIdx, const LhsEval& temperature)
    {
        // use henry's law for the solutes and the vapor pressure for
        // the solvent.
        switch (compIdx) {
        case H2OIdx: return H2O::vaporPressure(temperature);

            // the values of the Henry constant for the solutes have
            // are faked so far...
        case C1Idx: return 5.57601e+09;
        case C3Idx: return 1e10;
        case C6Idx: return 1e10;
        case C10Idx: return 1e10;
        case C15Idx: return 1e10;
        case C20Idx: return 1e10;
        default: throw std::logic_error("Unknown component index "+std::to_string(compIdx));
        }
    }
};

template <class Scalar>
const Scalar Spe5FluidSystem<Scalar>::R = Constants<Scalar>::R;

} // namespace Opm

#endif