Common.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:
/*
  Copyright (C) 2010-2013 by Andreas Lauser
  Copyright (C) 2012 by Philipp Nuske
  Copyright (C) 2010 by Felix Bode

  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/>.
*/
#ifndef OPM_IAPWS_COMMON_HPP
#define OPM_IAPWS_COMMON_HPP

#include <opm/material/Constants.hpp>
#include <opm/material/common/MathToolbox.hpp>

#include <cmath>

namespace Opm {
namespace IAPWS {

template <class Scalar>
class Common
{
public:
    static const Scalar molarMass;

    static const Scalar Rs;

    static const Scalar criticalTemperature;

    static const Scalar criticalPressure;

    static const Scalar criticalDensity;

    static const Scalar criticalMolarVolume;

    static const Scalar acentricFactor;

    static const Scalar tripleTemperature;

    static const Scalar triplePressure;

    template <class Evaluation>
    static Evaluation viscosity(const Evaluation& temperature, const Evaluation& rho)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        Evaluation rhoBar = rho/322.0;
        Evaluation TBar = temperature/criticalTemperature;

        // muBar = muBar_1
        const Scalar Hij[6][7] = {
            { 5.20094e-1, 2.22531e-1,-2.81378e-1, 1.61913e-1,-3.25372e-2, 0, 0 },
            { 8.50895e-2, 9.99115e-1,-9.06851e-1, 2.57399e-1, 0, 0, 0 },
            { -1.08374, 1.88797 ,-7.72479e-1, 0, 0, 0, 0 },
            { -2.89555e-1, 1.26613 ,-4.89837e-1, 0, 6.98452e-2, 0,-4.35673e-3 },
            { 0, 0,-2.57040e-1, 0, 0, 8.72102e-3, 0 },
            { 0, 1.20573e-1, 0, 0, 0, 0,-5.93264e-4 }
        };

        Evaluation tmp, tmp2, tmp3 = 1;
        Evaluation muBar = 0;
        for (int i = 0; i <= 5; ++i) {
            tmp = 0;
            tmp2 = 1;
            for (int j = 0; j <= 6; ++j) {
                tmp += Hij[i][j]*tmp2;
                tmp2 *= (rhoBar - 1);
            };
            muBar += tmp3 * tmp;
            tmp3 *= 1.0/TBar - 1;
        };
        muBar *= rhoBar;
        muBar = Toolbox::exp(muBar);

        // muBar *= muBar_0
        muBar  *= 100*Toolbox::sqrt(TBar);
        const Scalar H[4] = {
            1.67752, 2.20462, 0.6366564, -0.241605
        };

        tmp = 0, tmp2 = 1;
        for (int i = 0; i < 4; ++i) {
            tmp += H[i]/tmp2;
            tmp2 *= TBar;
        };
        muBar /= tmp;

        return 1e-6*muBar;
    }

    template <class Evaluation>
    static Evaluation thermalConductivityIAPWS(const Evaluation& T, const Evaluation& rho)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        static const Scalar thcond_tstar = 647.26 ;
        static const Scalar thcond_rhostar = 317.7 ;
        /*static const Scalar thcond_kstar = 1.0 ;*/

        static const Scalar thcond_b0 = -0.397070 ;
        static const Scalar thcond_b1 = 0.400302 ;
        static const Scalar thcond_b2 = 1.060000 ;
        static const Scalar thcond_B1 = -0.171587 ;
        static const Scalar thcond_B2 = 2.392190 ;

        static const Scalar thcond_c1 = 0.642857 ;
        static const Scalar thcond_c2 = -4.11717 ;
        static const Scalar thcond_c3 = -6.17937 ;
        static const Scalar thcond_c4 = 0.00308976 ;
        static const Scalar thcond_c5 = 0.0822994 ;
        static const Scalar thcond_c6 = 10.0932 ;

        static const Scalar thcond_d1 = 0.0701309 ;
        static const Scalar thcond_d2 = 0.0118520 ;
        static const Scalar thcond_d3 = 0.00169937 ;
        static const Scalar thcond_d4 = -1.0200 ;
        static const int thcond_a_count = 4;
        static const Scalar thcond_a[thcond_a_count] = {
            0.0102811
            ,0.0299621
            ,0.0156146
            ,-0.00422464
        };

        Evaluation Tbar = T / thcond_tstar;
        Evaluation rhobar = rho / thcond_rhostar;

        /* fast implementation... minimised calls to 'pow' routine... */
        Evaluation Troot = Toolbox::sqrt(Tbar);
        Evaluation Tpow = Troot;
        Evaluation lam = 0;

        for(int k = 0; k < thcond_a_count; ++k) {
            lam += thcond_a[k] * Tpow;
            Tpow *= Tbar;
        }

        lam +=
            thcond_b0 + thcond_b1
            * rhobar + thcond_b2
            * Toolbox::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));

        Evaluation DTbar = Toolbox::abs(Tbar - 1) + thcond_c4;
        Evaluation DTbarpow = Toolbox::pow(DTbar, 3./5);
        Evaluation Q = 2. + thcond_c5 / DTbarpow;

        Evaluation S;
        if(Tbar >= 1)
            S = 1. / DTbar;
        else
            S = thcond_c6 / DTbarpow;

        Evaluation rhobar18 = Toolbox::pow(rhobar, 1.8);
        Evaluation rhobarQ = Toolbox::pow(rhobar, Q);

        lam +=
            (thcond_d1 / Toolbox::pow(Tbar,10.0) + thcond_d2) * rhobar18 *
            Toolbox::exp(thcond_c1 * (1 - rhobar * rhobar18))
            + thcond_d3 * S * rhobarQ *
            Toolbox::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
            + thcond_d4 *
            Toolbox::exp(thcond_c2 * Toolbox::pow(Troot,3.0) + thcond_c3 / Toolbox::pow(rhobar,5.0));
        return /*thcond_kstar * */ lam;
    }
};

template <class Scalar>
const Scalar Common<Scalar>::molarMass = 18.01518e-3;
template <class Scalar>
const Scalar Common<Scalar>::Rs = Opm::Constants<Scalar>::R/molarMass;
template <class Scalar>
const Scalar Common<Scalar>::criticalTemperature = 647.096;
template <class Scalar>
const Scalar Common<Scalar>::criticalPressure = 22.064e6;
template <class Scalar>
const Scalar Common<Scalar>::criticalDensity = 322.0;
template <class Scalar>
const Scalar Common<Scalar>::criticalMolarVolume = molarMass/criticalDensity;
template <class Scalar>
const Scalar Common<Scalar>::acentricFactor = 0.344;
template <class Scalar>
const Scalar Common<Scalar>::tripleTemperature = 273.16;
template <class Scalar>
const Scalar Common<Scalar>::triplePressure = 611.657;

} // namespace IAPWS
} // namespace Opm

#endif