Xylene.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 Vishal Jambhekar
  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_XYLENE_HPP
#define OPM_XYLENE_HPP

#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Component.hpp>
#include <opm/material/Constants.hpp>

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

#include <cmath>

namespace Opm {
template <class Scalar>
class Xylene : public Component<Scalar, Xylene<Scalar> >
{
    typedef Constants<Scalar> Consts;
    typedef Opm::IdealGas<Scalar> IdealGas;

public:
    static const char *name()
    { return "xylene"; }

    static Scalar molarMass()
    { return 0.106; }

    static Scalar criticalTemperature()
    { return 617.1; }

    static Scalar criticalPressure()
    { return 35.4e5; }

    static Scalar tripleTemperature()
    { OPM_THROW(std::runtime_error, "Not implemented: tripleTemperature for xylene"); }

    static Scalar triplePressure()
    { OPM_THROW(std::runtime_error, "Not implemented: triplePressure for xylene"); }

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

        const Scalar A = 7.00909;;
        const Scalar B = 1462.266;;
        const Scalar C = 215.110;;

        return 100*1.334*Toolbox::pow(10.0, (A - (B/(temperature - 273.15 + C))));
    }

    template <class Evaluation>
    static Evaluation spHeatCapLiquidPhase(const Evaluation& temperature, const Evaluation& /*pressure*/)
    {
        Evaluation CH3,C6H5,H;
        // after Reid et al. : Missenard group contrib. method (s. example 5-8)
        // Xylene: C9H12  : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was too much!)
        // linear interpolation between table values [J/(mol K)]

        if(temperature < 298.0){                              // take care: extrapolation for Temperature<273
            H = 13.4 + 1.2*(temperature - 273.0)/25.0;        // 13.4 + 1.2 = 14.6 = H(T=298K) i.e. interpolation of table values 273<T<298
            CH3 = 40.0 + 1.6*(temperature - 273.0)/25.0;    // 40 + 1.6 = 41.6 = CH3(T=298K)
            C6H5 = 113.0 + 4.2*(temperature - 273.0)/25.0; // 113 + 4.2 = 117.2 = C6H5(T=298K)
        }
        else if(temperature < 323.0){
            H = 14.6 + 0.9*(temperature - 298.0)/25.0;        // i.e. interpolation of table values 298<T<323
            CH3 = 41.6 + 1.9*(temperature - 298.0)/25.0;
            C6H5 = 117.2 + 6.2*(temperature - 298.0)/25.0;
        }
        else if(temperature < 348.0){
            H = 15.5 + 1.2*(temperature - 323.0)/25.0;        // i.e. interpolation of table values 323<T<348
            CH3 = 43.5 + 2.3*(temperature - 323.0)/25.0;
            C6H5 = 123.4 + 6.3*(temperature - 323.0)/25.0;
        }
        else {
            H = 16.7 + 2.1*(temperature - 348.0)/25.0;         // i.e. interpolation of table values 348<T<373
            CH3 = 45.8 + 2.5*(temperature - 348.0)/25.0;        // take care: extrapolation for Temperature>373
            C6H5 = 129.7 + 6.3*(temperature - 348.0)/25.0;        // most likely leads to underestimation
        }

        return (C6H5 + 2*CH3 - H)/molarMass();// J/(mol K) -> J/(kg K)
    }


    template <class Evaluation>
    static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
    {
        // Gauss quadrature rule:
        // Interval: [0K; temperature (K)]
        // Gauss-Legendre-Integration with variable transformation:
        // \int_a^b f(T) dT  \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 )
        // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1
        // here: a=0, b=actual temperature in Kelvin
        // \leadsto h(T) = \int_0^T c_p(T) dT
        //                 \approx 0.5 T * (cp( (0.5-0.5*\sqrt(1/3)) T) + cp((0.5+0.5*\sqrt(1/3)) T))
        //                = 0.5 T * (cp(0.2113 T) + cp(0.7887 T) )

        // enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input
        return 0.5*temperature*(spHeatCapLiquidPhase(0.2113*temperature,pressure)
                                + spHeatCapLiquidPhase(0.7887*temperature,pressure));
    }

    static Scalar boilingTemperature()
    { return 412.3; }

    template <class Evaluation>
    static Evaluation heatVap(Evaluation temperature,
                              const Evaluation& /*pressure*/)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        temperature = Toolbox::min(temperature, criticalTemperature()); // regularization
        temperature = Toolbox::max(temperature, 0.0); // regularization

        const Scalar T_crit = criticalTemperature();
        const Scalar Tr1 = boilingTemperature()/criticalTemperature();
        const Scalar p_crit = criticalPressure();

        //        Chen method, eq. 7-11.4 (at boiling)
        const Scalar DH_v_boil = Consts::R * T_crit * Tr1
            * (3.978 * Tr1 - 3.958 + 1.555*std::log(p_crit * 1e-5 /*Pa->bar*/ ) )
            / (1.07 - Tr1); /* [J/mol] */

        /* Variation with temperature according to Watson relation eq 7-12.1*/
        const Evaluation& Tr2 = temperature/criticalTemperature();
        const Scalar n = 0.375;
        const Evaluation& DH_vap = DH_v_boil * Toolbox::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);

        return (DH_vap/molarMass());          // we need [J/kg]
    }

    template <class Evaluation>
    static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
    {
        return liquidEnthalpy(temperature, pressure) + heatVap(temperature, pressure);
    }

    template <class Evaluation>
    static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
    {
        return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
    }

    template <class Evaluation>
    static Evaluation molarGasDensity(const Evaluation& temperature, const Evaluation& pressure)
    {
        return gasDensity(temperature, pressure) / molarMass();
    }

    template <class Evaluation>
    static Evaluation molarLiquidDensity(Evaluation temperature, const Evaluation& /*pressure*/)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // saturated molar volume according to Lide, CRC Handbook of
        // Thermophysical and Thermochemical Data, CRC Press, 1994
        // valid for 245 < Temperature < 600
        temperature = Toolbox::min(temperature, 500.0); // regularization
        temperature = Toolbox::max(temperature, 250.0); // regularization

        const Scalar A1 = 0.25919;           // from table
        const Scalar A2 = 0.0014569;         // from table
        const Evaluation& expo = 1.0 + Toolbox::pow((1.0 - temperature/criticalTemperature()), (2.0/7.0));
        return 1.0/(A2*Toolbox::pow(A1, expo));
    }

    template <class Evaluation>
    static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
    { return molarLiquidDensity(temperature, pressure)*molarMass(); }

    static bool gasIsCompressible()
    { return true; }

    static bool gasIsIdeal()
    { return true; }

    static bool liquidIsCompressible()
    { return false; }

    template <class Evaluation>
    static Evaluation gasViscosity(Evaluation temperature, const Evaluation& /*pressure*/)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        temperature = Toolbox::min(temperature, 500.0); // regularization
        temperature = Toolbox::max(temperature, 250.0); // regularization

        const Evaluation& Tr = Toolbox::max(temperature/criticalTemperature(), 1e-10);
        const Scalar Fp0 = 1.0;
        const Scalar xi = 0.004623;
        const Evaluation& eta_xi = Fp0*(0.807*Toolbox::pow(Tr, 0.618)
                                   - 0.357*Toolbox::exp(-0.449*Tr)
                                   + 0.34*Toolbox::exp(-4.058*Tr)
                                   + 0.018);
        return eta_xi/xi / 1e7; // [Pa s]
    }

    template <class Evaluation>
    static Evaluation liquidViscosity(Evaluation temperature, const Evaluation& /*pressure*/)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        temperature = Toolbox::min(temperature, 500.0); // regularization
        temperature = Toolbox::max(temperature, 250.0); // regularization

        const Scalar A = -3.82;
        const Scalar B = 1027.0;
        const Scalar C = -6.38e-4;
        const Scalar D = 4.52e-7;

        return 1e-3*Toolbox::exp(A
                             + B/temperature
                             + C*temperature
                             + D*temperature*temperature); // in [cP]
    }
};

} // namespace Opm

#endif