Merge pull request #1178 from chleh/mat-ads

Material properties for adsorption

CMakeLists.txt

@@ -199,6 +199,7 @@ add_subdirectory( DataHolderLib )
 add_subdirectory( FileIO )
 add_subdirectory( GeoLib )
 add_subdirectory( InSituLib )
+add_subdirectory( MaterialsLib/Adsorption )
 add_subdirectory( MathLib )
 add_subdirectory( MeshLib )
 add_subdirectory( MeshGeoToolsLib )

MaterialsLib/Adsorption/Adsorption.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <logog/include/logog.hpp>
+#include "Adsorption.h"
+
+namespace {
+    const double k_rate = 6.0e-3; // to be specified
+
+    template <typename T>
+    T square(const T& v)
+    {
+        return v * v;
+    }
+}
+
+namespace Adsorption
+{
+
+//Saturation pressure for water used in Nunez
+double AdsorptionReaction::getEquilibriumVapourPressure(const double T_Ads)
+{
+    // critical T and p
+    const double Tc = 647.3; // K
+    const double pc = 221.2e5; // Pa
+    // dimensionless T
+    const double Tr = T_Ads/Tc;
+    const double theta = 1. - Tr;
+    // empirical constants
+    const double c[] = {-7.69123,-26.08023,-168.17065,64.23285,-118.96462,4.16717,20.97506,1.0e9,6.0};
+    const double K[] = {c[0]*theta + c[1]*pow(theta,2) + c[2]*pow(theta,3) + c[3]*pow(theta,4) + c[4]*pow(theta,5),
+                        1. + c[5]*theta + c[6]*pow(theta,2)};
+
+    const double exponent = K[0]/(K[1]*Tr) - theta/(c[7]*pow(theta,2) + c[8]);
+    return pc * exp(exponent); // in Pa
+}
+
+// Evaporation enthalpy of water from Nunez
+double AdsorptionReaction::getEvaporationEnthalpy(double T_Ads) // in kJ/kg
+{
+    T_Ads -= 273.15;
+    if (T_Ads <= 10.){
+        const double c[] = {2.50052e3,-2.1068,-3.57500e-1,1.905843e-1,-5.11041e-2,7.52511e-3,-6.14313e-4,2.59674e-5,-4.421e-7};
+        double hv = 0.;
+        for (size_t i=0; i< sizeof(c)/sizeof(c[0]);i++)
+            hv += c[i] * pow(T_Ads,i);
+        return hv;
+    } else if (T_Ads <= 300.){
+        const double c[] = {2.50043e3,-2.35209,1.91685e-4,-1.94824e-5,2.89539e-7,-3.51199e-9,2.06926e-11,-6.4067e-14,8.518e-17,1.558e-20,-1.122e-22};
+        double hv = 0.;
+        for (size_t i=0; i< sizeof(c)/sizeof(c[0]);i++)
+            hv += c[i] * pow(T_Ads,i);
+        return hv;
+    } else {
+        const double c[] = {2.99866e3,-3.1837e-3,-1.566964e1,-2.514e-6,2.045933e-2,1.0389e-8};
+        return ((c[0] + c[2]*T_Ads + c[4]*pow(T_Ads,2))/(1. + c[1]*T_Ads + c[3]*pow(T_Ads,2) + c[5]*pow(T_Ads,3)));
+    }
+}
+
+
+// evaluate specific heat capacity of adsorbate follwing Nunez
+double AdsorptionReaction::getSpecificHeatCapacity(const double T_Ads)
+{
+    const double c[] = {4.224,-3.716e-3,9.351e-5,-7.1786e-7,-9.1266e-9,2.69247e-10,-2.773104e-12,1.553177e-14,-4.982795e-17,8.578e-20,-6.12423e-23};
+    double cp = 0.;
+    for (unsigned i=0; i< sizeof(c)/sizeof(c[0]);i++)
+        cp += c[i] * pow(T_Ads,i);
+    return cp; // kJ/(kg*K)
+}
+
+
+double AdsorptionReaction::getMolarFraction(double xm, double M_this, double M_other)
+{
+    return M_other*xm/(M_other*xm + M_this*(1.0-xm));
+}
+
+
+double AdsorptionReaction::getMassFraction(double xn, double M_this, double M_other)
+{
+    return M_this*xn/(M_this*xn + M_other*(1.0-xn));
+}
+
+
+double AdsorptionReaction::dMolarFraction(double xm, double M_this, double M_other)
+{
+    return M_other * M_this
+            / square(M_other * xm + M_this * (1.0 - xm));
+}
+
+
+double AdsorptionReaction::getReactionRate(const double p_Ads, const double T_Ads,
+                                     const double M_Ads, const double loading) const
+{
+    const double A = getPotential(p_Ads, T_Ads, M_Ads);
+    double C_eq = getAdsorbateDensity(T_Ads) * characteristicCurve(A);
+    if (C_eq < 0.0) C_eq = 0.0;
+
+    return k_rate * (C_eq - loading); // scaled with mass fraction
+                                      // this the rate in terms of loading!
+}
+
+void AdsorptionReaction::getDReactionRate(const double p_Ads, const double T_Ads,
+                                     const double M_Ads, const double /*loading*/,
+                                     std::array<double, 3> &dqdr) const
+{
+    const double A = getPotential(p_Ads, T_Ads, M_Ads);
+    const double p_S = getEquilibriumVapourPressure(T_Ads);
+    const double dAdT = GAS_CONST * log(p_S/p_Ads) / (M_Ads*1.e3);
+    const double dAdp = - GAS_CONST * T_Ads / M_Ads / p_Ads;
+
+    const double W = characteristicCurve(A);
+    const double dWdA = dCharacteristicCurve(A);
+
+    const double rho_Ads = getAdsorbateDensity(T_Ads);
+    const double drhodT = - rho_Ads * getAlphaT(T_Ads);
+
+    dqdr = std::array<double, 3>{{
+        rho_Ads*dWdA*dAdp,
+        drhodT*W + rho_Ads*dWdA*dAdT,
+        -k_rate
+    }};
+}
+
+
+// Evaluate adsorbtion potential A
+double AdsorptionReaction::getPotential(const double p_Ads, double T_Ads, const double M_Ads) const
+{
+    double A = GAS_CONST * T_Ads * log(getEquilibriumVapourPressure(T_Ads)/p_Ads) / (M_Ads*1.e3); // in kJ/kg = J/g
+    return A;
+}
+
+
+double AdsorptionReaction::getLoading(const double rho_curr, const double rho_dry)
+{
+    return rho_curr / rho_dry - 1.0;
+}
+
+
+// Calculate sorption entropy
+double AdsorptionReaction::getEntropy(const double T_Ads, const double A) const
+{
+    const double epsilon = 1.0e-8;
+
+    //* // This change will also change simulation results.
+    const double W_p = characteristicCurve(A+epsilon);
+    const double W_m = characteristicCurve(A-epsilon);
+    const double dAdlnW = 2.0*epsilon/(log(W_p/W_m));
+    // */
+
+    if (W_p <= 0.0 || W_m <= 0.0)
+    {
+        ERR("characteristic curve in negative region (W-, W+): %g, %g", W_m, W_p);
+        return 0.0;
+    }
+
+    return dAdlnW * getAlphaT(T_Ads);
+}
+
+
+//Calculate sorption enthalpy
+double AdsorptionReaction::getEnthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const
+{
+    // TODO [CL] consider using A as obtained from current loading (needs inverse CC A(W)) instead of p_Vapour, T_Vapour
+    const double A = getPotential(p_Ads, T_Ads, M_Ads);
+
+    return (getEvaporationEnthalpy(T_Ads) + A - T_Ads * getEntropy(T_Ads,A))*1000.0; //in J/kg
+}
+
+
+double AdsorptionReaction::getEquilibriumLoading(
+        const double p_Ads, const double T_Ads, const double M_Ads) const
+{
+    const double A = getPotential(p_Ads, T_Ads, M_Ads);
+    return getAdsorbateDensity(T_Ads) * characteristicCurve(A);
+}
+
+
+
+} // namespace Ads

MaterialsLib/Adsorption/Adsorption.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef ADSORPTION_H
+#define ADSORPTION_H
+
+#include <array>
+#include <cmath>
+
+#include "Reaction.h"
+
+namespace Adsorption
+{
+
+const double GAS_CONST = 8.3144621;
+
+const double M_N2  = 0.028;
+const double M_H2O = 0.018;
+
+class AdsorptionReaction : public Reaction
+{
+public:
+    // TODO [CL] move those three methods to water properties class
+    static double getEvaporationEnthalpy(const double T_Ads);
+    static double getEquilibriumVapourPressure(const double T_Ads);
+    static double getSpecificHeatCapacity(const double T_Ads); // TODO [CL] why unused?
+
+    static double getMolarFraction(double xm, double M_this, double M_other);
+    static double getMassFraction(double xn, double M_this, double M_other);
+    static double dMolarFraction(double xm, double M_this, double M_other);
+
+    static double getLoading(const double rho_curr, const double rho_dry);
+
+    double getEquilibriumLoading(const double p_Ads, const double T_Ads, const double M_Ads)
+    const override;
+
+    virtual double getEnthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const override;
+    virtual double getReactionRate(const double p_Ads, const double T_Ads,
+                                   const double M_Ads, const double loading) const override;
+    /**
+     * @brief get_d_reaction_rate
+     * @param p_Ads
+     * @param T_ads
+     * @param M_Ads
+     * @param loading
+     * @param dqdr array containing the differentials wrt: p, T, C
+     */
+    virtual void getDReactionRate(const double p_Ads, const double T_Ads,
+                                     const double M_Ads, const double loading,
+                                     std::array<double, 3>& dqdr) const;
+
+protected:
+    virtual double getAdsorbateDensity(const double T_Ads) const = 0;
+    virtual double getAlphaT(const double T_Ads) const = 0;
+    virtual double characteristicCurve(const double A) const = 0;
+    virtual double dCharacteristicCurve(const double A) const = 0;
+
+private:
+    double getPotential(const double p_Ads, const double T_Ads, const double M_Ads) const;
+    double getEntropy(const double T_Ads, const double A) const;
+};
+
+
+inline double curvePolyfrac(const double* coeffs, const double x)
+{
+    // TODO use Horner scheme
+    return ( coeffs[0] + coeffs[2] * x + coeffs[4] * pow(x,2) + coeffs[6] * pow(x,3) )
+            / ( 1.0 + coeffs[1] * x + coeffs[3] * pow(x,2) + coeffs[5] * pow(x,3) );
+}
+
+inline double dCurvePolyfrac(const double* coeffs, const double x)
+{
+    const double x2 = x*x;
+    const double x3 = x2*x;
+    const double u  = coeffs[0] + coeffs[2] * x +     coeffs[4] * x2 +     coeffs[6] * x3;
+    const double du =             coeffs[2]     + 2.0*coeffs[4] * x  + 3.0*coeffs[6] * x2;
+    const double v  = 1.0 + coeffs[1] * x +     coeffs[3] * x2 +     coeffs[5] * x3;
+    const double dv =       coeffs[1]     + 2.0*coeffs[3] * x  + 3.0*coeffs[5] * x2;
+
+    return (du*v - u*dv) / v / v;
+}
+
+}
+
+#endif // ADSORPTION_H

MaterialsLib/Adsorption/CMakeLists.txt

+file(GLOB Adsorption_SOURCES *.cpp)
+file(GLOB Adsorption_HEADERS *.h)
+
+add_library(MaterialsLibAdsorption ${Adsorption_HEADERS} ${Adsorption_SOURCES} )

MaterialsLib/Adsorption/Density100MPa.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "Density100MPa.h"
+
+namespace
+{
+
+// NaX_HighP_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:57
+const double c[] = {
+    0.3490302932983226,        /* a0 */
+    -0.0014061345691831226,    /* a1 */
+    -0.0007399303393402753,    /* a2 */
+    5.129318840267485e-09,    /* a3 */
+    5.243619689772646e-07,    /* a4 */
+    6.347011955956523e-10,    /* a5 */
+    -9.919599580166727e-11    /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double Density100MPa::getAdsorbateDensity(const double T_Ads) const
+{
+    return -0.0013*T_Ads*T_Ads + 0.3529*T_Ads + 1049.2;
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double Density100MPa::getAlphaT(const double T_Ads) const
+{
+    const double rho    = -0.0013*T_Ads*T_Ads+0.3529*T_Ads+1049.2;
+    const double drhodT = -0.0026*T_Ads + 0.3529;
+
+    return - drhodT / rho;
+}
+
+// Characteristic curve. Return W (A)
+double Density100MPa::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double Density100MPa::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/Density100MPa.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITY100MPA_H
+#define MATERIALSLIB_ADSORPTION_DENSITY100MPA_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class Density100MPa : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITY100MPA_H

MaterialsLib/Adsorption/DensityConst.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityConst.h"
+#include "DensityHauer.h"
+
+namespace
+{
+
+// NaX_Constant_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:20:05
+const double c[] = {
+    0.3824098506898007,            /* a0 */
+    -0.001316857559708455,        /* a1 */
+    -0.0007935756090263691,        /* a2 */
+    -1.1600036977157845e-07,    /* a3 */
+    5.610354459181838e-07,        /* a4 */
+    7.113664938298873e-10,        /* a5 */
+    -1.0668790477629686e-10        /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityConst::getAdsorbateDensity(const double /*T_Ads*/) const
+{
+    return rhoWaterHauer(150.0+273.15);
+}
+
+double DensityConst::getAlphaT(const double /*T_Ads*/) const
+{
+    return 0.0;
+}
+
+// Characteristic curve. Return W (A)
+double DensityConst::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); //cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityConst::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityConst.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYCONST_H
+#define MATERIALSLIB_ADSORPTION_DENSITYCONST_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityConst : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYCONST_H

MaterialsLib/Adsorption/DensityCook.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityCook.h"
+
+namespace
+{
+
+// NaX_Dean_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:42
+const double c[] = {
+    0.3632627555646154,        /* a0 */
+    -0.0014090624975800715,    /* a1 */
+    -0.0007717609035743321,    /* a2 */
+    5.03634836561135e-09,    /* a3 */
+    5.478509959282738e-07,    /* a4 */
+    6.36458510620815e-10,    /* a5 */
+    -1.037977321231462e-10    /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityCook::getAdsorbateDensity(const double T_Ads) const
+{
+    return rhoWaterDean(T_Ads);
+}
+
+double DensityCook::getAlphaT(const double T_Ads) const
+{
+    return alphaTWaterDean(T_Ads);
+}
+
+// Characteristic curve. Return W (A)
+double DensityCook::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); //cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; //m^3/kg
+}
+
+double DensityCook::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityCook.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYCOOK_H
+#define MATERIALSLIB_ADSORPTION_DENSITYCOOK_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityCook : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+inline double rhoWaterDean(const double T_Ads)
+{
+    const double Tcel = T_Ads - 273.15;
+    const double b[] = { 999.9,2.03E-02,-6.16E-03,2.26E-05,-4.68E-08 };
+    if (Tcel <= 100.) {
+        return b[0] + Tcel * (b[1] + Tcel * (b[2] + Tcel * (b[3] + Tcel * b[4]) ) );
+    }
+    else {
+        const double rho_100 = b[0] + b[1]*1.e2 + b[2]*1.e4 + b[3]*1.e6 + b[4]*1.e8;
+        const double aT_100  = -1./rho_100 * (b[1] + 2.*b[2]*1.e2 + 3.*b[3]*1.e4 + 4.*b[4]*1.e6);
+        return rho_100 * (1. - aT_100*(Tcel-100.));
+    }
+}
+
+inline double alphaTWaterDean(const double T_Ads)
+{
+    const double Tcel = T_Ads - 273.15;
+    const double b[] = { 999.9,2.03E-02,-6.16E-03,2.26E-05,-4.68E-08 };
+    if (Tcel <= 100.) {
+        const double r = b[0] + Tcel * (b[1] + Tcel * (b[2] + Tcel * (b[3] + Tcel * b[4]) ) );
+        return -1.0/r * ( b[1] + Tcel * (2.0*b[2] + Tcel * (3.0*b[3] + Tcel * 4.0*b[4]) ) );
+    }
+    else {
+        const double rho_100 = b[0] + b[1]*1.e2 + b[2]*1.e4 + b[3]*1.e6 + b[4]*1.e8;
+        const double aT_100  = -1./rho_100 * (b[1] + 2.*b[2]*1.e2 + 3.*b[3]*1.e4 + 4.*b[4]*1.e6);
+        return aT_100 / (1. - aT_100*(Tcel-100.));
+    }
+}
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYCOOK_H

MaterialsLib/Adsorption/DensityDubinin.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityDubinin.h"
+#include "DensityCook.h"
+#include "Adsorption.h"
+
+namespace
+{
+
+// NaX_Dubinin_polyfrac_CC.pickle
+// date extracted 2015-06-23 16:47:50 file mtime 2015-06-23 16:47:23
+const double c[] = {
+    0.3635538371322433,        /* a0 */
+    -0.0014521033261199435,    /* a1 */
+    -0.0007855160157616825,    /* a2 */
+    4.385666000850253e-08,    /* a3 */
+    5.567776459188524e-07,    /* a4 */
+    6.026002134230559e-10,    /* a5 */
+    -1.0477401124006098e-10    /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityDubinin::getAdsorbateDensity(const double T_Ads) const
+{
+    const double Tb = 373.1;
+
+    if (T_Ads < Tb) {
+        return rhoWaterDean(T_Ads);
+    } else {
+        const double Tc = 647.3; // K
+        const double pc = 221.2e5; // Pa
+        // boiling point density
+        const double rhob = rhoWaterDean(Tb);
+        // state values
+        const double R = GAS_CONST;
+        const double M = M_H2O;
+        const double b = R * Tc/(8. * pc); // m^3/mol
+        const double rhom = M/b; // kg/m^3
+        const double rho = rhob - (rhob-rhom)/(Tc-Tb)*(T_Ads-Tb);
+        return rho;
+    }
+}
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityDubinin::getAlphaT(const double T_Ads) const
+{
+    const double Tb = 373.1;
+    if (T_Ads <= Tb) {
+        return alphaTWaterDean(T_Ads);
+    } else {
+        // critical T and p
+        const double Tc = 647.3; // K
+        // const double rhoc = 322.; // kg/m^3
+        const double pc = 221.2e5; // Pa
+        // boiling point density
+        const double rhob = rhoWaterDean(Tb);
+        // state values
+        const double R = GAS_CONST;
+        const double M = M_H2O;
+        const double b = R * Tc/(8. * pc); // m^3/mol
+        const double rhom = M/(b); // kg/m^3
+        const double rho = rhob - (rhob-rhom)/(Tc-Tb)*(T_Ads-Tb);
+        return ((rhob-rhom)/(Tc-Tb)*1./rho);
+    }
+}
+
+// Characteristic curve. Return W (A)
+double DensityDubinin::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityDubinin::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityDubinin.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H
+#define MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityDubinin : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H

MaterialsLib/Adsorption/DensityHauer.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityHauer.h"
+
+namespace
+{
+
+// NaX_Hauer_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:19
+const double c[] = {
+    0.36490158988356747,    /* a0 */
+    -0.0013723270478333963,    /* a1 */
+    -0.0007655780628099964,    /* a2 */
+    -3.353324854315774e-08,    /* a3 */
+    5.424357157710913e-07,    /* a4 */
+    6.613430586648678e-10,    /* a5 */
+    -1.0300151379421499e-10    /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityHauer::getAdsorbateDensity(const double T_Ads) const
+{
+    return rhoWaterHauer(T_Ads);
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityHauer::getAlphaT(const double T_Ads) const
+{
+    // data like in python script
+    const double T0 = 283.15, alpha0 = 3.781e-4; //K; 1/K
+
+    return alpha0/(1. - alpha0 * (T_Ads-T0)); //in 1/K
+}
+
+// Characteristic curve. Return W (A)
+double DensityHauer::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityHauer::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityHauer.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYHAUER_H
+#define MATERIALSLIB_ADSORPTION_DENSITYHAUER_H
+
+#include "Adsorption.h"
+#include "DensityCook.h"
+
+namespace Adsorption
+{
+
+class DensityHauer : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+inline double rhoWaterHauer(const double T_Ads)
+{
+    // data like in python script
+    const double T0 = 283.15, rho0 = rhoWaterDean(T0), alpha0 = 3.781e-4; // K; kg/m^3; 1/K
+
+    return rho0 * (1. - alpha0 * (T_Ads-T0)); // in kg/m^3
+}
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYHAUER_H

MaterialsLib/Adsorption/DensityLegacy.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityLegacy.h"
+
+namespace
+{
+
+//parameters from least squares fit (experimental data)
+const double c[] = { 0.34102920966608297,
+                     -0.0013106032830951296,
+                     -0.00060754147575378876,
+                     3.7843404172683339e-07,
+                     4.0107503869519016e-07,
+                     3.1274595098338057e-10,
+                     -7.610441241719489e-11
+                   };
+
+}
+
+namespace Adsorption
+{
+
+double DensityLegacy::getAdsorbateDensity(const double T_Ads) const
+{
+    //set reference state for adsorbate EOS in Hauer
+    const double T0 = 293.15, rho0 = 998.084, alpha0 = 2.06508e-4; // K; kg/m^3; 1/K
+
+    return (rho0 * (1. - alpha0 * (T_Ads-T0))); // in kg/m^3
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityLegacy::getAlphaT(const double T_Ads) const
+{
+    //set reference state for adsorbate EOS in Hauer
+    const double T0 = 293.15, alpha0 = 2.06508e-4; // K; 1/K
+
+    return (alpha0/(1. - alpha0 * (T_Ads-T0))); // in 1/K
+}
+
+// Characteristic curve. Return W (A)
+double DensityLegacy::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityLegacy::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityLegacy.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H
+#define MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityLegacy : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H

MaterialsLib/Adsorption/DensityMette.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityMette.h"
+#include "DensityCook.h"
+
+namespace
+{
+// NaX_Mette_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:26
+const double c[] = {
+    0.36340572890087813,    /* a0 */
+    -0.0013449597038375108,    /* a1 */
+    -0.0007581210111121073,    /* a2 */
+    -7.331279615575401e-08,    /* a3 */
+    5.365656973806218e-07,    /* a4 */
+    6.854673678427112e-10,    /* a5 */
+    -1.0197050219481966e-10    /* a6 */
+};
+}
+
+namespace Adsorption
+{
+
+double DensityMette::getAdsorbateDensity(const double T_Ads) const
+{
+    const double T0 = 293.15;
+    const double rho0 = rhoWaterDean(T0);
+    const double alpha20 = alphaTWaterDean(T0);
+    return rho0 / (1. + alpha20*(T_Ads-T0));
+}
+
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityMette::getAlphaT(const double T_Ads) const
+{
+    const double T0 = 293.15;
+    const double alpha20 = alphaTWaterDean(T0);
+    return alpha20 / (1. + alpha20 * (T_Ads-T0));
+}
+
+
+// Characteristic curve. Return W (A)
+double DensityMette::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityMette::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityMette.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYMETTE_H
+#define MATERIALSLIB_ADSORPTION_DENSITYMETTE_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityMette : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYMETTE_H

MaterialsLib/Adsorption/DensityNunez.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityNunez.h"
+
+namespace
+{
+
+// NaX_Nunez_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:34
+const double c[] = {
+    0.3631900485031771,        /* a0 */
+    -0.0014242280940080726,    /* a1 */
+    -0.0007751726942386291,    /* a2 */
+    2.1775655036811842e-08,    /* a3 */
+    5.488166913667265e-07,    /* a4 */
+    6.204064716725214e-10,    /* a5 */
+    -1.0345385018952998e-10    /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityNunez::getAdsorbateDensity(const double T_Ads) const
+{
+    // TODO admissable T range: 273.16 K <= T_Ads <= 633.15 K
+    const double a[] = { 1.0644e3,-8.01905,1.445348e-2,-4.19589e-6,-4.5294e-9 };
+    const double b[] = { -8.039e-3,1.8698e-5,-2.3015e-8,2.3809e-11,-1.388e-14 };
+    const double u = a[0] + T_Ads * (a[1] + T_Ads * (a[2] + T_Ads * (a[3] + T_Ads * a[4]) ) );
+    const double v = 1.0 + T_Ads * (b[0] + T_Ads * (b[1] + T_Ads * (b[2] + T_Ads * (b[3] + T_Ads * b[4]) ) ) );
+    return u/v;
+}
+
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityNunez::getAlphaT(const double T_Ads) const
+{
+    // TODO admissable T range: 273.16 K <= T_Ads <= 633.15 K
+    const double a[] = { 1.0644e3,-8.01905,1.445348e-2,-4.19589e-6,-4.5294e-9 };
+    const double b[] = { -8.039e-3,1.8698e-5,-2.3015e-8,2.3809e-11,-1.388e-14 };
+    const double u = a[0] + T_Ads * (a[1] + T_Ads * (a[2] + T_Ads * (a[3] + T_Ads * a[4]) ) );
+    const double v = 1.0 + T_Ads * (b[0] + T_Ads * (b[1] + T_Ads * (b[2] + T_Ads * (b[3] + T_Ads * b[4]) ) ) );
+    const double du = a[1] + T_Ads * (2.0*a[2] + T_Ads * (3.0*a[3] + T_Ads * 4.0*a[4]) );
+    const double dv = b[0] + T_Ads * (2.0*b[1] + T_Ads * (3.0*b[2] + T_Ads * (4.0*b[3] + T_Ads * 5.0*b[4]) ) );
+    return dv/v - du/u;
+}
+
+
+// Characteristic curve. Return W (A)
+double DensityNunez::characteristicCurve(const double A) const
+{
+    double W = curvePolyfrac(c, A); // cm^3/g
+
+    if (W < 0.0) {
+        W = 0.0; // TODO [CL] debug output
+    }
+
+    return W/1.e3; // m^3/kg
+}
+
+double DensityNunez::dCharacteristicCurve(const double A) const
+{
+    return dCurvePolyfrac(c, A);
+}
+
+}

MaterialsLib/Adsorption/DensityNunez.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H
+#define MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityNunez : public AdsorptionReaction
+{
+public:
+    double getAdsorbateDensity(const double T_Ads) const;
+    double getAlphaT(const double T_Ads) const;
+    double characteristicCurve(const double A) const;
+    double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H