H2OAirXyleneFluidSystem.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_XYLENE_FLUID_SYSTEM_HPP
#define OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HPP

#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Air.hpp>
#include <opm/material/components/H2O.hpp>
#include <opm/material/components/Xylene.hpp>
#include <opm/material/components/TabulatedComponent.hpp>

#include <opm/material/binarycoefficients/H2O_Air.hpp>
#include <opm/material/binarycoefficients/H2O_Xylene.hpp>
#include <opm/material/binarycoefficients/Air_Xylene.hpp>

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

#include <string_view>

namespace Opm {

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

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

    typedef ::Opm::H2O<Scalar> H2O;
    typedef Xylene<Scalar> NAPL;
    typedef ::Opm::Air<Scalar> Air;

    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()
    { }

    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 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 bool isCompressible(unsigned phaseIdx)
    {
        return
            (phaseIdx == gasPhaseIdx)
            // gases are always compressible
            ? true
            : (phaseIdx == waterPhaseIdx)
            // the water component decides for the water phase...
            ? H2O::liquidIsCompressible()
            // the NAPL component decides for the napl phase...
            : NAPL::liquidIsCompressible();
    }

    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)
            // gases are always compressible
            ? H2O::molarMass()
            : (compIdx == airIdx)
            // the water component decides for the water comp...
            ? Air::molarMass()
            // the NAPL component decides for the napl comp...
            : (compIdx == NAPLIdx)
            ? NAPL::molarMass()
            : 1e30;
    }

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

            // See: Ochs 2008
            // \todo: proper citation
            const LhsEval& rholH2O = H2O::liquidDensity(T, p);
            const LhsEval& clH2O = rholH2O/H2O::molarMass();

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

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

        const LhsEval& pH2O = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p;
        const LhsEval& pAir = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*p;
        const LhsEval& pNAPL = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*p;
        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 auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& 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::simpleGasViscosity(T, p),
            NAPL::gasViscosity(T, NAPL::vaporPressure(T))
        };
        // molar masses
        const Scalar M[numComponents] = {
            H2O::molarMass(),
            Air::molarMass(),
            NAPL::molarMass()
        };

        const auto& xgAir = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
        const auto& xgH2O = decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
        const auto& 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)
    {
        if (phaseIdx==gasPhaseIdx) {
            const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
            const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

            const LhsEval& diffAC = BinaryCoeff::Air_Xylene::gasDiffCoeff(T, p);
            const LhsEval& diffWC = BinaryCoeff::H2O_Xylene::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){
            Scalar diffACl = 1.e-9; // BinaryCoeff::Air_Xylene::liquidDiffCoeff(temperature, pressure);
            Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure);
            Scalar 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:
                return (1.- xww)/(xwa/diffAWl + xwc/diffWCl);
            case airIdx:
                return (1.- xwc)/(xww/diffWCl + xwa/diffACl);
            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;
    }

    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 auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));

        if (phaseIdx == waterPhaseIdx) {
            if (compIdx == H2OIdx)
                return H2O::vaporPressure(T)/p;
            else if (compIdx == airIdx)
                return BinaryCoeff::H2O_Air::henry(T)/p;
            else if (compIdx == NAPLIdx)
                return BinaryCoeff::H2O_Xylene::henry(T)/p;
        }

        // 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.
        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;
        }

        // 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 auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& 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
            const LhsEval& hgc = NAPL::gasEnthalpy(T, p);
            const LhsEval& hgw = H2O::gasEnthalpy(T, p);
            const LhsEval& hga = Air::gasEnthalpy(T, p);

            LhsEval result = 0;
            result += hgw * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
            result += hga * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
            result += hgc * decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));

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

private:
    template <class LhsEval>
    static LhsEval waterPhaseDensity_(const LhsEval& T,
                                      const LhsEval& pw,
                                      const LhsEval& xww,
                                      const LhsEval& xwa,
                                      const LhsEval& xwc)
    {
        const LhsEval& rholH2O = H2O::liquidDensity(T, pw);
        const LhsEval& clH2O = rholH2O/H2O::molarMass();

        // this assumes each dissolved molecule displaces exactly one
        // water molecule in the liquid
        return clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass());
    }

    template <class LhsEval>
    static LhsEval gasPhaseDensity_(const LhsEval& T,
                                    const LhsEval& pg,
                                    const LhsEval& xgw,
                                    const LhsEval& xga,
                                    const LhsEval& xgc)
    { return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc); }

    template <class LhsEval>
    static LhsEval NAPLPhaseDensity_(const LhsEval& T, const LhsEval& pn)
    { return NAPL::liquidDensity(T, pn); }
};

} // namespace Opm

#endif