BrineCO2FluidSystem.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_BRINE_CO2_SYSTEM_HPP
#define OPM_BRINE_CO2_SYSTEM_HPP

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

#include <opm/common/Exceptions.hpp>

#include <opm/material/IdealGas.hpp>

#include <opm/material/components/Brine.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/components/CO2Tables.hpp>
#include <opm/material/components/SimpleCO2.hpp>
#include <opm/material/components/SimpleH2O.hpp>
#include <opm/material/components/TabulatedComponent.hpp>

#include <opm/material/binarycoefficients/H2O_CO2.hpp>
#include <opm/material/binarycoefficients/Brine_CO2.hpp>
#include <opm/material/binarycoefficients/H2O_N2.hpp>

#include <string_view>

namespace Opm {

template <class Scalar>
class BrineCO2FluidSystem
    : public BaseFluidSystem<Scalar, BrineCO2FluidSystem<Scalar> >
{
    typedef ::Opm::H2O<Scalar> H2O_IAPWS;
    typedef ::Opm::Brine<Scalar, H2O_IAPWS> Brine_IAPWS;
    typedef TabulatedComponent<Scalar, H2O_IAPWS> H2O_Tabulated;
    typedef TabulatedComponent<Scalar, Brine_IAPWS> Brine_Tabulated;

    typedef H2O_Tabulated H2O;

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

    typedef Brine_Tabulated Brine;
    typedef ::Opm::CO2<Scalar> CO2;

    typedef BinaryCoeff::Brine_CO2<Scalar, H2O, CO2> BinaryCoeffBrineCO2;

    /****************************************
     * Fluid phase related static parameters
     ****************************************/

    static const int numPhases = 2;

    static const int liquidPhaseIdx = 0;
    static const int gasPhaseIdx = 1;

    static std::string_view phaseName(unsigned phaseIdx)
    {
        static const std::string_view name[] = {
            "liquid",
            "gas"
        };

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

    static bool isLiquid(unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        return phaseIdx != gasPhaseIdx;
    }

    static bool isIdealGas(unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        if (phaseIdx == gasPhaseIdx)
            return CO2::gasIsIdeal();
        return false;
    }

    static bool isIdealMixture([[maybe_unused]] unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        return true;
    }

    static bool isCompressible([[maybe_unused]] unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);

        return true;
    }

    /****************************************
     * Component related static parameters
     ****************************************/
    static const int numComponents = 2;

    static const int BrineIdx = 0;
    static const int CO2Idx = 1;

    static std::string_view componentName(unsigned compIdx)
    {
        static const std::string_view name[] = {
            Brine::name(),
            CO2::name(),
        };

        assert(compIdx < numComponents);
        return name[compIdx];
    }

    static Scalar molarMass(unsigned compIdx)
    {
        assert(compIdx < numComponents);
        return (compIdx==BrineIdx)
            ? Brine::molarMass()
            : CO2::molarMass();
    }

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

    static void init()
    {
        init(/*startTemp=*/273.15, /*endTemp=*/623.15, /*tempSteps=*/50,
             /*startPressure=*/1e4, /*endPressure=*/40e6, /*pressureSteps=*/50);
    }

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

        // set the salinity of brine to the one used by the CO2 tables
        Brine_IAPWS::salinity = CO2::brineSalinity;

        if (Brine::isTabulated) {
            Brine_Tabulated::init(tempMin, tempMax, nTemp,
                                  pressMin, pressMax, nPress);
        }
    }

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

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

        if (phaseIdx == liquidPhaseIdx) {
            // use normalized composition for to calculate the density
            // (the relations don't seem to take non-normalized
            // compositions too well...)
            LhsEval xlBrine = min(1.0, max(0.0, decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, BrineIdx))));
            LhsEval xlCO2 = min(1.0, max(0.0,  decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, CO2Idx))));
            LhsEval sumx = xlBrine + xlCO2;
            xlBrine /= sumx;
            xlCO2 /= sumx;

            LhsEval result = liquidDensity_(temperature,
                                            pressure,
                                            xlBrine,
                                            xlCO2);

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

        assert(phaseIdx == gasPhaseIdx);

        // use normalized composition for to calculate the density
        // (the relations don't seem to take non-normalized
        // compositions too well...)
        LhsEval xgBrine = min(1.0, max(0.0, decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, BrineIdx))));
        LhsEval xgCO2 = min(1.0, max(0.0,  decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, CO2Idx))));
        LhsEval sumx = xgBrine + xgCO2;
        xgBrine /= sumx;
        xgCO2 /= sumx;

        LhsEval result = gasDensity_(temperature,
                                     pressure,
                                     xgBrine,
                                     xgCO2);
        Valgrind::CheckDefined(result);
        return result;
    }

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

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

        if (phaseIdx == liquidPhaseIdx) {
            // assume pure brine for the liquid phase. TODO: viscosity
            // of mixture
            LhsEval result = Brine::liquidViscosity(temperature, pressure);
            Valgrind::CheckDefined(result);
            return result;
        }

        assert(phaseIdx == gasPhaseIdx);
        LhsEval result = CO2::gasViscosity(getTableInstance(), temperature, pressure);
        Valgrind::CheckDefined(result);
        return result;
    }

    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);

        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;

        const LhsEval& temperature = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& pressure = decay<LhsEval>(fluidState.pressure(phaseIdx));
        assert(temperature > 0);
        assert(pressure > 0);

        // Activity model for salt-out effect in brine-CO2 mutual solubility
        // 2 = Duan & Sun as modified in Spycher & Pruess, Trans. Porous Media, (2009)
        const int activityModel = 3;

        // calulate the equilibrium composition for the given
        // temperature and pressure. TODO: calculateMoleFractions()
        // could use some cleanup.
        LhsEval xlH2O, xgH2O;
        LhsEval xlCO2, xgCO2;
        BinaryCoeffBrineCO2::calculateMoleFractions(getTableInstance(),
                                                    temperature,
                                                    pressure,
                                                    LhsEval(Brine_IAPWS::salinity),
                                                    /*knownPhaseIdx=*/-1,
                                                    xlCO2,
                                                    xgH2O,
                                                    activityModel);

        // normalize the phase compositions
        xlCO2 = max(0.0, min(1.0, xlCO2));
        xgH2O = max(0.0, min(1.0, xgH2O));

        xlH2O = 1.0 - xlCO2;
        xgCO2 = 1.0 - xgH2O;

        if (compIdx == BrineIdx) {
            Scalar phigH2O = 1.0;
            return phigH2O * xgH2O / xlH2O;
        }
        else {
            assert(compIdx == CO2Idx);

            Scalar phigCO2 = 1.0;
            return phigCO2 * xgCO2 / xlCO2;
        };
    }

    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*/)
    {
        const LhsEval& temperature = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& pressure = decay<LhsEval>(fluidState.pressure(phaseIdx));
        if (phaseIdx == liquidPhaseIdx)
            return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure);

        assert(phaseIdx == gasPhaseIdx);
        return BinaryCoeffBrineCO2::gasDiffCoeff(getTableInstance(), temperature, pressure);
    }

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

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

        if (phaseIdx == liquidPhaseIdx) {
            const LhsEval& XlCO2 = decay<LhsEval>(fluidState.massFraction(phaseIdx, CO2Idx));
            const LhsEval& result = liquidEnthalpyBrineCO2_(temperature,
                                                            pressure,
                                                            Brine_IAPWS::salinity,
                                                            XlCO2);
            Valgrind::CheckDefined(result);
            return result;
        }
        else {
            const LhsEval& XCO2 = decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, CO2Idx));
            const LhsEval& XBrine = decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, BrineIdx));

            LhsEval result = 0;
            result += XBrine * Brine::gasEnthalpy(temperature, pressure);
            result += XCO2 * CO2::gasEnthalpy(getTableInstance(), temperature, pressure);
            Valgrind::CheckDefined(result);
            return result;
        }
    }

    template <class FluidState, class LhsEval = typename FluidState::ValueType, class ParamCacheEval = LhsEval>
    static LhsEval thermalConductivity(const FluidState& /*fluidState*/,
                                       const ParameterCache<ParamCacheEval>& /*paramCache*/,
                                       unsigned phaseIdx)
    {
        // TODO way too simple!
        if (phaseIdx == liquidPhaseIdx)
            return  0.6; // conductivity of water[W / (m K ) ]

        // gas phase
        return 0.025; // conductivity of air [W / (m K ) ]
    }

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

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

        if(phaseIdx == liquidPhaseIdx){
            return H2O::liquidHeatCapacity(temperature, pressure);
        }
        else{
            return CO2::gasHeatCapacity(temperature, pressure);
        }
    }

private:
    static const CO2Tables<double>& getTableInstance() {
        static CO2Tables<double> instance;
        return instance;
    }

    template <class LhsEval>
    static LhsEval gasDensity_(const LhsEval& T,
                               const LhsEval& pg,
                               const LhsEval& xgH2O,
                               const LhsEval& xgCO2)
    {
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(pg);
        Valgrind::CheckDefined(xgH2O);
        Valgrind::CheckDefined(xgCO2);

        return CO2::gasDensity(getTableInstance(), T, pg);
    }

    /***********************************************************************/
    /*                                                                     */
    /* Total brine density with dissolved CO2                              */
    /* rho_{b,CO2} = rho_w + contribution(salt) + contribution(CO2)        */
    /*                                                                     */
    /***********************************************************************/
    template <class LhsEval>
    static LhsEval liquidDensity_(const LhsEval& T,
                                  const LhsEval& pl,
                                  const LhsEval& xlH2O,
                                  const LhsEval& xlCO2)
    {
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(pl);
        Valgrind::CheckDefined(xlH2O);
        Valgrind::CheckDefined(xlCO2);

        auto cast = [](const auto d)
        {
#if HAVE_QUAD
            if constexpr (std::is_same_v<decltype(d), const quad>)
                return static_cast<double>(d);
            else
#endif
            return d;
        };

        auto tostring = [cast](const auto& val) -> std::string
                        {
                            if constexpr (DenseAd::is_evaluation<LhsEval>::value) {
                                return std::to_string(cast(getValue(val.value())));
                            } else {
                                return std::to_string(cast(val));
                            }
                        };

        if(T < 273.15) {
            const std::string msg =
                "Liquid density for Brine and CO2 is only "
                "defined above 273.15K (is +"
                + tostring(T)+ "K)";
            throw NumericalProblem(msg);
        }
        if(pl >= 2.5e8) {
            const std::string msg =
                "Liquid density for Brine and CO2 is only "
                "defined below 250MPa (is "
                + tostring(pl) + "Pa)";
            throw NumericalProblem(msg);
        }

        const LhsEval& rho_brine = Brine::liquidDensity(T, pl);
        const LhsEval& rho_pure = H2O::liquidDensity(T, pl);
        const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlH2O, xlCO2);
        const LhsEval& contribCO2 = rho_lCO2 - rho_pure;

        return rho_brine + contribCO2;
    }

    template <class LhsEval>
    static LhsEval liquidDensityWaterCO2_(const LhsEval& temperature,
                                          const LhsEval& pl,
                                          const LhsEval& /*xlH2O*/,
                                          const LhsEval& xlCO2)
    {
        Scalar M_CO2 = CO2::molarMass();
        Scalar M_H2O = H2O::molarMass();

        const LhsEval& tempC = temperature - 273.15;        /* tempC : temperature in °C */
        const LhsEval& rho_pure = H2O::liquidDensity(temperature, pl);
        // calculate the mole fraction of CO2 in the liquid. note that xlH2O is available
        // as a function parameter, but in the case of a pure gas phase the value of M_T
        // for the virtual liquid phase can become very large
        const LhsEval xlH2O = 1.0 - xlCO2;
        const LhsEval& M_T = M_H2O * xlH2O + M_CO2 * xlCO2;
        const LhsEval& V_phi =
            (37.51 +
             tempC*(-9.585e-2 +
                    tempC*(8.74e-4 -
                           tempC*5.044e-7))) / 1.0e6;
        return 1/ (xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T));
    }

    template <class LhsEval>
    static LhsEval liquidEnthalpyBrineCO2_(const LhsEval& T,
                                           const LhsEval& p,
                                           Scalar S, // salinity
                                           const LhsEval& X_CO2_w)
    {
        /* X_CO2_w : mass fraction of CO2 in brine */

        /* same function as enthalpy_brine, only extended by CO2 content */

        /*Numerical coefficents from PALLISER*/
        static const Scalar f[] = {
            2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
        };

        /*Numerical coefficents from MICHAELIDES for the enthalpy of brine*/
        static const Scalar a[4][3] = {
            { 9633.6, -4080.0, +286.49 },
            { +166.58, +68.577, -4.6856 },
            { -0.90963, -0.36524, +0.249667E-1 },
            { +0.17965E-2, +0.71924E-3, -0.4900E-4 }
        };

        LhsEval theta, h_NaCl;
        LhsEval h_ls1, d_h;
        LhsEval delta_h;
        LhsEval delta_hCO2, hg, hw;

        theta = T - 273.15;

        // Regularization
        Scalar scalarTheta = scalarValue(theta);
        Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3]));
        if (S > S_lSAT)
            S = S_lSAT;

        hw = H2O::liquidEnthalpy(T, p) /1E3; /* kJ/kg */

        /*DAUBERT and DANNER*/
        /*U=*/h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
                        +((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02; /* kJ/kg */

        Scalar m = 1E3/58.44 * S/(1-S);
        d_h = 0;

        for (int i = 0; i <= 3; ++i) {
            for (int j = 0; j <= 2; ++j) {
                d_h = d_h + a[i][j] * pow(theta, static_cast<Scalar>(i)) * std::pow(m, j);
            }
        }
        /* heat of dissolution for halite according to Michaelides 1971 */
        delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;

        /* Enthalpy of brine without CO2 */
        h_ls1 =(1-S)*hw + S*h_NaCl + S*delta_h; /* kJ/kg */

        /* heat of dissolution for CO2 according to Fig. 6 in Duan and Sun 2003. (kJ/kg)
           In the relevant temperature ranges CO2 dissolution is
           exothermal */
        delta_hCO2 = (-57.4375 + T * 0.1325) * 1000/44;

        /* enthalpy contribution of CO2 (kJ/kg) */
        hg = CO2::gasEnthalpy(getTableInstance(), T, p)/1E3 + delta_hCO2;

        /* Enthalpy of brine with dissolved CO2 */
        return (h_ls1 - X_CO2_w*hw + hg*X_CO2_w)*1E3; /*J/kg*/
    }
};

} // namespace Opm

#endif