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>
 
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 const char* phaseName(unsigned phaseIdx)
    {
        static const char* 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 const char* componentName(unsigned compIdx)
    {
        static const char* name[] = {
            H2O::name(),
            "C1",
            "C3",
            "C6",
            "C10",
            "C15",
            "C20"
        };
 
        assert(0 <= compIdx && 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::Scalar, 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::Scalar, 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::Scalar, 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