biomin.hh

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup Fluidsystems
 * \brief A fluid system with water and gas as phases and brine and CO2
 *        as components.
 */
#ifndef DUMUX_BIOMIN_FLUID_SYSTEM_HH
#define DUMUX_BIOMIN_FLUID_SYSTEM_HH

#include <dumux/common/parameters.hh>
#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidsystems/base.hh>

#include <dumux/material/components/brine.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/co2.hh>
#include <dumux/material/components/tabulatedcomponent.hh>

#include <dumux/material/components/calciumion.hh>
#include "../components/urea.hh"

#include <dumux/material/binarycoefficients/brine_co2.hh>

namespace Dumux {
namespace FluidSystems {

/*!
 * \ingroup Fluidsystems
 * \brief A compositional fluid with brine and carbon dioxide as
 *        components in both, the liquid and the gas (supercritical) phase,
 *        additional biomineralisation components (Ca and Urea) in the liquid phase
 *
 * This class provides acess to the Bio fluid system when no property system is used.
 * For Dumux users, using BioMinFluid<TypeTag> and the documentation therein is
 * recommended.
 *
 *  The user can provide their own material table for co2 properties.
 *  This fluidsystem is initialized as default with the tabulated version of
 *  water of the IAPWS-formulation, and the tabularized adapter to transfer
 *  this into brine.
 *  In the non-TypeTagged version, salinity information has to be provided with
 *  the init() methods.
 */
template <class Scalar,
          class CO2Table,
          class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
class BioMin
: public Base<Scalar, BioMin<Scalar, CO2Table, H2OType> >
{

    using ThisType = BioMin<Scalar, H2OType>;
    using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
    using Brine_CO2 = BinaryCoeff::Brine_CO2<Scalar, CO2Table>;

    using IdealGas = Dumux::IdealGas<Scalar>;

public:
    using CO2 = Components::CO2<Scalar, CO2Table>;
    using H2O = H2OType;
    using Ca = Components::CalciumIon<Scalar>;
    using Urea = Components::Urea<Scalar>;
    using Brine = Components::Brine<Scalar, H2O>;

    // the type of parameter cache objects. this fluid system does not
    // cache anything, so it uses Dumux::NullParameterCache
    using ParameterCache = Dumux::NullParameterCache;

    /****************************************
     * Fluid phase related static parameters
     ****************************************/
    static constexpr int numPhases = 2; // liquid and gas phases
    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
    /*!
     * \brief Return the human readable name of a fluid phase
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static std::string phaseName(int phaseIdx)
    {
        static std::string name[] =
        {
            "liquid",
            "gas"
        };

        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return name[phaseIdx];
    }

    /*!
     * \brief Returns whether the fluids are miscible
     */
    static constexpr bool isMiscible()
    { return true; }

    /*!
     * \brief Return whether a phase is liquid
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isLiquid(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx != gasPhaseIdx;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be an ideal mixture.
     *
     * We define an ideal mixture as a fluid phase where the fugacity
     * coefficients of all components times the pressure of the phase
     * are independent on the fluid composition. This assumption is true
     * if Henry's law and Raoult's law apply. If you are unsure what
     * this function should return, it is safe to return false. The
     * only damage done will be (slightly) increased computation times
     * in some cases.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isIdealMixture(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be compressible.
     *
     * Compressible means that the partial derivative of the density
     * to the fluid pressure is always larger than zero.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isCompressible(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be an ideal gas.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static bool isIdealGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        // let the fluids decide
        if (phaseIdx == gasPhaseIdx)
            return H2O::gasIsIdeal() && CO2::gasIsIdeal();
        return false; // not a gas
    }

    /****************************************
     * Component related static parameters
     ****************************************/
    static constexpr int numComponents = 4; // H2O/brine, CO2, Ca, urea
    static constexpr int H2OIdx = 0;
    static constexpr int CO2Idx = 1;
    static constexpr int CaIdx = 2;
    static constexpr int UreaIdx = 3;

    static constexpr int BrineIdx = H2OIdx;
    static constexpr int comp0Idx = BrineIdx;
    static constexpr int comp1Idx = CO2Idx;

    /*!
     * \brief Return the human readable name of a component
     *
     * \param compIdx The index of the component to consider
     */
    static std::string componentName(int compIdx)
    {
        switch (compIdx) {
        case BrineIdx: return Brine::name();
        case CO2Idx: return "CO2";
        case CaIdx: return Ca::name();
        case UreaIdx: return Urea::name();
        default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
        };
    }

    /*!
     * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
     *
     * \param compIdx The index of the component to consider
     */
    static Scalar molarMass(int compIdx)
    {
        switch (compIdx) {
        case H2OIdx: return H2O::molarMass();
        // actually, the molar mass of brine is only needed for diffusion
        // but since solutes are accounted for seperately
        // only the molar mass of water is returned.
        case CO2Idx: return CO2::molarMass();
        case CaIdx: return Ca::molarMass();
        case UreaIdx: return Urea::molarMass();
        default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
        };
    }

    /****************************************
     * thermodynamic relations
     ****************************************/

    static void init()
    {
        init(/*startTemp=*/295.15, /*endTemp=*/305.15, /*tempSteps=*/10,
             /*startPressure=*/1e4, /*endPressure=*/1e6, /*pressureSteps=*/200);

    }

    static void init(Scalar startTemp, Scalar endTemp, int tempSteps,
                     Scalar startPressure, Scalar endPressure, int pressureSteps)
    {
        std::cout << "Initializing tables for the pure-water properties.\n";
        H2O::init(startTemp, endTemp, tempSteps,
                  startPressure, endPressure, pressureSteps);
    }

    /*!
     * \brief Given a phase's composition, temperature, pressure, and
     *        the partial pressures of all components, return its
     *        density \f$\mathrm{[kg/m^3]}\f$.
     *
     * \param fluidState The fluid state
     * \param phaseIdx The index of the phase
     */
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          const ParameterCache &paramCache,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        switch (phaseIdx) {
            // assume pure brine for the liquid phase.
            case liquidPhaseIdx:
                return liquidDensity_(temperature,
                                      pressure,
                                      fluidState.moleFraction(liquidPhaseIdx, CO2Idx),
                                      fluidState.moleFraction(liquidPhaseIdx, H2OIdx),
                                      fluidState.massFraction(liquidPhaseIdx, CaIdx)); //consider density effect of dissolved calcium

            // assume pure CO2 for the gas phase.
            case gasPhaseIdx:
                return CO2::gasDensity(temperature, pressure);

            default:
                DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); break;
        }
    }

    /*!
     * \brief The molar density \f$\rho_{mol,\alpha}\f$
     *   of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
     *
     * The molar density is defined by the
     * mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
     *
     * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
     */
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState,
                               const ParameterCache &paramCache,
                               int phaseIdx)
    {
        const Scalar temperature = fluidState.temperature(phaseIdx);
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == liquidPhaseIdx)
        {
            return density(fluidState, paramCache, phaseIdx)
                   / fluidState.averageMolarMass(phaseIdx);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            // for the gas phase assume an ideal gas
            return CO2::gasMolarDensity(temperature, pressure);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

    /*!
     * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$.
     *
     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
     * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
     *
     * Equation given in: - Batzle & Wang (1992)
     *                    - cited by: Bachu & Adams (2002)
     *                    "Equations of State for basin geofluids"
     */
    template <class FluidState>
    static Scalar viscosity(const FluidState &fluidState,
                            const ParameterCache &paramCache,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == liquidPhaseIdx)
        {
            // assume brine with viscosity effect of Ca for the liquid phase.
            Scalar result = Brine::liquidViscosity(temperature, pressure);

            Valgrind::CheckDefined(result);
            return result;
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            // assume pure CO2 for the gas phase.
            Scalar result = CO2::gasViscosity(temperature, pressure);
            Valgrind::CheckDefined(result);
            return result;
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);

    }

    /*!
     * \brief Returns the fugacity coefficient [Pa] of a component in a
     *        phase.
     *
     * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of a
     * component \f$\kappa\f$ for a fluid phase \f$\alpha\f$ defines
     * the fugacity \f$f^\kappa_\alpha\f$ by the equation
     *
     * \f[
     f^\kappa_\alpha := \phi^\kappa_\alpha x^\kappa_\alpha p_\alpha\;.
     \f]
     *
     * The fugacity itself is just an other way to express the
     * chemical potential \f$\zeta^\kappa_\alpha\f$ of the component:
     *
     * \f[
     f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
     \f]
     * where \f$k_B\f$ is Boltzmann's constant.
     */
    template <class FluidState>
    static Scalar fugacityCoefficient(const FluidState &fluidState,
                                      const ParameterCache &paramCache,
                                      int phaseIdx,
                                      int compIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        assert(0 <= compIdx && compIdx < numComponents);

        if (phaseIdx == gasPhaseIdx)
            // use the fugacity coefficients of an ideal gas. the
            // actual value of the fugacity is not relevant, as long
            // as the relative fluid compositions are observed,
            return 1.0;

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);
        Scalar salinity = Brine::salinity(); // 0.1; //TODO major assumption in favor of runtime!
        //function is actually designed for use with NaCl not Ca.
        //Theoretically it should be: fluidState.massFraction(liquidPhaseIdx, CaIdx);

        assert(temperature > 0);
        assert(pressure > 0);
        // calulate the equilibrium composition for the given
        // temperature and pressure.
        Scalar xwH2O, xnH2O;
        Scalar xwCO2, xnCO2;
        Brine_CO2::calculateMoleFractions(temperature,
                                          pressure,
                                          salinity,
                                          /*knowgasPhaseIdx=*/-1,
                                          xwCO2,
                                          xnH2O);

        // normalize the phase compositions
        using std::min;
        using std::max;
        xwCO2 = max(0.0, min(1.0, xwCO2));
        xnH2O = max(0.0, min(1.0, xnH2O));

        xwH2O = 1.0 - xwCO2;
        xnCO2 = 1.0 - xnH2O;

        if (compIdx == BrineIdx)
        {
            Scalar phigH2O = 1.0;
            return phigH2O * xnH2O / xwH2O;
        }
        if (compIdx == CO2Idx)
        {
            Scalar phigCO2 = 1.0;
            return phigCO2 * xnCO2 / xwCO2;
        }
        else
            return 1/pressure; //all other components stay in the liquid phase
    }

    /*!
     * \brief Given the phase compositions, return the binary
     *        diffusion coefficient \f$\mathrm{[m^2/s]}\f$ of two components in a phase.
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     * \param compIIdx Index of the component i
     * \param compJIdx Index of the component j
     */
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             const ParameterCache &paramCache,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        assert(0 <= compIIdx && compIIdx < numComponents);
        assert(0 <= compJIdx && compJIdx < numComponents);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == liquidPhaseIdx)
        {
            assert(compIIdx == H2OIdx);
            Scalar result = 0.0;
            if(compJIdx == CO2Idx)
                result = Brine_CO2::liquidDiffCoeff(temperature, pressure);

            else if (compJIdx < numComponents) //Calcium and urea
                result = 1.587e-9;  //[m²/s]  educated guess, value for NaCl from J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
            else
                DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");

            Valgrind::CheckDefined(result);
            return result;
        }
        else
        {
            assert(phaseIdx == gasPhaseIdx);
            assert(compIIdx == CO2Idx);
            Scalar result = 0.0;
            if(compJIdx == H2OIdx)
                result = Brine_CO2::gasDiffCoeff(temperature, pressure);

            else if (compJIdx <numComponents) //Calcium and urea will stay in brine, no gaseous calcium or urea!
                result = 0.0;
            else
                DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");

            Valgrind::CheckDefined(result);
            return result;
        }
    };

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its specific enthalpy \f$\mathrm{[J/kg]}\f$.
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     *
     * See:
     * Class 2001:
     * Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien
     * Chapter 2.1.13 Innere Energie, Wäremekapazität, Enthalpie \cite A3:class:2001 <BR>
     *
     * Formula (2.42):
     * the specific enthalpy of a gasphase result from the sum of (enthalpies*mass fraction) of the components
     *
     */
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           const ParameterCache &paramCache,
                           int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure brine for the liquid phase.
            return Brine::liquidEnthalpy(temperature, pressure);
        }
        else
        {
            // assume pure CO2 for the gas phase.
            return CO2::gasEnthalpy(temperature, pressure);
        }
    };

private:
    //! calculate liquid density with respect to Water, CO2 and salt
    static Scalar liquidDensity_(Scalar T,
                                 Scalar pl,
                                 Scalar xwCO2,
                                 Scalar xwH2O,
                                 Scalar XlSal)
    {
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(pl);
        Valgrind::CheckDefined(XlSal);
        Valgrind::CheckDefined(xwCO2);

        if(T < 273.15)
        {
            DUNE_THROW(NumericalProblem,
                       "Liquid density for Brine and CO2 is only "
                       "defined above 273.15K (is" << T << ")");
        }
        if(pl >= 2.5e8)
        {
            DUNE_THROW(NumericalProblem,
                       "Liquid density for Brine and CO2 is only "
                       "defined below 250MPa (is" << pl << ")");
        }

        Scalar rho_brine = Brine::liquidDensity(T, pl);
        Scalar rho_pure = H2O::liquidDensity(T, pl);
        Scalar rho_lCO2 = liquidDensityWaterCO2_(T, pl, xwH2O, xwCO2);
        Scalar contribCO2 = rho_lCO2 - rho_pure;
        return rho_brine + contribCO2;
    }

    //! calculate liquid Density of water and CO2
    static Scalar liquidDensityWaterCO2_(Scalar temperature,
                                         Scalar pl,
                                         Scalar xwH2O,
                                         Scalar xwCO2)
    {
        const Scalar M_CO2 = CO2::molarMass();
        const Scalar M_H2O = H2O::molarMass();

        const Scalar tempC = temperature - 273.15;        /* tempC : temperature in °C */
        const Scalar rho_pure = H2O::liquidDensity(temperature, pl);
        xwH2O = 1.0 - xwCO2; // xwH2O is available, but in case of a pure gas phase
                             // the value of M_T for the virtual liquid phase can become very large
        const Scalar M_T = M_H2O * xwH2O + M_CO2 * xwCO2;
        const Scalar V_phi =
            (37.51 +
             tempC*(-9.585e-2 +
                    tempC*(8.74e-4 -
                           tempC*5.044e-7))) / 1.0e6;
        return 1 / (xwCO2 * V_phi/M_T + M_H2O * xwH2O / (rho_pure * M_T));
    }

};

} // end namespace FluidSystems
} // end namespace Dumux

#endif