region1.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_IAPWS_REGION1_HH
#define DUMUX_IAPWS_REGION1_HH
 
#include <cmath>
#include <iostream>
#include <dumux/common/exceptions.hh>
 
namespace Dumux {
namespace IAPWS {
template <class Scalar>
class Region1
{
public:
    static void checkValidityRange(Scalar temperature, Scalar pressure,
                                   const std::string& propertyName = "This property")
    {
        // actually this is:
        /* 273.15 <= temperature &&
           temperature <= 623.15 &&
           pressure >= vaporPressure(temperature) &&
           pressure <= 100e6 */
        if (temperature <= 623.15 && pressure <= 100e6)
            return;
 
        DUNE_THROW(NumericalProblem,
                   propertyName << " of water is only implemented for temperatures below 623.15K and "
                   "pressures below 100MPa. (T=" << temperature << ", p=" << pressure << ")");
    }
 
    static constexpr Scalar tau(Scalar temperature)
    { return 1386.0 / temperature; }
 
    static constexpr Scalar dTau_dt(Scalar temperature)
    { return - 1386.0 / (temperature*temperature); }
 
    static constexpr Scalar pi(Scalar pressure)
    { return pressure / 16.53e6; }
 
    static constexpr Scalar dPi_dp(Scalar pressure)
    { return 1.0 / 16.53e6; }
 
    static constexpr Scalar dp_dPi(Scalar pressure)
    { return 16.53e6; }
 
    static Scalar gamma(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        Scalar result = 0;
        for (int i = 0; i < 34; ++i) {
            result += n(i)*pow(7.1 - pi_, I(i))*pow(tau_ - 1.222, J(i));
        }
 
        return result;
    }
 
 
    static Scalar dGamma_dTau(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        using std::pow;
        Scalar result = 0.0;
        for (int i = 0; i < 34; i++) {
            result += n(i) *
                      pow(7.1 - pi_, I(i)) *
                      pow(tau_ - 1.222,  J(i)-1) *
                      J(i);
        }
 
        return result;
    }
 
    static Scalar dGamma_dPi(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        using std::pow;
        Scalar result = 0.0;
        for (int i = 0; i < 34; i++) {
            result += -n(i) *
                       I(i) *
                       pow(7.1 - pi_, I(i) - 1) *
                       pow(tau_ - 1.222, J(i));
        }
 
        return result;
    }
 
    static Scalar ddGamma_dTaudPi(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        using std::pow;
        Scalar result = 0.0;
        for (int i = 0; i < 34; i++) {
            result += -n(i) *
                       I(i) *
                       J(i) *
                       pow(7.1 - pi_, I(i) - 1) *
                       pow(tau_ - 1.222, J(i) - 1);
        }
 
        return result;
    }
 
    static Scalar ddGamma_ddPi(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        using std::pow;
        Scalar result = 0.0;
        for (int i = 0; i < 34; i++) {
            result += n(i) *
                      I(i) *
                      (I(i) - 1) *
                      pow(7.1 - pi_, I(i) - 2) *
                      pow(tau_ - 1.222, J(i));
        }
 
        return result;
    }
 
    static Scalar ddGamma_ddTau(Scalar temperature, Scalar pressure)
    {
        Scalar tau_ = tau(temperature);   /* reduced temperature */
        Scalar pi_ = pi(pressure);    /* reduced pressure */
 
        using std::pow;
        Scalar result = 0.0;
        for (int i = 0; i < 34; i++) {
            result += n(i) *
                      pow(7.1 - pi_, I(i)) *
                      J(i) *
                      (J(i) - 1) *
                      pow(tau_ - 1.222,  J(i) - 2);
        }
 
        return result;
    }
 
private:
    static Scalar n(int i)
    {
        constexpr Scalar n[34] = {
            0.14632971213167, -0.84548187169114, -0.37563603672040e1,
            0.33855169168385e1, -0.95791963387872, 0.15772038513228,
           -0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3,
           -0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1,
           -0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3,
           -0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5,
           -0.72694996297594e-15,-0.31679644845054e-4, -0.28270797985312e-5,
           -0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6,
           -0.14341729937924e-12,-0.40516996860117e-6, -0.12734301741641e-8,
           -0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19,
            0.26335781662795e-22,-0.11947622640071e-22, 0.18228094581404e-23,
           -0.93537087292458e-25
        };
        return n[i];
    }
 
    static short int I(int i)
    {
        constexpr short int I[34] = {
            0, 0, 0,
            0, 0, 0,
            0, 0, 1,
            1, 1, 1,
            1, 1, 2,
            2, 2, 2,
            2, 3, 3,
            3, 4, 4,
            4, 5, 8,
            8, 21, 23,
            29, 30, 31,
            32
        };
        return I[i];
    }
 
    static short int J(int i)
    {
        constexpr short int J[34] = {
             -2, -1, 0,
              1, 2, 3,
              4, 5, -9,
             -7, -1, 0,
              1, 3, -3,
              0, 1, 3,
             17, -4, 0,
              6, -5, -2,
             10, -8, -11,
             -6, -29, -31,
            -38, -39, -40,
            -41
        };
        return J[i];
    }
 
};
 
} // end namespace IAPWS
} // end namespace Dumux
 
#endif