components/brine.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:
/*****************************************************************************
 *   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 3 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/>.   *
 *****************************************************************************/
#ifndef DUMUX_BRINE_HH
#define DUMUX_BRINE_HH
 
#include <cmath>
 
#include <dune/common/math.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/nacl.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/base.hh>
#include <dumux/material/components/liquid.hh>
#include <dumux/material/components/gas.hh>
 
namespace Dumux::Components {

template <class Scalar,
          class H2O_Tabulated = Components::TabulatedComponent<Components::H2O<Scalar>>>
class Brine
: public Components::Base<Scalar, Brine<Scalar, H2O_Tabulated> >
, public Components::Liquid<Scalar, Brine<Scalar, H2O_Tabulated> >
, public Components::Gas<Scalar, Brine<Scalar, H2O_Tabulated> >
{
    using ThisType = Brine<Scalar, H2O_Tabulated>;
public:
    using H2O = Components::TabulatedComponent<Dumux::Components::H2O<Scalar>>;

    static constexpr Scalar R = Constants<Scalar>::R;

    static std::string name()
    { return "Brine"; }

    static Scalar salinity()
    {
        static const Scalar salinity = getParam<Scalar>("Brine.Salinity");
        return salinity;
    }

   static Scalar molarMass()
   {
       const Scalar M1 = H2O::molarMass();
       const Scalar M2 = Components::NaCl<Scalar>::molarMass(); // molar mass of NaCl [kg/mol]
       return M1*M2/(M2 + ThisType::salinity()*(M1 - M2));
   };

    static Scalar criticalTemperature()
    { return H2O::criticalTemperature(); }

    static Scalar criticalPressure()
    { return H2O::criticalPressure(); }

    static Scalar tripleTemperature()
    { return H2O::tripleTemperature(); }

    static Scalar triplePressure()
    { return H2O::triplePressure(); }

    static Scalar vaporPressure(Scalar temperature)
    {
        //calculate mole fraction
        const Scalar M1 = H2O::molarMass();
        const Scalar M2 = Components::NaCl<Scalar>::molarMass(); // molar mass of NaCl [kg/mol]
        const Scalar xNaClLiquid = - M1 * ThisType::salinity() / ((M2-M1) * ThisType::salinity() - M2);
 
        // Raoult's law, see Thomas Fetzer's Dissertation Eq. 2.11.
        return H2O::vaporPressure(temperature) * (1 - xNaClLiquid);
    }

    static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
    { return H2O::gasEnthalpy(temperature, pressure); }

    static const Scalar liquidEnthalpy(Scalar T, Scalar p)
    {
        /*Numerical coefficients from PALLISER*/
        static const Scalar f[] = {
            2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
        };
 
        /*Numerical coefficients 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 }
        };
 
        const Scalar theta = T - 273.15;
        const Scalar salSat = f[0] + f[1]*theta + f[2]*theta*theta + f[3]*theta*theta*theta;
 
        /*Regularization*/
        using std::min;
        using std::max;
        const Scalar salinity = min(max(ThisType::salinity(),0.0), salSat);
 
        const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
 
        /*DAUBERT and DANNER*/
        /*U=*/const Scalar 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 */
 
        const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
 
        using Dune::power;
        Scalar d_h = 0;
        for (int i = 0; i<=3; i++) {
            for (int j=0; j<=2; j++) {
                d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
            }
        }
 
        /* heat of dissolution for halite according to Michaelides 1971 */
        const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
 
        /* Enthalpy of brine without any dissolved gas */
        const Scalar h_ls1 =(1-salinity)*hw + salinity*h_NaCl + salinity*delta_h; /* kJ/kg */
        return h_ls1*1E3; /*J/kg*/
    }

    static const Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
    {
        const Scalar eps = temperature*1e-8;
        return (liquidEnthalpy(temperature + eps, pressure)- liquidEnthalpy(temperature, pressure))/eps;
    }

    static const Scalar gasHeatCapacity(Scalar temperature,
                                        Scalar pressure)
    {
        return H2O::gasHeatCapacity(temperature, pressure);
    }

    static const Scalar gasInternalEnergy(Scalar temperature,
                                          Scalar pressure)
    {
        return H2O::gasInternalEnergy(temperature, pressure);
    }

    static const Scalar liquidInternalEnergy(Scalar temperature,
                                             Scalar pressure)
    {
        return liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
    }

    static Scalar gasDensity(Scalar temperature, Scalar pressure)
    { return H2O::gasDensity(temperature, pressure); }

    static Scalar gasMolarDensity(Scalar temperature, Scalar pressure)
    { return H2O::gasMolarDensity(temperature, pressure); }

    static constexpr bool gasIsIdeal()
    { return H2O::gasIsIdeal(); }

    static constexpr bool gasIsCompressible()
    { return H2O::gasIsCompressible(); }

    static constexpr bool liquidIsCompressible()
    { return H2O::liquidIsCompressible(); }

    static Scalar liquidDensity(Scalar temperature, Scalar pressure)
    {
        using std::max;
        const Scalar TempC = temperature - 273.15;
        const Scalar pMPa = pressure/1.0E6;
        const Scalar salinity = max(0.0, ThisType::salinity());
 
        const Scalar rhow = H2O::liquidDensity(temperature, pressure);
 
        const Scalar  density = rhow +
                                1000*salinity*(
                                    0.668 +
                                    0.44*salinity +
                                    1.0E-6*(
                                        300*pMPa -
                                        2400*pMPa*salinity +
                                        TempC*(
                                            80.0 +
                                            3*TempC -
                                            3300*salinity -
                                            13*pMPa +
                                            47*pMPa*salinity)));
        assert(density > 0.0);
        return density;
    }

    static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
    { return liquidDensity(temperature, pressure)/molarMass(); }

   static Scalar gasPressure(Scalar temperature, Scalar density)
   { return H2O::gasPressure(temperature, density); }

   static Scalar liquidPressure(Scalar temperature, Scalar density)
   {
       // We use the Newton method for this. For the initial value we
       // assume the pressure to be 10% higher than the vapor
       // pressure
       Scalar pressure = 1.1*vaporPressure(temperature);
       const Scalar eps = pressure*1e-7;
 
       Scalar deltaP = pressure*2;
 
       using std::abs;
       for (int i = 0; i < 5 && abs(pressure*1e-9) < abs(deltaP); ++i) {
           Scalar f = liquidDensity(temperature, pressure) - density;
 
           Scalar df_dp;
           df_dp = liquidDensity(temperature, pressure + eps);
           df_dp -= liquidDensity(temperature, pressure - eps);
           df_dp /= 2*eps;
 
           deltaP = - f/df_dp;
 
           pressure += deltaP;
       }
       assert(pressure > 0.0);
       return pressure;
   }

    static Scalar gasViscosity(Scalar temperature, Scalar pressure)
    { return H2O::gasViscosity(temperature, pressure); };

    static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
    {
        // regularisation
        using std::max;
        temperature = max(temperature, 275.0);
        const Scalar salinity = max(0.0, ThisType::salinity());
 
        using std::pow;
        using Dune::power;
        using std::exp;
        const Scalar T_C = temperature - 273.15;
        const Scalar A = (0.42*power((pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
        const Scalar mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A); //[cP]
        assert(mu_brine > 0.0);
        return mu_brine/1000.0; //[Pa s]
    }

    static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
    {
            Scalar tempC = temperature-273.15;
            Scalar xNaCl = ThisType::salinity() * H2O::molarMass() / (ThisType::salinity() * H2O::molarMass() + (1-ThisType::salinity() )*Components::NaCl<Scalar>::molarMass()); // mole fraction of NaCl
            Scalar m = xNaCl/(H2O::molarMass()*(1- xNaCl)); // molality of NaCl
            Scalar S = 5844.3 * m / (1000 + 58.443 *m);
            Scalar contribNaClFactor = 1.0 - (2.3434e-3 - 7.924e-6*tempC + 3.924e-8*tempC*tempC)*S + (1.06e-5 - 2.0e-8*tempC + 1.2e-10*tempC*tempC)*S*S;
            return contribNaClFactor * H2O::liquidThermalConductivity(temperature, pressure);
    }
};
 
template <class Scalar, class H2O>
struct IsAqueous<Brine<Scalar, H2O>> : public std::true_type {};
 
} // end namespace Dumux::Components
 
#endif