h2oair.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_H2O_AIR_SYSTEM_HH
#define DUMUX_H2O_AIR_SYSTEM_HH
 
#include <cassert>
#include <iomanip>
 
#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidsystems/base.hh>
 
#include <dumux/material/binarycoefficients/h2o_air.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>
 
#include <dumux/common/exceptions.hh>
 
#include <dumux/io/name.hh>
 
namespace Dumux {
namespace FluidSystems {
template<bool fastButSimplifiedRelations = false>
struct H2OAirDefaultPolicy
{
    static constexpr bool useH2ODensityAsLiquidMixtureDensity() { return fastButSimplifiedRelations; }
    static constexpr bool useIdealGasDensity() { return fastButSimplifiedRelations; }
    static constexpr bool useAirViscosityAsGasMixtureViscosity() { return fastButSimplifiedRelations; }
};
 
template <class Scalar,
          class H2Otype = Components::TabulatedComponent<Components::H2O<Scalar> >,
          class Policy = H2OAirDefaultPolicy<>,
          bool useKelvinVaporPressure = false>
class H2OAir
: public Base<Scalar, H2OAir<Scalar, H2Otype, Policy> >
{
    using ThisType = H2OAir<Scalar,H2Otype, Policy>;
    using IdealGas = Dumux::IdealGas<Scalar>;
 
public:
    using H2O = H2Otype;
    using Air = Dumux::Components::Air<Scalar>;
 
    static constexpr int numPhases = 2; 
    static constexpr int numComponents = 2; 
 
    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 AirIdx = 1; 
    static constexpr int comp0Idx = H2OIdx; 
    static constexpr int comp1Idx = AirIdx; 
    static constexpr int liquidCompIdx = H2OIdx; 
    static constexpr int gasCompIdx = AirIdx; 
 
    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 isMiscible()
    { return true; }
 
    static constexpr bool isGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx == gasPhaseIdx;
    }
 
    static constexpr 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);
        // ideal gases are always compressible
        if (phaseIdx == gasPhaseIdx)
            return true;
        // the water component decides for the liquid phase...
        return H2O::liquidIsCompressible();
    }
 
    static constexpr bool viscosityIsConstant(int phaseIdx)
    {
        // water decides for the liquid phase
        if (phaseIdx == liquidPhaseIdx)
            return H2O::liquidViscosityIsConstant();
        // air decides if policy is enabled
        else if (phaseIdx == gasPhaseIdx && Policy::useAirViscosityAsGasMixtureViscosity())
            return Air::gasViscosityIsConstant();
        // in general it depends on the mixture
        else
            return false;
    }
 
    static constexpr bool isIdealGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
 
        // let the fluids decide
        if (phaseIdx == gasPhaseIdx)
            return H2O::gasIsIdeal() && Air::gasIsIdeal();
        return false; // not a gas
    }
 
    /****************************************
     * Component related static parameters
     ****************************************/
    static std::string componentName(int compIdx)
    {
        switch (compIdx)
        {
            case H2OIdx: return H2O::name();
            case AirIdx: return Air::name();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    static Scalar molarMass(int compIdx)
    {
        switch (compIdx)
        {
            case H2OIdx: return H2O::molarMass();
            case AirIdx: return Air::molarMass();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    static Scalar criticalTemperature(int compIdx)
    {
        static const Scalar TCrit[] = {
            H2O::criticalTemperature(),
            Air::criticalTemperature()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return TCrit[compIdx];
    }
 
    static Scalar criticalPressure(int compIdx)
    {
        static const Scalar pCrit[] = {
            H2O::criticalPressure(),
            Air::criticalPressure()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return pCrit[compIdx];
    }
 
    template <class FluidState>
    static Scalar vaporPressure(const FluidState &fluidState, int compIdx)
    {
        if (compIdx == H2OIdx)
        {
            const auto t = fluidState.temperature(H2OIdx);
            // cppcheck-suppress internalAstError
            if constexpr (!useKelvinVaporPressure)
                return H2O::vaporPressure(t);
            else
            {
                const auto pc =  (fluidState.wettingPhase() == (int) H2OIdx)
                                 ? fluidState.pressure(AirIdx)-fluidState.pressure(H2OIdx)
                                 : fluidState.pressure(H2OIdx)-fluidState.pressure(AirIdx);
                return H2O::vaporPressure(t)*exp( -pc * molarMass(H2OIdx)
                                                      / density(fluidState, H2OIdx)
                                                      / (Dumux::Constants<Scalar>::R*t) );
            }
        }
        else if (compIdx == AirIdx)
            // return Air::vaporPressure(fluidState.temperature(AirIdx));
            DUNE_THROW(Dune::NotImplemented, "Air::vaporPressure(t)");
        else
            DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
    }
 
    static Scalar criticalMolarVolume(int compIdx)
    {
        DUNE_THROW(Dune::NotImplemented,
                   "H2OAirFluidSystem::criticalMolarVolume()");
    }
 
    static Scalar acentricFactor(int compIdx)
    {
        static const Scalar accFac[] = {
            H2O::acentricFactor(),
            Air::acentricFactor()
        };
 
        assert(0 <= compIdx && compIdx < numComponents);
        return accFac[compIdx];
    }
 
    /****************************************
     * thermodynamic relations
     ****************************************/
 
    static void init()
    {
        init(/*tempMin=*/273.15,
             /*tempMax=*/623.15,
             /*numTemp=*/100,
             /*pMin=*/-10.,
             /*pMax=*/20e6,
             /*numP=*/200);
    }
 
    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
                     Scalar pressMin, Scalar pressMax, unsigned nPress)
    {
        std::cout << "The H2O-air 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 air viscosity as gas mixture viscosity: " << std::boolalpha << Policy::useAirViscosityAsGasMixtureViscosity() << std::endl;
 
        if (H2O::isTabulated)
        {
            H2O::init(tempMin, tempMax, nTemp,
                      pressMin, pressMax, nPress);
        }
    }
 
    using Base<Scalar, ThisType>::density;
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          const int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == phase0Idx)
        {
            if (Policy::useH2ODensityAsLiquidMixtureDensity())
                // assume pure water
                return H2O::liquidDensity(T, p);
            else
            {
                // See: Eq. (7) in Class et al. (2002a)
                // This assumes each gas molecule displaces exactly one
                // molecule in the liquid.
                return H2O::liquidMolarDensity(T, p)
                       * (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
                          + Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx));
            }
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            if (Policy::useIdealGasDensity())
                // for the gas phase assume an ideal gas
            {
                const Scalar averageMolarMass = fluidState.averageMolarMass(gasPhaseIdx);
                return IdealGas::density(averageMolarMass, T, p);
            }
 
            return H2O::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, H2OIdx))
                   + Air::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, AirIdx));
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::molarDensity;
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == phase0Idx)
        {
            // assume pure water or that each gas molecule displaces exactly one
            // molecule in the liquid.
            return H2O::liquidMolarDensity(T, p);
        }
        else if (phaseIdx == phase1Idx)
        {
            if (Policy::useIdealGasDensity())
                // for the gas phase assume an ideal gas
            { return IdealGas::molarDensity(T, p); }
 
            return H2O::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx))
                   + Air::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, AirIdx));
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::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);
 
        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure water for the liquid phase
            return H2O::liquidViscosity(T, p);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            if(Policy::useAirViscosityAsGasMixtureViscosity()){
                return Air::gasViscosity(T, p);
            }
            else //using a complicated version of this fluid system
            {
                // Wilke method (Reid et al.):
                Scalar muResult = 0;
                const Scalar mu[numComponents] = {
                    h2oGasViscosityInMixture(T, p),
                    Air::gasViscosity(T, p)
                };
 
                // molar masses
                const Scalar M[numComponents] =  {
                    H2O::molarMass(),
                    Air::molarMass()
                };
 
                for (int i = 0; i < numComponents; ++i)
                {
                    Scalar divisor = 0;
                    using std::sqrt;
                    for (int j = 0; j < numComponents; ++j)
                    {
                        // 1 + (mu[i]/mu[j]^1/2 * (M[i]/M[j])^1/4)
                        Scalar phiIJ = 1 + sqrt(mu[i]/mu[j] * sqrt(M[j]/M[i]));
                        phiIJ = phiIJ * phiIJ / sqrt(8*(1 + M[i]/M[j]));
                        divisor += fluidState.moleFraction(phaseIdx, j)*phiIJ;
                    }
                    muResult += fluidState.moleFraction(phaseIdx, i)*mu[i] / divisor;
                }
                return muResult;
            }
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::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);
 
        if (phaseIdx == liquidPhaseIdx) {
            if (compIdx == H2OIdx)
                return vaporPressure(fluidState, compIdx)/p;
            return BinaryCoeff::H2O_Air::henry(T)/p;
        }
 
        // for the gas phase, assume an ideal gas when it comes to
        // fugacity (-> fugacity == partial pressure)
        return 1.0;
    }
 
    template <class FluidState>
    static Scalar relativeHumidity(const FluidState &fluidState)
    {
        return fluidState.partialPressure(gasPhaseIdx, comp0Idx)
               / H2O::vaporPressure(fluidState.temperature(gasPhaseIdx));
    }
 
    using Base<Scalar, ThisType>::diffusionCoefficient;
    template <class FluidState>
    static Scalar diffusionCoefficient(const FluidState &fluidState,
                                       int phaseIdx,
                                       int compIdx)
    {
        DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAir::diffusionCoefficient()");
    }
 
    using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
    {
        using std::swap;
        if (compIIdx > compJIdx)
            swap(compIIdx, compJIdx);
 
        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);
 
        // we are in the liquid phase
        if (phaseIdx == liquidPhaseIdx)
        {
            if (compIIdx == H2OIdx && compJIdx == AirIdx)
                return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p);
            else
                DUNE_THROW(Dune::InvalidStateException,
                           "Binary diffusion coefficient of components "
                            << compIIdx << " and " << compJIdx
                            << " in phase " << phaseIdx << " is undefined!\n");
        }
 
        // we are in the gas phase
        else if (phaseIdx == gasPhaseIdx)
        {
            if (compIIdx == H2OIdx && compJIdx == AirIdx)
                return BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
            else
                DUNE_THROW(Dune::InvalidStateException,
                           "Binary diffusion coefficient of components "
                           << compIIdx << " and " << compJIdx
                           << " in phase " << phaseIdx << " is undefined!\n");
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::enthalpy;
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           int phaseIdx)
    {
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
            return H2O::liquidEnthalpy(T, p);
 
        else if (phaseIdx == gasPhaseIdx)
            return H2O::gasEnthalpy(T, p)*fluidState.massFraction(gasPhaseIdx, H2OIdx)
                   + Air::gasEnthalpy(T, p)*fluidState.massFraction(gasPhaseIdx, AirIdx);
 
        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 == liquidPhaseIdx)
        {
            // the liquid enthalpy is constant
            return H2O::liquidEnthalpy(T, p);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            if (componentIdx == H2OIdx)
            {
                return H2O::gasEnthalpy(T, p);
            }
            else if (componentIdx == AirIdx)
            {
                return Air::gasEnthalpy(T, p);
            }
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::thermalConductivity;
    template <class FluidState>
    static Scalar thermalConductivity(const FluidState &fluidState,
                                      int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
        const Scalar temperature  = fluidState.temperature(phaseIdx) ;
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == liquidPhaseIdx)
        {
            return H2O::liquidThermalConductivity(temperature, pressure);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            return Air::gasThermalConductivity(temperature, pressure);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::heatCapacity;
    template <class FluidState>
    static Scalar heatCapacity(const FluidState &fluidState,
                               int phaseIdx)
    {
        const Scalar temperature  = fluidState.temperature(phaseIdx);
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == liquidPhaseIdx)
        {
            // influence of air is neglected
            return H2O::liquidHeatCapacity(temperature, pressure);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            return Air::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(gasPhaseIdx, AirIdx)
                   + H2O::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(gasPhaseIdx, H2OIdx);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
};
 
} // end namespace FluidSystems
} // end namespace Dumux
 
#endif