spe5parametercache.hh

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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef SPE5_PARAMETER_CACHE_HH
#define SPE5_PARAMETER_CACHE_HH
 
#include <cassert>
 
#include <dumux/material/components/h2o.hh>
#include <dumux/material/fluidsystems/parametercachebase.hh>
 
#include <dumux/material/eos/pengrobinson.hh>
#include <dumux/material/eos/pengrobinsonparamsmixture.hh>
 
namespace Dumux {

template <class Scalar, class FluidSystem>
class Spe5ParameterCache
    : public ParameterCacheBase<Spe5ParameterCache<Scalar, FluidSystem> >
{
    using ThisType = Spe5ParameterCache<Scalar, FluidSystem>;
    using ParentType = ParameterCacheBase<ThisType>;
 
    using PengRobinson = Dumux::PengRobinson<Scalar>;
 
    enum { numPhases = FluidSystem::numPhases };
 
    enum { wPhaseIdx = FluidSystem::wPhaseIdx };
    enum { oPhaseIdx = FluidSystem::oPhaseIdx };
    enum { gPhaseIdx = FluidSystem::gPhaseIdx };
 
public:
    // types of the parameter objects for each phase
    using OilPhaseParams = PengRobinsonParamsMixture<Scalar, FluidSystem, oPhaseIdx, /*useSpe5=*/true>;
    using GasPhaseParams = PengRobinsonParamsMixture<Scalar, FluidSystem, gPhaseIdx, /*useSpe5=*/true>;

    Spe5ParameterCache()
    {
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            VmUpToDate_[phaseIdx] = false;
            Valgrind::SetUndefined(Vm_[phaseIdx]);
        }
    }

    template <class FluidState>
    void updatePhase(const FluidState &fs,
                     int phaseIdx,
                     int except = ParentType::None)
    {
        updateEosParams(fs, phaseIdx, except);
 
        // if we don't need to recalculate the molar volume, we exit
        // here
        if (VmUpToDate_[phaseIdx])
            return;
 
        // update the phase's molar volume
        updateMolarVolume_(fs, phaseIdx);
    }

    template <class FluidState>
    void updateSingleMoleFraction(const FluidState &fs,
                                  int phaseIdx,
                                  int compIdx)
    {
        if (phaseIdx == oPhaseIdx)
            oilPhaseParams_.updateSingleMoleFraction(fs, compIdx);
        else if (phaseIdx == gPhaseIdx)
            gasPhaseParams_.updateSingleMoleFraction(fs, compIdx);
 
        // update the phase's molar volume
        updateMolarVolume_(fs, phaseIdx);
    }

    Scalar a(int phaseIdx) const
    {
        switch (phaseIdx)
        {
            case oPhaseIdx: return oilPhaseParams_.a();
            case gPhaseIdx: return gasPhaseParams_.a();
            default:
                DUNE_THROW(Dune::InvalidStateException,
                           "The a() parameter is only defined for "
                           "oil and gas phases");
        }
    }

    Scalar b(int phaseIdx) const
    {
        switch (phaseIdx)
        {
            case oPhaseIdx: return oilPhaseParams_.b();
            case gPhaseIdx: return gasPhaseParams_.b();
            default:
                DUNE_THROW(Dune::InvalidStateException,
                           "The b() parameter is only defined for "
                           "oil and gas phases");
        }
    }

    Scalar aPure(int phaseIdx, int compIdx) const
    {
        switch (phaseIdx)
        {
            case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a();
            case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
            default:
                DUNE_THROW(Dune::InvalidStateException,
                           "The a() parameter is only defined for "
                           "oil and gas phases");
        }
    }

    Scalar bPure(int phaseIdx, int compIdx) const
    {
        switch (phaseIdx)
        {
            case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b();
            case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
            default:
                DUNE_THROW(Dune::InvalidStateException,
                           "The b() parameter is only defined for "
                           "oil and gas phases");
        }
    }

    Scalar molarVolume(int phaseIdx) const
    {
        assert(VmUpToDate_[phaseIdx]);
        return Vm_[phaseIdx];
    }
 

    const OilPhaseParams &oilPhaseParams() const
    { return oilPhaseParams_; }

    const GasPhaseParams &gasPhaseParams() const
    { return gasPhaseParams_; }

    template <class FluidState>
    void updateEosParams(const FluidState &fs,
                         int phaseIdx,
                         int exceptQuantities = ParentType::None)
    {
        if (!(exceptQuantities & ParentType::Temperature))
        {
            updatePure_(fs, phaseIdx);
            updateMix_(fs, phaseIdx);
            VmUpToDate_[phaseIdx] = false;
        }
        else if (!(exceptQuantities & ParentType::Composition))
        {
            updateMix_(fs, phaseIdx);
            VmUpToDate_[phaseIdx] = false;
        }
        else if (!(exceptQuantities & ParentType::Pressure)) {
            VmUpToDate_[phaseIdx] = false;
        }
    }
 
protected:
    template <class FluidState>
    void updatePure_(const FluidState &fs, int phaseIdx)
    {
        Scalar T = fs.temperature(phaseIdx);
        Scalar p = fs.pressure(phaseIdx);
 
        switch (phaseIdx)
        {
            case oPhaseIdx: oilPhaseParams_.updatePure(T, p); break;
            case gPhaseIdx: gasPhaseParams_.updatePure(T, p); break;
        }
    }

    template <class FluidState>
    void updateMix_(const FluidState &fs, int phaseIdx)
    {
        Valgrind::CheckDefined(fs.averageMolarMass(phaseIdx));
        switch (phaseIdx)
        {
            case oPhaseIdx:
                oilPhaseParams_.updateMix(fs);
                break;
            case gPhaseIdx:
                gasPhaseParams_.updateMix(fs);
                break;
            case wPhaseIdx:
                break;
        }
    }
 
    template <class FluidState>
    void updateMolarVolume_(const FluidState &fs,
                            int phaseIdx)
    {
        VmUpToDate_[phaseIdx] = true;
 
        // calculate molar volume of the phase (we will need this for the
        // fugacity coefficients and the density anyway)
        switch (phaseIdx) {
            case gPhaseIdx: {
                // calculate molar volumes for the given composition. although
                // this isn't a Peng-Robinson parameter strictly speaking, the
                // molar volume appears in basically every quantity the fluid
                // system can get queried, so it is okay to calculate it
                // here...
                Vm_[gPhaseIdx] =
                    PengRobinson::computeMolarVolume(fs,
                                                     *this,
                                                     phaseIdx,
                                                     /*isGasPhase=*/true);
                break;
            }
            case oPhaseIdx: {
                // calculate molar volumes for the given composition. although
                // this isn't a Peng-Robinson parameter strictly speaking, the
                // molar volume appears in basically every quantity the fluid
                // system can get queried, so it is okay to calculate it
                // here...
                Vm_[oPhaseIdx] =
                    PengRobinson::computeMolarVolume(fs,
                                                     *this,
                                                     phaseIdx,
                                                     /*isGasPhase=*/false);
                break;
            }
            case wPhaseIdx: {
                // Density of water in the stock tank (i.e. atmospheric
                // pressure) is specified as 62.4 lb/ft^3 by the SPE-5
                // paper. Also 1 lb = 0.4535923 and 1 ft = 0.3048 m.
                const Scalar stockTankWaterDensity = 62.4 * 0.45359237 / 0.028316847;
                // Water compressibility is specified as 3.3e-6 per psi
                // overpressure, where 1 psi = 6894.7573 Pa
                Scalar overPressure = fs.pressure(wPhaseIdx) - 1.013e5; // [Pa]
                Scalar waterDensity =
                    stockTankWaterDensity * (1 + 3.3e-6*overPressure/6894.7573);
 
                // convert water density [kg/m^3] to molar volume [m^3/mol]
                Vm_[wPhaseIdx] = fs.averageMolarMass(wPhaseIdx)/waterDensity;
                break;
            }
            default:
                DUNE_THROW(Dune::InvalidStateException, "invalid phaseIdx " << phaseIdx);
        }
    }
 
    bool VmUpToDate_[numPhases];
    Scalar Vm_[numPhases];
 
    OilPhaseParams oilPhaseParams_;
    GasPhaseParams gasPhaseParams_;
};
 
} // end namespace Dumux
 
#endif