brineair.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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_BRINE_AIR_FLUID_SYSTEM_HH
#define DUMUX_BRINE_AIR_FLUID_SYSTEM_HH
 
#include <array>
#include <cassert>
#include <iomanip>
 
#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidstates/adapter.hh>
#include <dumux/material/fluidsystems/base.hh>
#include <dumux/material/fluidsystems/brine.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/nacl.hh>
#include <dumux/material/binarycoefficients/h2o_air.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
 
#include <dumux/common/exceptions.hh>
 
#include <dumux/io/name.hh>
 
#include "brine.hh"
 
namespace Dumux::FluidSystems {
 
template<bool fastButSimplifiedRelations = false>
struct BrineAirDefaultPolicy
{
    static constexpr bool useBrineDensityAsLiquidMixtureDensity() { return fastButSimplifiedRelations;}
    static constexpr bool useIdealGasDensity() { return fastButSimplifiedRelations; }
    static constexpr bool useKelvinVaporPressure() { return false; }
};
 
template <class Scalar,
          class H2Otype = Components::TabulatedComponent<Components::H2O<Scalar>>,
          class Policy = BrineAirDefaultPolicy<>>
class BrineAir
: public Base<Scalar, BrineAir<Scalar, H2Otype, Policy>>
{
    using ThisType = BrineAir<Scalar, H2Otype, Policy>;
    using IdealGas = Dumux::IdealGas<Scalar>;
 
public:
    using H2O = H2Otype;
    using Air = Components::Air<Scalar>;
    using NaCl = Components::NaCl<Scalar>;
 
    using Brine = Dumux::FluidSystems::Brine<Scalar, H2Otype>;
 
    using H2O_Air = BinaryCoeff::H2O_Air;
 
    using ParameterCache = NullParameterCache;
 
    /****************************************
     * Fluid phase related static parameters
     ****************************************/
    static constexpr int numPhases = 2;     // one liquid and one gas phase
    static constexpr int numComponents = 3; // H2O, Air, NaCl
 
    static constexpr int liquidPhaseIdx = 0; // index of the liquid phase
    static constexpr int gasPhaseIdx = 1;    // index of the gas phase
 
    static constexpr int phase0Idx = liquidPhaseIdx; // index of the first phase
    static constexpr int phase1Idx = gasPhaseIdx;    // index of the second phase
 
    // export component indices to indicate the main component
    // of the corresponding phase at atmospheric pressure 1 bar
    // and room temperature 20°C:
    static constexpr int H2OIdx = 0;
    static constexpr int AirIdx = 1;
    static constexpr int NaClIdx = 2;
    static constexpr int comp0Idx = H2OIdx;
    static constexpr int comp1Idx = AirIdx;
    static constexpr int comp2Idx = NaClIdx;
 
private:
    struct BrineAdapterPolicy
    {
        using FluidSystem = Brine;
 
        static constexpr int phaseIdx(int brinePhaseIdx) { return liquidPhaseIdx; }
        static constexpr int compIdx(int brineCompIdx)
        {
            switch (brineCompIdx)
            {
                case Brine::H2OIdx: return H2OIdx;
                case Brine::NaClIdx: return NaClIdx;
                default: return 0; // this will never be reached, only needed to suppress compiler warning
            }
        }
    };
 
    template<class FluidState>
    using BrineAdapter = FluidStateAdapter<FluidState, BrineAdapterPolicy>;
 
public:
 
    /****************************************
     * 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 isMiscible()
    { return true; }
 
    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);
        // ideal gases are always compressible
        if (phaseIdx == gasPhaseIdx)
            return true;
        // let brine decide for the liquid phase...
        return Brine::isCompressible(Brine::liquidPhaseIdx);
    }
 
    static 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
    }
 
    static constexpr int getMainComponent(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        if (phaseIdx == liquidPhaseIdx)
            return H2OIdx;
        else
            return AirIdx;
    }
 
    /****************************************
     * Component related static parameters
     ****************************************/
    static std::string componentName(int compIdx)
    {
        assert(0 <= compIdx && compIdx < numComponents);
        switch (compIdx)
        {
            case H2OIdx: return H2O::name();
            case AirIdx: return Air::name();
            case NaClIdx: return NaCl::name();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    static Scalar molarMass(int compIdx)
    {
        assert(0 <= compIdx && compIdx < numComponents);
        switch (compIdx)
        {
            case H2OIdx: return H2O::molarMass();
            case AirIdx: return Air::molarMass();
            case NaClIdx: return NaCl::molarMass();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }
 
    template <class FluidState>
    static Scalar vaporPressure(const FluidState& fluidState, int compIdx)
    {
        // The vapor pressure of the water is affected by the
        // salinity, thus, we forward to the interface of Brine here
        if (compIdx == H2OIdx)
        {
            if (!Policy::useKelvinVaporPressure())
                return Brine::vaporPressure(BrineAdapter<FluidState>(fluidState), Brine::H2OIdx);
            else
            {
                //additionally taking the influence of the capillary pressure on the vapor pressure, described by Kelvin's equation, into account.
                const auto t = fluidState.temperature(liquidPhaseIdx);
                const auto pc =  (fluidState.wettingPhase() == (int) liquidPhaseIdx)
                                 ? fluidState.pressure(gasPhaseIdx)-fluidState.pressure(liquidPhaseIdx)
                                 : fluidState.pressure(liquidPhaseIdx)-fluidState.pressure(gasPhaseIdx);
                //influence of the capillary pressure on the vapor pressure, the influence of the salt concentration is already taken into account in the brine vapour pressure
                return Brine::vaporPressure(BrineAdapter<FluidState>(fluidState), Brine::H2OIdx)
                        *exp( -pc / (Brine::molarDensity(fluidState, liquidPhaseIdx) * (Dumux::Constants<Scalar>::R*t)));
            }
        }
        else if (compIdx == NaClIdx)
            DUNE_THROW(Dune::NotImplemented, "NaCl::vaporPressure(t)");
        else
            DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
    }
 
    /****************************************
     * thermodynamic relations
     ****************************************/
    static void init()
    {
        init(/*tempMin=*/273.15,
             /*tempMax=*/800.0,
             /*numTemptempSteps=*/200,
             /*startPressure=*/-10,
             /*endPressure=*/20e6,
             /*pressureSteps=*/200);
    }
 
    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
                      Scalar pressMin, Scalar pressMax, unsigned nPress)
    {
        std::cout << "The brine-air fluid system was configured with the following policy:\n";
        std::cout << " - use brine density as liquid mixture density: " << std::boolalpha << Policy::useBrineDensityAsLiquidMixtureDensity() << "\n";
        std::cout << " - use ideal gas density: " << std::boolalpha << Policy::useIdealGasDensity() << 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, int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
 
        const auto T = fluidState.temperature(phaseIdx);
        const auto p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure brine
            if (Policy::useBrineDensityAsLiquidMixtureDensity())
                return Brine::density(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
 
            // assume one molecule of gas replaces one "brine" molecule
            else
                return Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx)
                       *(H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
                         + NaCl::molarMass()*fluidState.moleFraction(liquidPhaseIdx, NaClIdx)
                         + Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx));
        }
        else if (phaseIdx == phase1Idx)
        {
            // for the gas phase assume an ideal gas
            if (Policy::useIdealGasDensity())
                return IdealGas::density(fluidState.averageMolarMass(phase1Idx), T, p);
 
            // if useComplexRelations = true, compute density. NaCl is assumed
            // not to be present in gas phase, NaCl has only solid interfaces implemented
            return H2O::gasDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx))
                   + Air::gasDensity(T, fluidState.partialPressure(phase1Idx, AirIdx));
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::molarDensity;
    template <class FluidState>
    static Scalar molarDensity(const FluidState& fluidState, int phaseIdx)
    {
        if (phaseIdx == liquidPhaseIdx)
            return Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
        else if (phaseIdx == phase1Idx)
        {
            const Scalar T = fluidState.temperature(phaseIdx);
 
            // for the gas phase assume an ideal gas
            if (Policy::useIdealGasDensity())
                return IdealGas::molarDensity(T, fluidState.pressure(phaseIdx));
 
            // if useComplexRelations = true, compute density. NaCl is assumed
            // not to be present in gas phase, NaCl has only solid interfaces implemented
            return H2O::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx))
                   + Air::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, AirIdx));
        }
        else
            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);
 
        if (phaseIdx == liquidPhaseIdx)
            return Brine::viscosity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
        else
            return Air::gasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(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 == gasPhaseIdx)
            return 1.0;
 
        else if (phaseIdx == liquidPhaseIdx)
        {
            // TODO: Should we use the vapor pressure of the mixture (brine) here?
            //       The presence of NaCl lowers the vapor pressure, thus, we would
            //       expect the fugacity coefficient to be lower as well. However,
            //       with the fugacity coefficient being dependent on the salinity,
            //       the equation system for the phase equilibria becomes non-linear
            //       and our constraint solvers assume linear system of equations.
            if (compIdx == H2OIdx)
                return H2O::vaporPressure(T)/p;
 
            else if (compIdx == AirIdx)
                return BinaryCoeff::H2O_Air::henry(T)/p;
 
            // we assume nacl always stays in the liquid phase
            else
                return 0.0;
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::diffusionCoefficient;
    template <class FluidState>
    static Scalar diffusionCoefficient(const FluidState &fluidState,
                                       int phaseIdx,
                                       int compIdx)
    {
        DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
    }
 
    using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState& fluidState,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        assert(0 <= compIIdx && compIIdx < numComponents);
        assert(0 <= compJIdx && compJIdx < numComponents);
 
        const auto T = fluidState.temperature(phaseIdx);
        const auto p = fluidState.pressure(phaseIdx);
 
        if (compIIdx > compJIdx)
            std::swap(compIIdx, compJIdx);
 
        if (phaseIdx == liquidPhaseIdx)
        {
            if(compIIdx == H2OIdx && compJIdx == AirIdx)
                return H2O_Air::liquidDiffCoeff(T, p);
            else if (compIIdx == H2OIdx && compJIdx == NaClIdx)
                return Brine::binaryDiffusionCoefficient(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx, Brine::H2OIdx, Brine::NaClIdx);
            else
                DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
                                                 << compIIdx << " and " << compJIdx
                                                 << " in phase " << phaseIdx);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            if (compIIdx == H2OIdx && compJIdx == AirIdx)
                return H2O_Air::gasDiffCoeff(T, p);
 
            // NaCl is expected to never be present in the gas phase. we need to
            // return a diffusion coefficient that does not case numerical problems.
            // We choose a very small value here.
            else if (compIIdx == AirIdx && compJIdx == NaClIdx)
                return 1e-12;
 
            else
                DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
                                                 << compIIdx << " and " << compJIdx
                                                 << " in phase " << phaseIdx);
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
 
    using Base<Scalar, ThisType>::enthalpy;
    template <class FluidState>
    static Scalar enthalpy(const FluidState& fluidState, int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
 
        const Scalar T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
            return Brine::enthalpy(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
        else if  (phaseIdx == gasPhaseIdx)
        {
            // This assumes NaCl not to be present in the gas phase
            const Scalar XAir = fluidState.massFraction(gasPhaseIdx, AirIdx);
            const Scalar XH2O = fluidState.massFraction(gasPhaseIdx, H2OIdx);
            return XH2O*H2O::gasEnthalpy(T, p) + XAir*Air::gasEnthalpy(T, p);
        }
 
        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(gasPhaseIdx);
        const Scalar p = fluidState.pressure(gasPhaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
            DUNE_THROW(Dune::NotImplemented, "The component enthalpies in the liquid phase are not implemented.");
 
        else if (phaseIdx == gasPhaseIdx)
        {
            if (componentIdx == H2OIdx)
                return H2O::gasEnthalpy(T, p);
            else if (componentIdx == AirIdx)
                return Air::gasEnthalpy(T, p);
            else if (componentIdx == NaClIdx)
                DUNE_THROW(Dune::InvalidStateException, "Implementation assumes NaCl not to be present in gas phase");
            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)
    {
        if (phaseIdx == liquidPhaseIdx)
            return Brine::thermalConductivity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
        // This assumes NaCl not to be present in the gas phase
        else if (phaseIdx == gasPhaseIdx)
            return Air::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)) * fluidState.massFraction(gasPhaseIdx, AirIdx)
            + H2O::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)) * fluidState.massFraction(gasPhaseIdx, H2OIdx);
 
        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 T = fluidState.temperature(phaseIdx);
        const Scalar p = fluidState.pressure(phaseIdx);
 
        if (phaseIdx == liquidPhaseIdx)
            return Brine::heatCapacity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
 
        // We assume NaCl not to be present in the gas phase here
        else if (phaseIdx == gasPhaseIdx)
            return Air::gasHeatCapacity(T, p)*fluidState.massFraction(gasPhaseIdx, AirIdx)
                   + H2O::gasHeatCapacity(T, p)*fluidState.massFraction(gasPhaseIdx, H2OIdx);
 
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }
};
 
} // end namespace Dumux::FluidSystems
 
#endif