H2OAirMesityleneFluidSystem.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_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HPP
#define OPM_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HPP

#include "BaseFluidSystem.hpp"
#include "NullParameterCache.hpp"

#include <opm/material/IdealGas.hpp>
#include <opm/material/components/N2.hpp>
#include <opm/material/components/Air.hpp>
#include <opm/material/components/H2O.hpp>
#include <opm/material/components/SimpleH2O.hpp>
#include <opm/material/components/Mesitylene.hpp>
#include <opm/material/components/TabulatedComponent.hpp>
#include <opm/material/binarycoefficients/H2O_Air.hpp>
#include <opm/material/binarycoefficients/H2O_N2.hpp>
#include <opm/material/binarycoefficients/H2O_Mesitylene.hpp>
#include <opm/material/binarycoefficients/Air_Mesitylene.hpp>

#include <string_view>

namespace Opm {

template <class Scalar>
class H2OAirMesityleneFluidSystem
    : public BaseFluidSystem<Scalar, H2OAirMesityleneFluidSystem<Scalar> >
{
    typedef H2OAirMesityleneFluidSystem<Scalar> ThisType;
    typedef BaseFluidSystem<Scalar, ThisType> Base;

    typedef ::Opm::H2O<Scalar> IapwsH2O;
    typedef TabulatedComponent<Scalar, IapwsH2O, /*alongVaporPressure=*/false> TabulatedH2O;

public:
    template <class Evaluation>
    struct ParameterCache : public NullParameterCache<Evaluation>
    {};

    typedef Mesitylene<Scalar> NAPL;

    typedef ::Opm::Air<Scalar> Air;

    //typedef SimpleH2O H2O;
    typedef TabulatedH2O H2O;
    //typedef IapwsH2O H2O;

    static const int numPhases = 3;
    static const int numComponents = 3;

    static const int waterPhaseIdx = 0;
    static const int naplPhaseIdx = 1;
    static const int gasPhaseIdx = 2;

    static const int H2OIdx = 0;
    static const int NAPLIdx = 1;
    static const int airIdx = 2;

    static void init()
    {
        init(/*tempMin=*/273.15,
             /*tempMax=*/623.15,
             /*numTemp=*/50,
             /*pMin=*/0.0,
             /*pMax=*/20e6,
             /*numP=*/50);
    }

    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
                     Scalar pressMin, Scalar pressMax, unsigned nPress)
    {
        if (H2O::isTabulated) {
            TabulatedH2O::init(tempMin, tempMax, nTemp,
                               pressMin, pressMax, nPress);
        }
    }

    static bool isLiquid(unsigned phaseIdx)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx != gasPhaseIdx;
    }

    static bool isIdealGas(unsigned phaseIdx)
    { return phaseIdx == gasPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }

    static bool isCompressible(unsigned phaseIdx)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        // gases are always compressible
        return (phaseIdx == gasPhaseIdx)
            ? true
            : (phaseIdx == waterPhaseIdx)
            ? H2O::liquidIsCompressible()
            : NAPL::liquidIsCompressible();
    }

    static bool isIdealMixture(unsigned /*phaseIdx*/)
    {
        //assert(0 <= phaseIdx && phaseIdx < numPhases);
        // we assume Henry's and Rault'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 std::string_view phaseName(unsigned phaseIdx)
    {
        switch (phaseIdx) {
        case waterPhaseIdx: return "water";
        case naplPhaseIdx: return "napl";
        case gasPhaseIdx: return "gas";
        };
        throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
    }

    static std::string_view componentName(unsigned compIdx)
    {
        switch (compIdx) {
        case H2OIdx: return H2O::name();
        case airIdx: return Air::name();
        case NAPLIdx: return NAPL::name();
        };
        throw std::logic_error("Invalid component index "+std::to_string(compIdx));
    }

    static Scalar molarMass(unsigned compIdx)
    {
        return
            (compIdx == H2OIdx)
            ? H2O::molarMass()
            : (compIdx == airIdx)
            ? Air::molarMass()
            : (compIdx == NAPLIdx)
            ? NAPL::molarMass()
            : -1e10;
        //throw std::logic_error("Invalid component index "+std::to_string(compIdx));
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval density(const FluidState& fluidState,
                           const ParameterCache<ParamCacheEval>& /*paramCache*/,
                           unsigned phaseIdx)
    {
        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));

        if (phaseIdx == waterPhaseIdx) {
            // See: Ochs 2008
            const LhsEval& p =
                H2O::liquidIsCompressible()
                ? decay<LhsEval>(fluidState.pressure(phaseIdx))
                : 1e30;

            const LhsEval& rholH2O = H2O::liquidDensity(T, p);
            const LhsEval& clH2O = rholH2O/H2O::molarMass();

            // this assumes each dissolved molecule displaces exactly one
            // water molecule in the liquid
            return
                clH2O*(H2O::molarMass()*decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) +
                       Air::molarMass()*decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx)) +
                       NAPL::molarMass()*decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)));
        }
        else if (phaseIdx == naplPhaseIdx) {
            // assume pure NAPL for the NAPL phase
            const LhsEval& p =
                NAPL::liquidIsCompressible()
                ? decay<LhsEval>(fluidState.pressure(phaseIdx))
                : 1e30;
            return NAPL::liquidDensity(T, p);
        }

        assert (phaseIdx == gasPhaseIdx);
        const LhsEval& pg = decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
        const LhsEval& pH2O = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg;
        const LhsEval& pAir = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg;
        const LhsEval& pNAPL = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg;
        return
            H2O::gasDensity(T, pH2O) +
            Air::gasDensity(T, pAir) +
            NAPL::gasDensity(T, pNAPL);
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval viscosity(const FluidState& fluidState,
                             const ParameterCache<ParamCacheEval>& /*paramCache*/,
                             unsigned phaseIdx)
    {
        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

        if (phaseIdx == waterPhaseIdx) {
            // assume pure water viscosity

            return H2O::liquidViscosity(T,
                                        p);
        }
        else if (phaseIdx == naplPhaseIdx) {
            // assume pure NAPL viscosity
            return NAPL::liquidViscosity(T, p);
        }

        assert (phaseIdx == gasPhaseIdx);

        /* Wilke method. See:
         *
         * See: R. Reid, et al.: The Properties of Gases and Liquids,
         * 4th edition, McGraw-Hill, 1987, 407-410
         * 5th edition, McGraw-Hill, 20001, p. 9.21/22
         *
         * in this case, we use a simplified version in order to avoid
         * computationally costly evaluation of sqrt and pow functions and
         * divisions
         * -- compare e.g. with Promo Class p. 32/33
         */
        const LhsEval mu[numComponents] = {
            H2O::gasViscosity(T, H2O::vaporPressure(T)),
            Air::gasViscosity(T, p),
            NAPL::gasViscosity(T, NAPL::vaporPressure(T))
        };
        // molar masses
        const Scalar M[numComponents] = {
            H2O::molarMass(),
            Air::molarMass(),
            NAPL::molarMass()
        };

        const LhsEval& xgAir = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
        const LhsEval& xgH2O = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
        const LhsEval& xgNapl = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
        const LhsEval& xgAW = xgAir + xgH2O;
        const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/xgAW;
        const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW;

        Scalar phiCAW = 0.3; // simplification for this particular system
        /* actually like this
         * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2)
         *                 / std::sqrt(8.*(1.+M[NAPLIdx]/MAW));
         */
        const LhsEval& phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);

        return (xgAW*muAW)/(xgAW + xgNapl*phiAWC) + (xgNapl*mu[NAPLIdx])/(xgNapl + xgAW*phiCAW);
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval diffusionCoefficient(const FluidState& /*fluidState*/,
                                        const ParameterCache<ParamCacheEval>& /*paramCache*/,
                                        unsigned /*phaseIdx*/,
                                        unsigned /*compIdx*/)
    {
        return 0;
#if 0
        typedef Opm::MathToolbox<typename FluidState::ValueType> FsToolbox;

        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        LhsEval diffCont;

        if (phaseIdx==gasPhaseIdx) {
            const LhsEval& diffAC = BinaryCoeff::Air_Mesitylene::gasDiffCoeff(T, p);
            const LhsEval& diffWC = BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p);
            const LhsEval& diffAW = BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);

            const LhsEval& xga = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
            const LhsEval& xgw = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
            const LhsEval& xgc = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));

            if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC);
            else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC);
            else if (compIdx==airIdx) throw std::logic_error("Diffusivity of air in the gas phase "
                                                             "is constraint by sum of diffusive fluxes = 0 !\n");
        }
        else if (phaseIdx==waterPhaseIdx){
            const LhsEval& diffACl = 1.e-9; // BinaryCoeff::Air_Mesitylene::liquidDiffCoeff(temperature, pressure);
            const LhsEval& diffWCl = 1.e-9; // BinaryCoeff::H2O_Mesitylene::liquidDiffCoeff(temperature, pressure);
            const LhsEval& diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);

            const LhsEval& xwa = decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
            const LhsEval& xww = decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
            const LhsEval& xwc = decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));

            switch (compIdx) {
            case NAPLIdx:
                diffCont = (1.- xww)/(xwa/diffAWl + xwc/diffWCl);
                return diffCont;
            case airIdx:
                diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl);
                return diffCont;
            case H2OIdx:
                throw std::logic_error("Diffusivity of water in the water phase "
                                       "is constraint by sum of diffusive fluxes = 0 !\n");
            };
        }
        else if (phaseIdx==naplPhaseIdx) {
            throw std::logic_error("Diffusion coefficients of "
                                   "substances in liquid phase are undefined!\n");
        }
        return 0;
#endif
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval fugacityCoefficient(const FluidState& fluidState,
                                       const ParameterCache<ParamCacheEval>& /*paramCache*/,
                                       unsigned phaseIdx,
                                       unsigned compIdx)
    {
        assert(phaseIdx < numPhases);
        assert(compIdx < numComponents);

        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(p);

        if (phaseIdx == waterPhaseIdx) {
            if (compIdx == H2OIdx)
                return H2O::vaporPressure(T)/p;
            else if (compIdx == airIdx)
                return BinaryCoeff::H2O_N2::henry(T)/p;
            else if (compIdx == NAPLIdx)
                return BinaryCoeff::H2O_Mesitylene::henry(T)/p;
            assert(false);
        }
        // for the NAPL phase, we assume currently that nothing is
        // dissolved. this means that the affinity of the NAPL
        // component to the NAPL phase is much higher than for the
        // other components, i.e. the fugacity cofficient is much
        // smaller.
        else if (phaseIdx == naplPhaseIdx) {
            const LhsEval& phiNapl = NAPL::vaporPressure(T)/p;
            if (compIdx == NAPLIdx)
                return phiNapl;
            else if (compIdx == airIdx)
                return 1e6*phiNapl;
            else if (compIdx == H2OIdx)
                return 1e6*phiNapl;
            assert(false);
        }

        // for the gas phase, assume an ideal gas when it comes to
        // fugacity (-> fugacity == partial pressure)
        assert(phaseIdx == gasPhaseIdx);
        return 1.0;
    }


    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval enthalpy(const FluidState& fluidState,
                            const ParameterCache<ParamCacheEval>& /*paramCache*/,
                            unsigned phaseIdx)
    {
        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

        if (phaseIdx == waterPhaseIdx) {
            return H2O::liquidEnthalpy(T, p);
        }
        else if (phaseIdx == naplPhaseIdx) {
            return NAPL::liquidEnthalpy(T, p);
        }
        else if (phaseIdx == gasPhaseIdx) {
            // gas phase enthalpy depends strongly on composition
            LhsEval result = 0;
            result += H2O::gasEnthalpy(T, p) * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
            result += NAPL::gasEnthalpy(T, p) * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
            result += Air::gasEnthalpy(T, p) * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));

            return result;
        }
        throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval thermalConductivity(const FluidState& fluidState,
                                       const ParameterCache<ParamCacheEval>& /*paramCache*/,
                                       unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        if (phaseIdx == waterPhaseIdx){ // water phase
            const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
            const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

            return H2O::liquidThermalConductivity(T, p);
        }
        else if (phaseIdx == gasPhaseIdx) { // gas phase
            const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
            const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

            return Air::gasThermalConductivity(T, p);
        }

        assert(phaseIdx == naplPhaseIdx);

        // Taken from:
        //
        // D. K. H. Briggs: "Thermal Conductivity of Liquids",
        // Ind. Eng. Chem., 1957, 49 (3), pp 418–421
        //
        // Convertion to SI units:
        // 344e-6 cal/(s cm K) = 0.0143964 J/(s m K)
        return 0.0143964;
    }
};
} // namespace Opm

#endif