diff --git a/MaterialsLib/adsorption/CMakeLists.txt b/MaterialsLib/adsorption/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5889ae34c4c433e862f53b023545e258a33eb498
--- /dev/null
+++ b/MaterialsLib/adsorption/CMakeLists.txt
@@ -0,0 +1,23 @@
+file(GLOB adsorption_SOURCES *.cpp)
+file(GLOB adsorption_HEADERS *.h)
+
+# set( HEADERS
+# adsorption.h
+# density_legacy.h
+# density_hauer.h
+# density_mette.h
+# )
+#
+# set( SOURCES
+# adsorption.cpp
+# density_legacy.cpp
+# density_hauer.cpp
+# density_mette.cpp
+# )
+
+
+# include_directories(
+# ${CMAKE_SOURCE_DIR}/FEM
+# )
+
+add_library(MaterialsLib_adsorption STATIC ${adsorption_HEADERS} ${adsorption_SOURCES} )
diff --git a/MaterialsLib/adsorption/adsorption.cpp b/MaterialsLib/adsorption/adsorption.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5c19af34bcbedd45f0044ebb89e92ae4ce1e57e0
--- /dev/null
+++ b/MaterialsLib/adsorption/adsorption.cpp
@@ -0,0 +1,262 @@
+
+#include "logog/include/logog.hpp"
+
+#include "adsorption.h"
+
+#include "density_legacy.h"
+#include "density_100MPa.h"
+#include "density_const.h"
+#include "density_cook.h"
+#include "density_dubinin.h"
+#include "density_hauer.h"
+#include "density_mette.h"
+#include "density_nunez.h"
+#include "density_inert.h"
+
+namespace {
+ const double k_rate = 6.0e-3; //to be specified
+
+ template <typename T>
+ T square(const T& v)
+ {
+ return v * v;
+ }
+}
+
+namespace Ads
+{
+
+//Saturation pressure for water used in Nunez
+double Adsorption::get_equilibrium_vapour_pressure(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 Adsorption::get_evaporation_enthalpy(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 Adsorption::get_specific_heat_capacity(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 Adsorption::get_molar_fraction(double xm, double M_this, double M_other)
+{
+ return M_other*xm/(M_other*xm + M_this*(1.0-xm));
+}
+
+
+double Adsorption::get_mass_fraction(double xn, double M_this, double M_other)
+{
+ return M_this*xn/(M_this*xn + M_other*(1.0-xn));
+}
+
+
+double Adsorption::d_molar_fraction(double xm, double M_this, double M_other)
+{
+ return M_other * M_this
+ / square(M_other * xm + M_this * (1.0 - xm));
+}
+
+
+double Adsorption::get_reaction_rate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading) const
+{
+ // const double k_rate = 3.0e-3; //to be specified
+
+ const double A = get_potential(p_Ads, T_Ads, M_Ads);
+ double C_eq = get_adsorbate_density(T_Ads) * characteristic_curve(A);
+ if (C_eq < 0.0) C_eq = 0.0;
+
+ // return 0.0; // TODO [CL] for testing only
+
+ return k_rate * (C_eq - loading); //scaled with mass fraction
+ // this the rate in terms of loading!
+}
+
+void Adsorption::get_d_reaction_rate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double /*loading*/,
+ std::array<double, 3> &dqdr) const
+{
+ const double A = get_potential(p_Ads, T_Ads, M_Ads);
+ const double p_S = get_equilibrium_vapour_pressure(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 = characteristic_curve(A);
+ const double dWdA = d_characteristic_curve(A);
+
+ const double rho_Ads = get_adsorbate_density(T_Ads);
+ const double drhodT = - rho_Ads * get_alphaT(T_Ads);
+
+ dqdr = std::array<double, 3>{{
+ rho_Ads*dWdA*dAdp,
+ drhodT*W + rho_Ads*dWdA*dAdT,
+ -k_rate
+ }};
+}
+
+
+//Evaluate adsorbtion potential A
+double Adsorption::get_potential(const double p_Ads, double T_Ads, const double M_Ads) const
+{
+ double A = GAS_CONST * T_Ads * log(get_equilibrium_vapour_pressure(T_Ads)/p_Ads) / (M_Ads*1.e3); //in kJ/kg = J/g
+ if (A < 0.0) {
+ // vapour partial pressure > saturation pressure
+ // A = 0.0; // TODO [CL] debug output
+ }
+ return A;
+}
+
+
+double Adsorption::get_loading(const double rho_curr, const double rho_dry)
+{
+ return rho_curr / rho_dry - 1.0;
+}
+
+
+//Calculate sorption entropy
+double Adsorption::get_entropy(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 = characteristic_curve(A+epsilon);
+ const double W_m = characteristic_curve(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;
+ }
+
+ // const double dAdlnW = 2.0*epsilon/(log(characteristic_curve(A+epsilon)) - log(characteristic_curve(A-epsilon)));
+ return dAdlnW * get_alphaT(T_Ads);
+}
+
+
+//Calculate sorption enthalpy
+double Adsorption::get_enthalpy(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 = get_potential(p_Ads, T_Ads, M_Ads);
+
+ // return 0.0; // TODO [CL] for testing only
+ return (get_evaporation_enthalpy(T_Ads) + A - T_Ads * get_entropy(T_Ads,A))*1000.0; //in J/kg
+}
+
+
+double Adsorption::get_equilibrium_loading(const double p_Ads, const double T_Ads, const double M_Ads)
+{
+ const double A = get_potential(p_Ads, T_Ads, M_Ads);
+ return get_adsorbate_density(T_Ads) * characteristic_curve(A);
+}
+
+
+bool
+stringToReactiveSystem(std::string const& name, SolidReactiveSystem& rsys_out)
+{
+ if (name == "Z13XBF")
+ rsys_out = SolidReactiveSystem::Z13XBF;
+ else if (name == "Z13XBF_100MPa")
+ rsys_out = SolidReactiveSystem::Z13XBF_100MPa;
+ else if (name == "Z13XBF_Const")
+ rsys_out = SolidReactiveSystem::Z13XBF_Const;
+ else if (name == "Z13XBF_Cook")
+ rsys_out = SolidReactiveSystem::Z13XBF_Cook;
+ else if (name == "Z13XBF_Dubinin")
+ rsys_out = SolidReactiveSystem::Z13XBF_Dubinin;
+ else if (name == "Z13XBF_Hauer")
+ rsys_out = SolidReactiveSystem::Z13XBF_Hauer;
+ else if (name == "Z13XBF_Mette")
+ rsys_out = SolidReactiveSystem::Z13XBF_Mette;
+ else if (name == "Z13XBF_Nunez")
+ rsys_out = SolidReactiveSystem::Z13XBF_Nunez;
+ else if (name == "Inert")
+ rsys_out = SolidReactiveSystem::Inert;
+ else return false;
+
+ return true;
+}
+
+
+Adsorption* Adsorption::newInstance(std::string const& rsys)
+{
+ SolidReactiveSystem r;
+ if (stringToReactiveSystem(rsys, r)) {
+ return newInstance(r);
+ } else {
+ ERR("unknown reactive system: %s", rsys);
+ return nullptr;
+ }
+}
+
+
+Adsorption* Adsorption::newInstance(const SolidReactiveSystem rsys)
+{
+ switch (rsys)
+ {
+ case SolidReactiveSystem::Z13XBF:
+ return new DensityLegacy();
+ case SolidReactiveSystem::Z13XBF_100MPa:
+ return new Density100MPa();
+ case SolidReactiveSystem::Z13XBF_Const:
+ return new DensityConst();
+ case SolidReactiveSystem::Z13XBF_Cook:
+ return new DensityCook();
+ case SolidReactiveSystem::Z13XBF_Dubinin:
+ return new DensityDubinin();
+ case SolidReactiveSystem::Z13XBF_Hauer:
+ return new DensityHauer();
+ case SolidReactiveSystem::Z13XBF_Mette:
+ return new DensityMette();
+ case SolidReactiveSystem::Z13XBF_Nunez:
+ return new DensityNunez();
+ case SolidReactiveSystem::Inert:
+ return new DensityInert();
+ }
+
+ return nullptr;
+}
+
+} // namespace Ads
diff --git a/MaterialsLib/adsorption/adsorption.h b/MaterialsLib/adsorption/adsorption.h
new file mode 100644
index 0000000000000000000000000000000000000000..8c953ca58061fe02fafa39e8f0cf23a2ba12dcc8
--- /dev/null
+++ b/MaterialsLib/adsorption/adsorption.h
@@ -0,0 +1,112 @@
+#ifndef ADSORPTION_H
+#define ADSORPTION_H
+
+#include <array>
+#include <cmath>
+
+namespace Ads
+{
+
+const double GAS_CONST = 8.3144621;
+
+const double M_N2 = 0.028;
+const double M_H2O = 0.018;
+
+enum class SolidReactiveSystem
+{
+ Z13XBF,
+ Z13XBF_Const,
+ Z13XBF_Hauer,
+ Z13XBF_Mette,
+ Z13XBF_Nunez,
+ Z13XBF_Cook,
+ Z13XBF_Dubinin,
+ Z13XBF_100MPa,
+ Inert
+};
+
+class Adsorption
+{
+public:
+// static members
+ static Adsorption* newInstance(SolidReactiveSystem rsys);
+ static Adsorption* newInstance(std::string const& rsys);
+
+ // TODO [CL] move those four methods to water properties class
+ static double get_evaporation_enthalpy(const double T_Ads);
+ static double get_equilibrium_vapour_pressure(const double T_Ads);
+ static double get_specific_heat_capacity(const double T_Ads); // TODO [CL] why unused?
+
+ static double get_molar_fraction(double xm, double M_this, double M_other);
+ static double get_mass_fraction(double xn, double M_this, double M_other);
+ static double d_molar_fraction(double xm, double M_this, double M_other);
+
+ static double get_loading(const double rho_curr, const double rho_dry);
+
+// non-virtual members
+ double get_equilibrium_loading(const double p_Ads, const double T_Ads, const double M_Ads);
+
+// virtual members:
+ virtual ~Adsorption() {}
+
+ virtual double get_enthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const;
+ virtual double get_reaction_rate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading) const;
+ /**
+ * @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 get_d_reaction_rate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading,
+ std::array<double, 3>& dqdr) const;
+
+protected:
+ virtual double get_adsorbate_density(const double T_Ads) const = 0;
+ virtual double get_alphaT(const double T_Ads) const = 0;
+ virtual double characteristic_curve(const double A) const = 0;
+ virtual double d_characteristic_curve(const double A) const = 0;
+
+private:
+// non-virtual members
+ double get_potential(const double p_Ads, const double T_Ads, const double M_Ads) const;
+ double get_entropy(const double T_Ads, const double A) const;
+};
+
+
+inline double curve_polyfrac(const double* coeffs, const double x)
+{
+ 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) );
+
+ // Apparently, even pow and std::pow are different
+ // return ( coeffs[0] + coeffs[2] * x + coeffs[4] * std::pow(x,2) + coeffs[6] * std::pow(x,3) )
+ // / ( 1.0 + coeffs[1] * x + coeffs[3] * std::pow(x,2) + coeffs[5] * std::pow(x,3) );
+
+ // Analytically the same, but numerically quite different
+ // return ( coeffs[0] + x * ( coeffs[2] + x * (coeffs[4] + x * coeffs[6] ) ) )
+ // / ( 1.0 + x * ( coeffs[1] + x * (coeffs[3] + x * coeffs[5] ) ) );
+
+ // Analytically the same, but numerically quite different
+ // return ( coeffs[0] + x * coeffs[2] + x*x * coeffs[4] + x*x*x * coeffs[6] )
+ // / ( 1.0 + x * coeffs[1] + x*x * coeffs[3] + x*x*x * coeffs[5] );
+}
+
+inline double d_curve_polyfrac(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
diff --git a/MaterialsLib/adsorption/density_100MPa.cpp b/MaterialsLib/adsorption/density_100MPa.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f49eb1f221beb7b3b85ed66ca80d91ac3f40e656
--- /dev/null
+++ b/MaterialsLib/adsorption/density_100MPa.cpp
@@ -0,0 +1,56 @@
+#include "density_100MPa.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 Ads
+{
+
+double Density100MPa::get_adsorbate_density(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::get_alphaT(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::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_100MPa.h b/MaterialsLib/adsorption/density_100MPa.h
new file mode 100644
index 0000000000000000000000000000000000000000..54ae755290bc7a8ecb29616d5c62c124eaf7937c
--- /dev/null
+++ b/MaterialsLib/adsorption/density_100MPa.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class Density100MPa : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_const.cpp b/MaterialsLib/adsorption/density_const.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..29164e7fc72492f5d9b0d9a4b238936deeb94680
--- /dev/null
+++ b/MaterialsLib/adsorption/density_const.cpp
@@ -0,0 +1,55 @@
+#include "density_const.h"
+
+#include "density_hauer.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 Ads
+{
+
+double DensityConst::get_adsorbate_density(const double /*T_Ads*/) const
+{
+ return rho_water_Hauer(150.0+273.15);
+}
+
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityConst::get_alphaT(const double /*T_Ads*/) const
+{
+ return 0.0;
+}
+
+
+//Characteristic curve. Return W (A)
+double DensityConst::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_const.h b/MaterialsLib/adsorption/density_const.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d64ec9df4993415ec0d589b096bcda180445bd3
--- /dev/null
+++ b/MaterialsLib/adsorption/density_const.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityConst : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_cook.cpp b/MaterialsLib/adsorption/density_cook.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..96aeb959ae858343b38d19564584f11fcb0fda82
--- /dev/null
+++ b/MaterialsLib/adsorption/density_cook.cpp
@@ -0,0 +1,53 @@
+#include "density_cook.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 Ads
+{
+
+double DensityCook::get_adsorbate_density(const double T_Ads) const
+{
+ return rho_water_Dean(T_Ads);
+}
+
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityCook::get_alphaT(const double T_Ads) const
+{
+ return alphaT_water_Dean(T_Ads);
+}
+
+
+//Characteristic curve. Return W (A)
+double DensityCook::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_cook.h b/MaterialsLib/adsorption/density_cook.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ef429833c430c7fe3c1e031f0269c977cd7e7b0
--- /dev/null
+++ b/MaterialsLib/adsorption/density_cook.h
@@ -0,0 +1,47 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityCook : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+
+inline double rho_water_Dean(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 alphaT_water_Dean(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.));
+ }
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_dubinin.cpp b/MaterialsLib/adsorption/density_dubinin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..91a5248f026490f12bde98d18d385f53402da1ce
--- /dev/null
+++ b/MaterialsLib/adsorption/density_dubinin.cpp
@@ -0,0 +1,89 @@
+#include "density_dubinin.h"
+#include "density_cook.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 Ads
+{
+
+double DensityDubinin::get_adsorbate_density(const double T_Ads) const
+{
+ const double Tb = 373.1;
+
+ if (T_Ads < Tb) {
+ return rho_water_Dean(T_Ads);
+ } else {
+ 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 = rho_water_Dean(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::get_alphaT(const double T_Ads) const
+{
+ const double Tb = 373.1;
+ if (T_Ads <= Tb) {
+ return alphaT_water_Dean(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 = rho_water_Dean(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::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_dubinin.h b/MaterialsLib/adsorption/density_dubinin.h
new file mode 100644
index 0000000000000000000000000000000000000000..678aea58c167702dd1d874c9a9b624b95edc6af7
--- /dev/null
+++ b/MaterialsLib/adsorption/density_dubinin.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityDubinin : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_hauer.cpp b/MaterialsLib/adsorption/density_hauer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e9211f552545cd40964ab572a71686aaae004409
--- /dev/null
+++ b/MaterialsLib/adsorption/density_hauer.cpp
@@ -0,0 +1,56 @@
+#include "density_hauer.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 Ads
+{
+
+double DensityHauer::get_adsorbate_density(const double T_Ads) const
+{
+ return rho_water_Hauer(T_Ads);
+}
+
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityHauer::get_alphaT(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::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_hauer.h b/MaterialsLib/adsorption/density_hauer.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d65991d0f129f8f302e0615ecebd34f39a20536
--- /dev/null
+++ b/MaterialsLib/adsorption/density_hauer.h
@@ -0,0 +1,26 @@
+#pragma once
+
+#include "adsorption.h"
+#include "density_cook.h"
+
+namespace Ads
+{
+
+class DensityHauer : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+inline double rho_water_Hauer(const double T_Ads)
+{
+ // data like in python script
+ const double T0 = 283.15, rho0 = rho_water_Dean(T0), alpha0 = 3.781e-4; //K; kg/m^3; 1/K
+
+ return rho0 * (1. - alpha0 * (T_Ads-T0)); //in kg/m^3
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_inert.h b/MaterialsLib/adsorption/density_inert.h
new file mode 100644
index 0000000000000000000000000000000000000000..e16f9c93b79c96772f52676f345e6d16f2fee779
--- /dev/null
+++ b/MaterialsLib/adsorption/density_inert.h
@@ -0,0 +1,42 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityInert final : public Adsorption
+{
+public:
+ double get_enthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+ const double /*M_Ads*/) const override
+ {
+ return 0.0;
+ }
+ double get_reaction_rate(const double /*p_Ads*/, const double /*T_Ads*/,
+ const double /*M_Ads*/, const double /*loading*/
+ ) const override
+ {
+ return 0.0;
+ }
+
+protected:
+ double get_adsorbate_density(const double /*T_Ads*/) const override
+ {
+ return 1.0;
+ }
+ double get_alphaT(const double /*T_Ads*/) const override
+ {
+ return 0.0;
+ }
+ double characteristic_curve(const double A) const override
+ {
+ return A;
+ }
+ double d_characteristic_curve(const double /*A*/) const override
+ {
+ return 1.0;
+ }
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_legacy.cpp b/MaterialsLib/adsorption/density_legacy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e8b07e86f2c6c7b9ff6b2c1e1b78c72986fa181c
--- /dev/null
+++ b/MaterialsLib/adsorption/density_legacy.cpp
@@ -0,0 +1,59 @@
+#include "density_legacy.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 Ads
+{
+
+double DensityLegacy::get_adsorbate_density(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::get_alphaT(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 (alpha0/(1. - alpha0 * (T_Ads-T0))); //in 1/K
+}
+
+
+//Characteristic curve. Return W (A)
+double DensityLegacy::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_legacy.h b/MaterialsLib/adsorption/density_legacy.h
new file mode 100644
index 0000000000000000000000000000000000000000..d9123fc6a1a293ed24d68676b07cbcf32f2784de
--- /dev/null
+++ b/MaterialsLib/adsorption/density_legacy.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityLegacy : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_mette.cpp b/MaterialsLib/adsorption/density_mette.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7cf6f36e29887a71781e0990ed5c82e8fa06b4cb
--- /dev/null
+++ b/MaterialsLib/adsorption/density_mette.cpp
@@ -0,0 +1,57 @@
+#include "density_mette.h"
+#include "density_cook.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 Ads
+{
+
+double DensityMette::get_adsorbate_density(const double T_Ads) const
+{
+ const double T0 = 293.15;
+ const double rho0 = rho_water_Dean(T0);
+ const double alpha20 = alphaT_water_Dean(T0);
+ return rho0 / (1. + alpha20*(T_Ads-T0));
+}
+
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityMette::get_alphaT(const double T_Ads) const
+{
+ const double T0 = 293.15;
+ const double alpha20 = alphaT_water_Dean(T0);
+ return alpha20 / (1. + alpha20 * (T_Ads-T0));
+}
+
+
+//Characteristic curve. Return W (A)
+double DensityMette::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_mette.h b/MaterialsLib/adsorption/density_mette.h
new file mode 100644
index 0000000000000000000000000000000000000000..b66b8ca214905ac7a70f6816a7515051a8e84701
--- /dev/null
+++ b/MaterialsLib/adsorption/density_mette.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityMette : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}
diff --git a/MaterialsLib/adsorption/density_nunez.cpp b/MaterialsLib/adsorption/density_nunez.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac81d53f6973bead70e502fc226a81d920b7955f
--- /dev/null
+++ b/MaterialsLib/adsorption/density_nunez.cpp
@@ -0,0 +1,75 @@
+#include "density_nunez.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 Ads
+{
+
+double DensityNunez::get_adsorbate_density(const double T_Ads) const
+{
+ if (T_Ads < 273.16 or T_Ads > 633.15) {
+ // print('Value outside admissible range for rho.');
+ // return -1;
+ }
+
+ 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::get_alphaT(const double T_Ads) const
+{
+ if (T_Ads < 273.16 or T_Ads > 633.15) {
+ // print('Value outside admissible range for rho.');
+ // return -1;
+ }
+
+ 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 - (du*v - dv*u) / u / v;
+ // return (dv*u - du*v) / u / v;
+ return dv/v - du/u;
+}
+
+
+//Characteristic curve. Return W (A)
+double DensityNunez::characteristic_curve(const double A) const
+{
+ double W = curve_polyfrac(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::d_characteristic_curve(const double A) const
+{
+ return d_curve_polyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/adsorption/density_nunez.h b/MaterialsLib/adsorption/density_nunez.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d05d83feba29a514388eedbf8da4ebcc33bcea2
--- /dev/null
+++ b/MaterialsLib/adsorption/density_nunez.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+class DensityNunez : public Adsorption
+{
+public:
+ double get_adsorbate_density(const double T_Ads) const;
+ double get_alphaT(const double T_Ads) const;
+ double characteristic_curve(const double A) const;
+ double d_characteristic_curve(const double A) const;
+};
+
+}