Region4.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) 2009-2013 by Andreas Lauser
  Copyright (C) 2012 by Benjamin Faigle
  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_REGION4_HPP
#define OPM_IAPWS_REGION4_HPP

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

#include <cmath>

namespace Opm {
namespace IAPWS {

template <class Scalar>
class Region4
{
public:
    template <class Evaluation>
    static Evaluation saturationPressure(const Evaluation& temperature)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        static const Scalar n[10] = {
            0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
            0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
            -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
            0.65017534844798e3
        };

        const Evaluation& sigma = temperature + n[8]/(temperature - n[9]);

        const Evaluation& A = (sigma + n[0])*sigma + n[1];
        const Evaluation& B = (n[2]*sigma + n[3])*sigma + n[4];
        const Evaluation& C = (n[5]*sigma + n[6])*sigma + n[7];

        Evaluation tmp = 2*C/(Toolbox::sqrt(B*B - 4*A*C) - B);
        tmp *= tmp;
        tmp *= tmp;

        return 1e6*tmp;
    }

    template <class Evaluation>
    static Evaluation vaporTemperature(const Evaluation& pressure)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        static const Scalar n[10] = {
            0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
            0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
            -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
            0.65017534844798e3
        };
        const Evaluation& beta = pow((pressure/1e6 /*from Pa to MPa*/), (1./4.));
        const Evaluation& beta2 = pow(beta, 2.);
        const Evaluation& E = beta2 + n[2] * beta + n[5];
        const Evaluation& F = n[0]*beta2 + n[3]*beta + n[6];
        const Evaluation& G = n[1]*beta2 + n[4]*beta + n[7];

        const Evaluation& D = ( 2.*G)/(-F -Toolbox::sqrt(pow(F,2.) - 4.*E*G));

        const Evaluation& temperature = (n[9] + D - Toolbox::sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;

        return temperature;
    }
};

} // namespace IAPWS
} // namespace Opm

#endif