h2oheavyoil.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_HEAVYOIL_FLUID_SYSTEM_HH
#define DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
 
#include <dumux/material/idealgas.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/heavyoil.hh>
 
#include <dumux/material/binarycoefficients/h2o_heavyoil.hh>
 
#include <dumux/material/fluidsystems/base.hh>
 
#include <dumux/io/name.hh>
 
namespace Dumux {
namespace FluidSystems {
 
template <class Scalar,
          class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar> > >
class H2OHeavyOil
    : public Base<Scalar, H2OHeavyOil<Scalar, H2OType> >
{
    using ThisType = H2OHeavyOil<Scalar, H2OType>;
 
public:
    using HeavyOil = Dumux::Components::HeavyOil<Scalar>;
    using H2O = H2OType;
 
 
    static const int numPhases = 3;
    static const int numComponents = 2;
 
    static const int wPhaseIdx = 0; // index of the water phase
    static const int nPhaseIdx = 1; // index of the NAPL phase
    static const int gPhaseIdx = 2; // index of the gas phase
 
    static const int H2OIdx = 0;
    static const int NAPLIdx = 1;
 
    // export component indices to indicate the main component
    // of the corresponding phase at atmospheric pressure 1 bar
    // and room temperature 20°C:
    static const int wCompIdx = H2OIdx;
    static const int nCompIdx = NAPLIdx;
 
    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)
    {
        if (H2O::isTabulated)
        {
            H2O::init(tempMin, tempMax, nTemp,
                      pressMin, pressMax, nPress);
        }
    }
 
    static constexpr int getMainComponent(int phaseIdx)
    {
        // For the gas phase, choosing a main component appears to be
        // rather arbitrary. Motivated by the fact that the thermal conductivity
        // of the gas phase is set to the thermal conductivity of pure water,
        // water is chosen for now.
        if (phaseIdx == nPhaseIdx)
            return nCompIdx;
        else
            return wCompIdx;
    }
 
    static constexpr bool isMiscible()
    { return true; }
 
    static constexpr bool isGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx == gPhaseIdx;
    }
 
    static bool isIdealGas(int phaseIdx)
    { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && HeavyOil::gasIsIdeal(); }
 
    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 == gPhaseIdx)
            return true;
        else if (phaseIdx == wPhaseIdx)
            // the water component decides for the water phase...
            return H2O::liquidIsCompressible();
 
        // the NAPL component decides for the napl phase...
        return HeavyOil::liquidIsCompressible();
    }
 
    static std::string phaseName(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        switch (phaseIdx)
        {
            case wPhaseIdx: return IOName::aqueousPhase();
            case nPhaseIdx: return IOName::naplPhase();
            case gPhaseIdx: return IOName::gaseousPhase();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    static std::string componentName(int compIdx)
    {
        switch (compIdx) {
            case H2OIdx: return H2O::name();
            case NAPLIdx: return HeavyOil::name();
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    static Scalar molarMass(int compIdx)
    {
        switch (compIdx) {
            case H2OIdx: return H2O::molarMass();
            case NAPLIdx: return HeavyOil::molarMass();
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    using Base<Scalar, ThisType>::density;
    template <class FluidState>
    static Scalar density(const FluidState &fluidState, int phaseIdx)
    {
        if (phaseIdx == wPhaseIdx) {
            // See: doctoral thesis of Steffen Ochs 2007
            // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
            // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
 
            // This assumes each gas molecule displaces exactly one
            // molecule in the liquid.
 
            return H2O::liquidMolarDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx))
                   * (H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
                      + HeavyOil::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
        }
        else if (phaseIdx == nPhaseIdx) {
            // assume pure NAPL for the NAPL phase
            Scalar pressure = HeavyOil::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
            return HeavyOil::liquidDensity(fluidState.temperature(phaseIdx), pressure);
        }
 
        assert (phaseIdx == gPhaseIdx);
        Scalar pH2O =
            fluidState.moleFraction(gPhaseIdx, H2OIdx)  *
            fluidState.pressure(gPhaseIdx);
        Scalar pNAPL =
            fluidState.moleFraction(gPhaseIdx, NAPLIdx)  *
            fluidState.pressure(gPhaseIdx);
        return H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O)
               + HeavyOil::gasDensity(fluidState.temperature(phaseIdx), pNAPL);
    }
 
    using Base<Scalar, ThisType>::molarDensity;
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
    {
        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == nPhaseIdx)
        {
            return HeavyOil::liquidMolarDensity(temperature, pressure);
        }
        else if (phaseIdx == wPhaseIdx)
        {   // This assumes each gas molecule displaces exactly one
            // molecule in the liquid.
            return H2O::liquidMolarDensity(temperature, pressure);
        }
        else
        {
            return H2O::gasMolarDensity(temperature, fluidState.partialPressure(gPhaseIdx, H2OIdx))
                   + HeavyOil::gasMolarDensity(temperature, fluidState.partialPressure(gPhaseIdx, NAPLIdx));
        }
    }
 
    using Base<Scalar, ThisType>::viscosity;
    template <class FluidState>
    static Scalar viscosity(const FluidState &fluidState,
                            int phaseIdx)
    {
        if (phaseIdx == wPhaseIdx) {
            // assume pure water viscosity
            return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
                                        fluidState.pressure(phaseIdx));
        }
        else if (phaseIdx == nPhaseIdx) {
            // assume pure NAPL viscosity
            return HeavyOil::liquidViscosity(fluidState.temperature(phaseIdx),
                                             fluidState.pressure(phaseIdx));
        }
 
        assert (phaseIdx == gPhaseIdx);
 
        /* 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 Scalar mu[numComponents] = {
            h2oGasViscosityInMixture(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)),
            HeavyOil::gasViscosity(fluidState.temperature(phaseIdx), HeavyOil::vaporPressure(fluidState.temperature(phaseIdx)))
        };
 
        return mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
               + mu[NAPLIdx]*fluidState.moleFraction(gPhaseIdx, NAPLIdx);
    }
 
    using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
 
    {
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        // liquid phase
        if (phaseIdx == wPhaseIdx)
            return BinaryCoeff::H2O_HeavyOil::liquidDiffCoeff(T, p);
 
        // gas phase
        else if (phaseIdx == gPhaseIdx)
            return BinaryCoeff::H2O_HeavyOil::gasDiffCoeff(T, p);
        else
            DUNE_THROW(Dune::InvalidStateException,
                       "Non-existent binary diffusion coefficient for phase index "
                       << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::diffusionCoefficient;
    template <class FluidState>
    static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx)
    {
        // liquid phase
        if (phaseIdx == wPhaseIdx)
            return binaryDiffusionCoefficient(fluidState, phaseIdx, H2OIdx, NAPLIdx);
        // gas phase
        else if (phaseIdx == gPhaseIdx)
            return binaryDiffusionCoefficient(fluidState, phaseIdx, NAPLIdx, H2OIdx);
        else
            DUNE_THROW(Dune::InvalidStateException,
                       "Non-existent diffusion coefficient for phase index "<< phaseIdx);
    }
 
    template <class FluidState>
    static Scalar henryCoefficient(const FluidState &fluidState,
                                   int phaseIdx,
                                   int compIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
        assert(0 <= compIdx  && compIdx < numComponents);
 
        const Scalar T = fluidState.temperature(phaseIdx);
 
        if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx)
            return Dumux::BinaryCoeff::H2O_HeavyOil::henryOilInWater(T);
 
        else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx)
            return Dumux::BinaryCoeff::H2O_HeavyOil::henryWaterInOil(T);
 
        else
            DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
                                                     << " and component index " << compIdx);
    }
 
    template <class FluidState>
    static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
                                      int compIdx)
    {
        assert(0 <= compIdx  && compIdx < numComponents);
 
        const Scalar T = fluidState.temperature(phaseIdx);
        if (compIdx == NAPLIdx)
            return HeavyOil::vaporPressure(T);
        else if (compIdx == H2OIdx)
            return H2O::vaporPressure(T);
        else
            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
    }
 
    template <class FluidState>
    static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
                                            int phaseIdx,
                                            int compIdx)
    {
        assert(0 <= compIdx  && compIdx < numComponents);
 
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (compIdx == NAPLIdx)
            return HeavyOil::vaporTemperature(pressure);
        else if (compIdx == H2OIdx)
            return H2O::vaporTemperature(pressure);
        else
            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
    }
 
    using Base<Scalar, ThisType>::enthalpy;
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           int phaseIdx)
    {
        if (phaseIdx == wPhaseIdx) {
            return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
        }
        else if (phaseIdx == nPhaseIdx) {
            return HeavyOil::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
        }
        else if (phaseIdx == gPhaseIdx) {  // gas phase enthalpy depends strongly on composition
            Scalar hgc = HeavyOil::gasEnthalpy(fluidState.temperature(phaseIdx),
                                           fluidState.pressure(phaseIdx));
            Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx),
                                          fluidState.pressure(phaseIdx));
 
            Scalar result = 0;
            result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
            result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
 
            return result;
        }
        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 == wPhaseIdx)
        {
            if (componentIdx == H2OIdx)
                return H2O::liquidEnthalpy(T, p);
            else if (componentIdx == NAPLIdx)
                DUNE_THROW(Dune::NotImplemented, "The component enthalpy for NAPL in water is not implemented.");
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        else if (phaseIdx == nPhaseIdx)
        {
            if (componentIdx == H2OIdx)
                DUNE_THROW(Dune::NotImplemented, "The component enthalpy for water in NAPL is not implemented.");
            else if (componentIdx == NAPLIdx)
                return HeavyOil::liquidEnthalpy(T, p);
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        else if (phaseIdx == gPhaseIdx)
        {
            if (componentIdx == H2OIdx)
                return H2O::gasEnthalpy(T, p);
            else if (componentIdx == NAPLIdx)
                return HeavyOil::gasEnthalpy(T, p);
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::heatCapacity;
    template <class FluidState>
    static Scalar heatCapacity(const FluidState &fluidState,
                               int phaseIdx)
    {
        DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()");
    }
 
    using Base<Scalar, ThisType>::thermalConductivity;
    template <class FluidState>
    static Scalar thermalConductivity(const FluidState &fluidState,
                                      int phaseIdx)
    {
        const Scalar temperature  = fluidState.temperature(phaseIdx) ;
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == wPhaseIdx)
        {
            return H2O::liquidThermalConductivity(temperature, pressure);
        }
        else if (phaseIdx == nPhaseIdx)
        {
            return HeavyOil::liquidThermalConductivity(temperature, pressure);
        }
        else if (phaseIdx == gPhaseIdx)
        {
            return H2O::gasThermalConductivity(temperature, pressure);
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
};
} // end namespace FluidSystems
} // end namespace Dumux
 
#endif