h2on2o2.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 DUMUX_H2O_N2_O2_FLUID_SYSTEM_HH
#define DUMUX_H2O_N2_O2_FLUID_SYSTEM_HH
 
#include <cassert>
#include <iomanip>
 
#include <dumux/common/exceptions.hh>
 
#include <dumux/material/fluidsystems/base.hh>
#include <dumux/material/idealgas.hh>
#include <dumux/material/constants.hh>
 
#include <dumux/material/components/n2.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/o2.hh>
 
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/binarycoefficients/h2o_n2.hh>
#include <dumux/material/binarycoefficients/h2o_o2.hh>
#include <dumux/material/binarycoefficients/n2_o2.hh>
 
#include <dumux/io/name.hh>
 
namespace Dumux {
namespace FluidSystems {
template<bool fastButSimplifiedRelations = false>
struct H2ON2O2DefaultPolicy
{
    static constexpr bool useH2ODensityAsLiquidMixtureDensity() { return fastButSimplifiedRelations; }
    static constexpr bool useIdealGasDensity() { return fastButSimplifiedRelations; }
    static constexpr bool useN2ViscosityAsGasMixtureViscosity() { return fastButSimplifiedRelations; }
    static constexpr bool useN2HeatConductivityAsGasMixtureHeatConductivity() { return fastButSimplifiedRelations; }
    static constexpr bool useIdealGasHeatCapacities() { return fastButSimplifiedRelations; }
};
 
template <class Scalar, class Policy = H2ON2O2DefaultPolicy<>>
class H2ON2O2
    : public Base<Scalar, H2ON2O2<Scalar, Policy> >
{
    using ThisType = H2ON2O2<Scalar, Policy>;
    using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
 
    using IdealGas = Dumux::IdealGas<Scalar>;
    using Constants = Dumux::Constants<Scalar>;
    using TabulatedH2O = Components::TabulatedComponent<Dumux::Components::H2O<Scalar> >;
    using SimpleN2 = Dumux::Components::N2<Scalar>;
    using O2 = Dumux::Components::O2<Scalar>;
 
    using H2O = TabulatedH2O;
 
    using N2 = SimpleN2;
 
public:
    static constexpr int numPhases = 2; 
    static constexpr int numComponents = 3; 
    static constexpr int numSPhases = 0; // TODO: Remove
 
    static constexpr int liquidPhaseIdx = 0; 
    static constexpr int gasPhaseIdx = 1; 
    static constexpr int phase0Idx = liquidPhaseIdx; 
    static constexpr int phase1Idx = gasPhaseIdx; 
 
    static constexpr int H2OIdx = 0;
    static constexpr int N2Idx = 1;
    static constexpr int O2Idx = 2;
 
    static constexpr int comp0Idx = H2OIdx; // first major component
    static constexpr int comp1Idx = N2Idx; // second major component
    static constexpr int comp2Idx = O2Idx; // secondary component
 
    // main component at 20°C and 1 bar
    static constexpr int liquidPhaseMainCompIdx = H2OIdx;
    static constexpr int gasPhaseMainCompIdx = N2Idx;
 
    /****************************************
     * Fluid phase related static parameters
     ****************************************/
    static std::string phaseName(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        switch (phaseIdx)
        {
            case liquidPhaseIdx: return IOName::liquidPhase();
            case gasPhaseIdx: return IOName::gaseousPhase();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    static constexpr bool isGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx == gasPhaseIdx;
    }
 
    static bool isIdealMixture(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        // we assume Henry's and Raoult's laws for the water phase and
        // and no interaction between gas molecules of different
        // components, so all phases are ideal mixtures!
        return true;
    }
 
    static constexpr bool isCompressible(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        // gases are always compressible
        if (phaseIdx == gasPhaseIdx)
            return true;
        // the water component decides for the liquid phase...
        return H2O::liquidIsCompressible();
    }
 
    static bool isIdealGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        if (phaseIdx == gasPhaseIdx)
            // let the components decide
            return H2O::gasIsIdeal() && N2::gasIsIdeal() && O2::gasIsIdeal();
        return false; // not a gas
    }
 
    static constexpr bool isMiscible()
    { return true; }
 
    /****************************************
     * Component related static parameters
     ****************************************/
    static std::string componentName(int compIdx)
    {
        switch (compIdx)
        {
            case H2OIdx: return H2O::name();
            case N2Idx: return N2::name();
            case O2Idx: return O2::name();
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    static Scalar molarMass(int compIdx)
    {
        static const Scalar M[] = {
            H2O::molarMass(),
            N2::molarMass(),
            O2::molarMass()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return M[compIdx];
    }
 
    static Scalar criticalTemperature(int compIdx)
    {
        static const Scalar Tcrit[] = {
            H2O::criticalTemperature(),
            N2::criticalTemperature(),
            O2::criticalTemperature()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return Tcrit[compIdx];
    }
 
    static Scalar criticalPressure(int compIdx)
    {
        static const Scalar pcrit[] = {
            H2O::criticalPressure(),
            N2::criticalPressure(),
            O2::criticalPressure()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return pcrit[compIdx];
    }
 
    static Scalar criticalMolarVolume(int compIdx)
    {
        DUNE_THROW(Dune::NotImplemented,
                   "H2ON2O2FluidSystem::criticalMolarVolume()");
    }
 
    static Scalar acentricFactor(int compIdx)
    {
        static const Scalar accFac[] = {
            H2O::acentricFactor(),
            N2::acentricFactor(),
            O2::acentricFactor()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return accFac[compIdx];
    }
 
    template <class FluidState>
    static Scalar kelvinVaporPressure(const FluidState &fluidState,
                                      const int phaseIdx,
                                      const int compIdx,
                                      const Scalar radius)
    {
        assert(0 <= phaseIdx  && phaseIdx == liquidPhaseIdx);
        assert(0 <= compIdx  && compIdx == liquidPhaseMainCompIdx);
 
        Scalar T = fluidState.temperature(phaseIdx);
 
        Scalar vaporPressure = H2O::vaporPressure(T);
        Scalar exponent = molarMass(compIdx)/(density(fluidState, phaseIdx) * Constants::R * T);
        exponent *= (2 * surfaceTension(fluidState) / radius);
        using std::exp;
        Scalar kelvinVaporPressure = vaporPressure * exp(exponent);
 
        return kelvinVaporPressure;
    }
 
    template <class FluidState>
    static Scalar kelvinVaporPressure(const FluidState &fluidState,
                                      const int phaseIdx,
                                      const int compIdx)
    {
        assert(compIdx == liquidPhaseMainCompIdx && phaseIdx == liquidPhaseIdx);
 
        using std::exp;
        return fugacityCoefficient(fluidState, phaseIdx, compIdx)
               * fluidState.pressure(phaseIdx)
               * exp(-(fluidState.pressure(gasPhaseIdx)-fluidState.pressure(liquidPhaseIdx))
                      / density(fluidState, phaseIdx)
                      / (Dumux::Constants<Scalar>::R / molarMass(compIdx))
                      / fluidState.temperature());
    }
 
    template <class FluidState>
    static Scalar surfaceTension(const FluidState &fluidState)
    {
        const Scalar T = fluidState.temperature(); //K
        const Scalar B   = 0.2358 ; // [N/m]
        const Scalar T_c = H2O::criticalTemperature(); //K
        const Scalar mu  = 1.256;
        const Scalar b   = -0.625;
        //Equation to calculate surface Tension of Water According to IAPWS Release on Surface Tension from September 1994
        using std::pow;
        const Scalar surfaceTension = B*pow((1.-(T/T_c)),mu)*(1.+b*(1.-(T/T_c)));
        return surfaceTension; //surface Tension [N/m]
    }
    /****************************************
     * thermodynamic relations
     ****************************************/
 
    static void init()
    {
        init(/*tempMin=*/273.15,
             /*tempMax=*/623.15,
             /*numTemp=*/100,
             /*pMin=*/0.0,
             /*pMax=*/20e6,
             /*numP=*/200);
    }
 
    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
                     Scalar pressMin, Scalar pressMax, unsigned nPress)
    {
        std::cout << "The H2O-N2-O2 fluid system was configured with the following policy:\n";
        std::cout << " - use H2O density as liquid mixture density: " << std::boolalpha << Policy::useH2ODensityAsLiquidMixtureDensity() << "\n";
        std::cout << " - use ideal gas density: " << std::boolalpha << Policy::useIdealGasDensity() << "\n";
        std::cout << " - use N2 viscosity as gas mixture viscosity: " << std::boolalpha << Policy::useN2ViscosityAsGasMixtureViscosity() << "\n";
        std::cout << " - use N2 heat conductivity as gas mixture heat conductivity: " << std::boolalpha << Policy::useN2HeatConductivityAsGasMixtureHeatConductivity() << "\n";
        std::cout << " - use ideal gas heat capacities: " << std::boolalpha << Policy::useIdealGasHeatCapacities() << std::endl;
 
        if (H2O::isTabulated)
        {
            TabulatedH2O::init(tempMin, tempMax, nTemp,
                               pressMin, pressMax, nPress);
        }
    }
 
    using Base::density;
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure water
            if (Policy::useH2ODensityAsLiquidMixtureDensity())
                return H2O::liquidDensity(T, p);
 
            // See: Eq. (7) in Class et al. (2002a)
            // This assumes each gas molecule displaces exactly one
            // molecule in the liquid.
            else
                return H2O::liquidMolarDensity(T, p)
                       * (fluidState.moleFraction(liquidPhaseIdx, H2OIdx)*H2O::molarMass()
                          + fluidState.moleFraction(liquidPhaseIdx, N2Idx)*N2::molarMass()
                          + fluidState.moleFraction(liquidPhaseIdx, O2Idx)*O2::molarMass());
        }
 
        // gas phase
        else if (phaseIdx == gasPhaseIdx)
        {
 
            // for the gas phase assume an ideal gas
            using std::max;
            if (Policy::useIdealGasDensity())
                return IdealGas::molarDensity(T, p) * fluidState.averageMolarMass(gasPhaseIdx);
 
            // assume ideal mixture: steam, nitrogen and oxygen don't "see" each other
            else
                return H2O::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, H2OIdx))
                       + N2::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, N2Idx))
                       + O2::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, O2Idx));
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Unknown phase index " << phaseIdx);
    }
 
    using Base::molarDensity;
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
    {
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure water or that each gas molecule displaces exactly one
            // molecule in the liquid.
            return H2O::liquidMolarDensity(T, p);
        }
        else
        {
            if (Policy::useIdealGasDensity())
            {   //assume ideal gas
                return IdealGas::molarDensity(T,p);
            }
 
            return H2O::gasMolarDensity(T, fluidState.partialPressure(gasPhaseIdx, H2OIdx))
                   + N2::gasMolarDensity(T, fluidState.partialPressure(gasPhaseIdx, N2Idx))
                   + O2::gasMolarDensity(T, fluidState.partialPressure(gasPhaseIdx, O2Idx));
        }
    }
 
    using Base::viscosity;
    template <class FluidState>
    static Scalar viscosity(const FluidState &fluidState,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == liquidPhaseIdx) {
            // assume pure water for the liquid phase
            return H2O::liquidViscosity(T, p);
        }
 
        // gas phase
        if (Policy::useN2ViscosityAsGasMixtureViscosity())
        {
            // assume pure nitrogen for the gas phase
            return N2::gasViscosity(T, p);
        }
        else
        {
            // Wilke method (Reid et al.):
            Scalar muResult = 0;
            const Scalar mu[numComponents] = {
                h2oGasViscosityInMixture(T, p),
                N2::gasViscosity(T, p),
                O2::gasViscosity(T, p)
            };
 
            Scalar sumx = 0.0;
            using std::max;
            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                sumx += fluidState.moleFraction(phaseIdx, compIdx);
            sumx = max(1e-10, sumx);
 
            for (int i = 0; i < numComponents; ++i) {
                Scalar divisor = 0;
                using std::pow;
                using std::sqrt;
                for (int j = 0; j < numComponents; ++j) {
                    Scalar phiIJ = 1 + sqrt(mu[i]/mu[j]) * pow(molarMass(j)/molarMass(i), 1/4.0);
                    phiIJ *= phiIJ;
                    phiIJ /= sqrt(8*(1 + molarMass(i)/molarMass(j)));
                    divisor += fluidState.moleFraction(phaseIdx, j)/sumx * phiIJ;
                }
                muResult += fluidState.moleFraction(phaseIdx, i)/sumx * mu[i] / divisor;
            }
            return muResult;
        }
    }
 
    using Base::fugacityCoefficient;
    template <class FluidState>
    static Scalar fugacityCoefficient(const FluidState &fluidState,
                                      int phaseIdx,
                                      int compIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
        assert(0 <= compIdx  && compIdx < numComponents);
 
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == liquidPhaseIdx)
        {
            switch(compIdx){
                case H2OIdx: return H2O::vaporPressure(T)/p;
                case N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
                case O2Idx: return BinaryCoeff::H2O_O2::henry(T)/p;
            };
        }
 
        // for the gas phase, assume an ideal gas when it comes to
        // fugacity (-> fugacity == partial pressure)
        return 1.0;
    }
 
    using Base::diffusionCoefficient;
    template <class FluidState>
    static Scalar diffusionCoefficient(const FluidState &fluidState,
                                       int phaseIdx,
                                       int compIdx)
    {
        DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
    }
 
    using Base::binaryDiffusionCoefficient;
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
 
    {
        if (compIIdx > compJIdx)
        {
            using std::swap;
            swap(compIIdx, compJIdx);
        }
 
#ifndef NDEBUG
        if (compIIdx == compJIdx ||
            phaseIdx > numPhases - 1 ||
            compJIdx > numComponents - 1)
        {
            DUNE_THROW(Dune::InvalidStateException,
                       "Binary diffusion coefficient of components "
                       << compIIdx << " and " << compJIdx
                       << " in phase " << phaseIdx << " is undefined!\n");
        }
#endif
 
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == liquidPhaseIdx) {
            if (compIIdx == H2OIdx && compJIdx == N2Idx)
                return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
            if (compIIdx == H2OIdx && compJIdx == O2Idx)
                return BinaryCoeff::H2O_O2::liquidDiffCoeff(T, p);
            DUNE_THROW(Dune::InvalidStateException,
                       "Binary diffusion coefficient of components "
                       << compIIdx << " and " << compJIdx
                       << " in phase " << phaseIdx << " is undefined!\n");
        }
        // gas phase
        if (phaseIdx == gasPhaseIdx) {
            if (compIIdx == H2OIdx && compJIdx == N2Idx)
                return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
            if (compIIdx == H2OIdx && compJIdx == O2Idx)
                return BinaryCoeff::H2O_O2::gasDiffCoeff(T, p);
            if(compIIdx == N2Idx && compJIdx == O2Idx)
                return BinaryCoeff::N2_O2::gasDiffCoeff(T, p);
            DUNE_THROW(Dune::InvalidStateException,
                       "Binary diffusion coefficient of components "
                       << compIIdx << " and " << compJIdx
                       << " in phase " << phaseIdx << " is undefined!\n");
        }
 
        DUNE_THROW(Dune::InvalidStateException,
                  "Binary diffusion coefficient of components "
                  << compIIdx << " and " << compJIdx
                  << " in phase " << phaseIdx << " is undefined!\n");
    }
 
    using Base::enthalpy;
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           int phaseIdx)
    {
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == liquidPhaseIdx) {
            return H2O::liquidEnthalpy(T, p);
        }
        // gas phase
        else if (phaseIdx == gasPhaseIdx)
        {
            // assume ideal mixture: which means
            // that the total specific enthalpy is the sum of the
            // "partial specific enthalpies" of the components.
            Scalar hH2O =
                fluidState.massFraction(gasPhaseIdx, H2OIdx)
                * H2O::gasEnthalpy(T, p);
            Scalar hN2 =
                fluidState.massFraction(gasPhaseIdx, N2Idx)
                * N2::gasEnthalpy(T,p);
            Scalar hO2 =
                fluidState.massFraction(gasPhaseIdx, O2Idx)
                * O2::gasEnthalpy(T,p);
            return hH2O + hN2 + hO2;
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
    template <class FluidState>
    static Scalar componentEnthalpy(const FluidState &fluidState,
                                    int phaseIdx,
                                    int componentIdx)
    {
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == phase0Idx)
        {
            if (componentIdx == H2OIdx)
                return H2O::liquidEnthalpy(T, p);
            else if (componentIdx == N2Idx)
                DUNE_THROW(Dune::NotImplemented, "Component enthalpy of nitrogen in liquid phase");
            else if (componentIdx == O2Idx)
                DUNE_THROW(Dune::NotImplemented, "Component enthalpy of oxygen in liquid phase");
            else
                DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        else if (phaseIdx == phase1Idx)
        {
            if (componentIdx == H2OIdx)
                return H2O::gasEnthalpy(T, p);
            else if (componentIdx == N2Idx)
                return N2::gasEnthalpy(T, p);
            else if (componentIdx == O2Idx)
                return O2::gasEnthalpy(T, p);
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base::thermalConductivity;
    template <class FluidState>
    static Scalar thermalConductivity(const FluidState &fluidState,
                                      const int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
        Scalar temperature  = fluidState.temperature(phaseIdx) ;
        Scalar pressure = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
        {
            return H2O::liquidThermalConductivity(temperature, pressure);
        }
        else
        {
            Scalar lambdaPureN2 = N2::gasThermalConductivity(temperature, pressure);
            Scalar lambdaPureO2 = O2::gasThermalConductivity(temperature, pressure);
            if (!Policy::useN2HeatConductivityAsGasMixtureHeatConductivity())
            {
                Scalar xN2 = fluidState.moleFraction(phaseIdx, N2Idx);
                Scalar xO2 = fluidState.moleFraction(phaseIdx, O2Idx);
                Scalar xH2O = fluidState.moleFraction(phaseIdx, H2OIdx);
                Scalar lambdaN2 = xN2 * lambdaPureN2;
                Scalar lambdaO2 = xO2 * lambdaPureO2;
                Scalar partialPressure  = pressure * xH2O;
                Scalar lambdaH2O = xH2O * H2O::gasThermalConductivity(temperature, partialPressure);
                return lambdaN2 + lambdaH2O + lambdaO2;
            }
            else
                return lambdaPureN2;
        }
    }
 
    using Base::heatCapacity;
    template <class FluidState>
    static Scalar heatCapacity(const FluidState &fluidState,
                               int phaseIdx)
    {
        if (phaseIdx == liquidPhaseIdx) {
            return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
                                           fluidState.pressure(phaseIdx));
        }
 
        Scalar c_pN2;
        Scalar c_pO2;
        Scalar c_pH2O;
        // let the water and nitrogen components do things their own way
        if (!Policy::useIdealGasHeatCapacities()) {
            c_pN2 = N2::gasHeatCapacity(fluidState.temperature(phaseIdx),
                                        fluidState.pressure(phaseIdx)
                                        * fluidState.moleFraction(phaseIdx, N2Idx));
 
            c_pH2O = H2O::gasHeatCapacity(fluidState.temperature(phaseIdx),
                                          fluidState.pressure(phaseIdx)
                                          * fluidState.moleFraction(phaseIdx, H2OIdx));
            c_pO2 = O2::gasHeatCapacity(fluidState.temperature(phaseIdx),
                                        fluidState.pressure(phaseIdx)
                                        * fluidState.moleFraction(phaseIdx, O2Idx));
        }
        else {
            // assume an ideal gas for both components. See:
            //
            //http://en.wikipedia.org/wiki/Heat_capacity
            Scalar c_vN2molar = Constants::R*2.39;
            Scalar c_pN2molar = Constants::R + c_vN2molar;
 
            Scalar c_vO2molar = Constants::R*2.43;
            Scalar c_pO2molar = Constants::R + c_vO2molar;
 
            Scalar c_vH2Omolar = Constants::R*3.37; // <- correct??
            Scalar c_pH2Omolar = Constants::R + c_vH2Omolar;
 
            c_pN2 = c_pN2molar/molarMass(N2Idx);
            c_pO2 = c_pO2molar/molarMass(O2Idx);
            c_pH2O = c_pH2Omolar/molarMass(H2OIdx);
        }
 
        // mangle all components together
        return
            c_pH2O*fluidState.massFraction(gasPhaseIdx, H2OIdx)
            + c_pN2*fluidState.massFraction(gasPhaseIdx, N2Idx)
            + c_pO2*fluidState.massFraction(gasPhaseIdx, O2Idx);
    }
 
};
 
} // end namespace FluidSystems
} // end namespace Dumux
 
#endif