spe5.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:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_SPE5_FLUID_SYSTEM_HH
#define DUMUX_SPE5_FLUID_SYSTEM_HH
 
#include "spe5parametercache.hh"
 
#include <dumux/common/spline.hh>
#include <dumux/material/eos/pengrobinsonmixture.hh>
#include <dumux/io/name.hh>
 
namespace Dumux::FluidSystems {
 
template <class Scalar>
class Spe5
{
    using ThisType = FluidSystems::Spe5<Scalar>;
 
    using PengRobinsonMixture = Dumux::PengRobinsonMixture<Scalar, ThisType>;
    using PengRobinson = Dumux::PengRobinson<Scalar>;
 
public:
    using ParameterCache = Spe5ParameterCache<Scalar, ThisType>;
 
    /****************************************
     * Fluid phase parameters
     ****************************************/
 
    static const int numPhases = 3;
 
    static const int gPhaseIdx = 0;
    static const int wPhaseIdx = 1;
    static const int oPhaseIdx = 2;
 
    using H2O = Dumux::Components::H2O<Scalar>;
 
    static std::string phaseName(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        switch (phaseIdx)
        {
            case gPhaseIdx: return IOName::gaseousPhase();
            case wPhaseIdx: return IOName::aqueousPhase();
            case oPhaseIdx: return IOName::naplPhase();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    static constexpr bool isMiscible()
    { return true; }
 
    static constexpr bool isGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx == gPhaseIdx;
    }
 
    static constexpr bool isCompressible(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }
 
    static bool isIdealMixture(int phaseIdx)
    {
        // always use the reference oil for the fugacity coefficients,
        // so they cannot be dependent on composition and they the
        // phases thus always an ideal mixture
        return phaseIdx == wPhaseIdx;
    }
 
    static bool isIdealGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return false; // gas is not ideal here!
    }
 
    /****************************************
     * 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 componentName(int compIdx)
    {
        static std::string name[] = {
            H2O::name(),
            std::string("C1"),
            std::string("C3"),
            std::string("C6"),
            std::string("C10"),
            std::string("C15"),
            std::string("C20")
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return name[compIdx];
    }
 
    static Scalar molarMass(int compIdx)
    {
        static const Scalar M[] = {
            H2O::molarMass(),
            16.04e-3, // C1
            44.10e-3, // C3
            86.18e-3, // C6
            142.29e-3, // C10
            206.00e-3, // C15
            282.00e-3 // C20
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return M[compIdx];
    }
 
    static Scalar criticalTemperature(int compIdx)
    {
        static const Scalar Tcrit[] = {
            H2O::criticalTemperature(), // H2O
            343.0*5/9, // C1
            665.7*5/9, // C3
            913.4*5/9, // C6
            1111.8*5/9, // C10
            1270.0*5/9, // C15
            1380.0*5/9 // C20
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return Tcrit[compIdx];
    }
 
    static Scalar criticalPressure(int compIdx)
    {
        static const Scalar pcrit[] = {
            H2O::criticalPressure(), // H2O
            667.8 * 6894.7573, // C1
            616.3 * 6894.7573, // C3
            436.9 * 6894.7573, // C6
            304.0 * 6894.7573, // C10
            200.0 * 6894.7573, // C15
            162.0 * 6894.7573  // C20
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return pcrit[compIdx];
    }
 
    static Scalar criticalMolarVolume(int compIdx)
    {
        static const Scalar R = 8.314472;
        static const Scalar vcrit[] = {
            H2O::criticalMolarVolume(), // H2O
            0.290*R*criticalTemperature(1)/criticalPressure(1), // C1
            0.277*R*criticalTemperature(2)/criticalPressure(2), // C3
            0.264*R*criticalTemperature(3)/criticalPressure(3), // C6
            0.257*R*criticalTemperature(4)/criticalPressure(4), // C10
            0.245*R*criticalTemperature(5)/criticalPressure(5), // C15
            0.235*R*criticalTemperature(6)/criticalPressure(6) // C20
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return vcrit[compIdx];
    }
 
    static Scalar acentricFactor(int compIdx)
    {
        static const Scalar accFac[] = {
            0.344, // H2O (from Reid, et al.)
            0.0130, // C1
            0.1524, // C3
            0.3007, // C6
            0.4885, // C10
            0.6500, // C15
            0.8500  // C20
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return accFac[compIdx];
    }
 
    static Scalar interactionCoefficient(int comp1Idx, int comp2Idx)
    {
        using std::min;
        using std::max;
        int i = min(comp1Idx, comp2Idx);
        int j = 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()
    {
        PengRobinsonParamsMixture<Scalar, ThisType, gPhaseIdx, /*useSpe5=*/true> prParams;
 
        // find envelopes of the 'a' and 'b' parameters for the range
        // 273.15K <= T <= 373.15K and 10e3 Pa <= p <= 100e6 Pa. 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 minT = 273.15, maxT = 373.15;
        Scalar minP = 1e4, maxP = 100e6;
 
        Scalar minA = 1e100, maxA = -1e100;
        Scalar minB = 1e100, maxB = -1e100;
 
        prParams.updatePure(minT, minP);
        using std::min;
        using std::max;
        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = min(prParams.pureParams(compIdx).a(), minA);
            maxA = max(prParams.pureParams(compIdx).a(), maxA);
            minB = min(prParams.pureParams(compIdx).b(), minB);
            maxB = max(prParams.pureParams(compIdx).b(), maxB);
        }
 
        prParams.updatePure(maxT, minP);
        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = min(prParams.pureParams(compIdx).a(), minA);
            maxA = max(prParams.pureParams(compIdx).a(), maxA);
            minB = min(prParams.pureParams(compIdx).b(), minB);
            maxB = max(prParams.pureParams(compIdx).b(), maxB);
        }
 
        prParams.updatePure(minT, maxP);
        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = min(prParams.pureParams(compIdx).a(), minA);
            maxA = max(prParams.pureParams(compIdx).a(), maxA);
            minB = min(prParams.pureParams(compIdx).b(), minB);
            maxB = max(prParams.pureParams(compIdx).b(), maxB);
        }
 
        prParams.updatePure(maxT, maxP);
        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
            minA = min(prParams.pureParams(compIdx).a(), minA);
            maxA = max(prParams.pureParams(compIdx).a(), maxA);
            minB = min(prParams.pureParams(compIdx).b(), minB);
            maxB = max(prParams.pureParams(compIdx).b(), maxB);
        }
 
        PengRobinson::init(/*aMin=*/minA, /*aMax=*/maxA, /*na=*/100,
                           /*bMin=*/minB, /*bMax=*/maxB, /*nb=*/200);
    }
 
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          const ParameterCache &paramCache,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
    }
 
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState, const ParameterCache &paramCache, int phaseIdx)
    {
        return 1.0/paramCache.molarVolume(phaseIdx);
    }
 
    template <class FluidState>
    static Scalar viscosity(const FluidState &fs,
                            const ParameterCache &paramCache,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
        if (phaseIdx == gPhaseIdx) {
            // given by SPE-5 in table on page 64. we use a constant
            // viscosity, though...
            return 0.0170e-2 * 0.1;
        }
        else if (phaseIdx == wPhaseIdx)
            // given by SPE-5: 0.7 centi-Poise  = 0.0007 Pa s
            return 0.7e-2 * 0.1;
        else {
            assert(phaseIdx == oPhaseIdx);
            // given by SPE-5 in table on page 64. we use a constant
            // viscosity, though...
            return 0.208e-2 * 0.1;
        }
    }
 
    template <class FluidState>
    static Scalar fugacityCoefficient(const FluidState &fs,
                                      const ParameterCache &paramCache,
                                      int phaseIdx,
                                      int compIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
        assert(0 <= compIdx  && compIdx <= numComponents);
 
        if (phaseIdx == oPhaseIdx || phaseIdx == gPhaseIdx)
            return PengRobinsonMixture::computeFugacityCoefficient(fs,
                                                                   paramCache,
                                                                   phaseIdx,
                                                                   compIdx);
        else {
            assert(phaseIdx == wPhaseIdx);
            return henryCoeffWater_(compIdx, fs.temperature(wPhaseIdx))
                   / fs.pressure(wPhaseIdx);
        }
    }
 
 
    template <class FluidState>
    static Scalar diffusionCoefficient(const FluidState &fs,
                                       const ParameterCache &paramCache,
                                       int phaseIdx,
                                       int compIdx)
    { DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); }
 
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             const ParameterCache &paramCache,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
    { DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficients"); }
 
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fs,
                           const ParameterCache &paramCache,
                           int phaseIdx)
    { DUNE_THROW(Dune::NotImplemented, "Enthalpies"); }
 
    template <class FluidState>
    static Scalar thermalConductivity(const FluidState &fluidState,
                                      const ParameterCache &paramCache,
                                      int phaseIdx)
    { DUNE_THROW(Dune::NotImplemented, "Thermal conductivities"); }
 
    template <class FluidState>
    static Scalar heatCapacity(const FluidState &fluidState,
                               const ParameterCache &paramCache,
                               int phaseIdx)
    { DUNE_THROW(Dune::NotImplemented, "Heat capacities"); }
 
 
private:
    static Scalar henryCoeffWater_(int compIdx, Scalar 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: DUNE_THROW(Dune::InvalidStateException, "Unknown component index " << compIdx);
        }
    }
};
 
} // end namespace Dumux::FluidSystems
 
#endif