[Mat] changes to adsorption

MaterialsLib/adsorption/adsorption.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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"
 
-#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
 
@@ -186,77 +183,13 @@ double Adsorption::get_enthalpy(const double p_Ads, const double T_Ads, const do
 }
 
 
-double Adsorption::get_equilibrium_loading(const double p_Ads, const double T_Ads, const double M_Ads)
+double Adsorption::get_equilibrium_loading(const double p_Ads, const double T_Ads,
+										   const double M_Ads) const
 {
 	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

MaterialsLib/adsorption/adsorption.h

@@ -4,35 +4,20 @@
 #include <array>
 #include <cmath>
 
+#include "reaction.h"
+
 namespace Ads
 {
 
-const double GAS_CONST = 8.3144621;
+constexpr double GAS_CONST = 8.3144621;
 
-const double M_N2  = 0.028;
-const double M_H2O = 0.018;
+constexpr double M_N2  = 0.028;
+constexpr 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
+class Adsorption : public Reaction
 {
 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
+	// TODO [CL] move those three 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?
@@ -44,14 +29,15 @@ public:
 	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);
+	double get_equilibrium_loading(const double p_Ads, const double T_Ads, const double M_Ads)
+	const override;
 
 // virtual members:
-	virtual ~Adsorption() {}
+	virtual ~Adsorption() = default;
 
-	virtual double get_enthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const;
+	virtual double get_enthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const override;
 	virtual double get_reaction_rate(const double p_Ads, const double T_Ads,
-									 const double M_Ads, const double loading) const;
+									 const double M_Ads, const double loading) const override;
 	/**
 	 * @brief get_d_reaction_rate
 	 * @param p_Ads

MaterialsLib/adsorption/reaction.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 "reaction.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 "reaction_inert.h"
+#include "reaction_sinusoidal.h"
+#include "reaction_CaOH2.h"
+
+
+namespace Ads
+{
+
+std::unique_ptr<Reaction>
+Reaction::
+newInstance(BaseLib::ConfigTreeNew const& conf)
+{
+	auto const type = conf.getConfParam<std::string>("type");
+
+	if (type == "Z13XBF")
+		return std::unique_ptr<Reaction>(new DensityLegacy);
+	else if (type == "Z13XBF_100MPa")
+		return std::unique_ptr<Reaction>(new Density100MPa);
+	else if (type == "Z13XBF_Const")
+		return std::unique_ptr<Reaction>(new DensityConst);
+	else if (type == "Z13XBF_Cook")
+		return std::unique_ptr<Reaction>(new DensityCook);
+	else if (type == "Z13XBF_Dubinin")
+		return std::unique_ptr<Reaction>(new DensityDubinin);
+	else if (type == "Z13XBF_Hauer")
+		return std::unique_ptr<Reaction>(new DensityHauer);
+	else if (type == "Z13XBF_Mette")
+		return std::unique_ptr<Reaction>(new DensityMette);
+	else if (type == "Z13XBF_Nunez")
+		return std::unique_ptr<Reaction>(new DensityNunez);
+	else if (type == "Inert")
+		return std::unique_ptr<Reaction>(new ReactionInert);
+	else if (type == "Sinusoidal")
+		return std::unique_ptr<Reaction>(new ReactionSinusoidal(conf));
+	else if (type == "CaOH2")
+		return std::unique_ptr<Reaction>(new ReactionCaOH2(conf));
+
+	ERR("Unknown reactive system: %s.", type.c_str());
+	std::abort();
+
+	return nullptr;
+}
+
+} // namespace Ads

MaterialsLib/adsorption/reaction.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include<memory>
+
+#include "BaseLib/ConfigTreeNew.h"
+
+namespace Ads
+{
+
+class Reaction
+{
+public:
+	static std::unique_ptr<Reaction> newInstance(BaseLib::ConfigTreeNew const& rsys);
+
+	virtual double get_enthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const = 0;
+	virtual double get_reaction_rate(const double p_Ads, const double T_Ads,
+									 const double M_Ads, const double loading) const = 0;
+
+	// TODO: get rid of
+	virtual double get_equilibrium_loading(const double, const double, const double) const {
+		return -1.0;
+	}
+
+    virtual ~Reaction() = default;
+};
+
+}

MaterialsLib/adsorption/reaction_CaOH2.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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<cassert>
+
+#include <logog/include/logog.hpp>
+
+#include "reaction_CaOH2.h"
+
+#include "adsorption.h"
+
+namespace Ads
+{
+
+constexpr double ReactionCaOH2::tol_l;
+constexpr double ReactionCaOH2::tol_u;
+
+
+double
+ReactionCaOH2::get_enthalpy(const double, const double, const double) const
+{
+    return - reaction_enthalpy/M_react;
+}
+
+double
+ReactionCaOH2::get_reaction_rate(const double, const double, const double, const double) const
+{
+    ERR("get_reaction_rate do not call directly");
+    std::abort();
+    // TODO: error
+    return -1.0;
+}
+
+
+void ReactionCaOH2::eval(double /*t*/,
+						 const BaseLib::ArrayRef<const double, 1>& y,
+						 BaseLib::ArrayRef<double, 1>& dydx)
+{
+	assert( y.size() == dydx.size() );
+
+	rho_s = y[0];
+	calculate_qR();
+	dydx[0] = qR;
+}
+
+void ReactionCaOH2::update_param(
+		double T_solid,
+		double p_gas,
+		double x_react,
+		double rho_s_initial)
+{
+	T_s     = T_solid;
+	this->p_gas = p_gas / 1e5; // convert Pa to bar
+	this->x_react = x_react;
+	rho_s   = rho_s_initial;
+}
+
+void ReactionCaOH2::calculate_qR()
+{
+	//Convert mass fraction into mole fraction
+	// const double mol_frac_react = get_mole_fraction(x_react);
+	const double mol_frac_react = Ads::Adsorption::get_molar_fraction(x_react, M_react, M_carrier);
+
+	p_r_g = std::max(mol_frac_react * p_gas, 1.0e-3); //avoid illdefined log
+	set_chemical_equilibrium();
+	const double dXdt = Ca_hydration();
+	qR = (rho_up - rho_low) * dXdt;
+}
+
+//determine equilibrium temperature and pressure according to van't Hoff
+void ReactionCaOH2::set_chemical_equilibrium()
+{
+	X_D = (rho_s - rho_up - tol_rho)/(rho_low - rho_up - 2.0*tol_rho) ;
+	X_D = (X_D < 0.5) ? std::max(tol_l,X_D) : std::min(X_D,tol_u); //constrain to interval [tol_l;tol_u]
+
+	X_H = 1.0 - X_D;
+
+	//calculate equilibrium
+	// using the p_eq to calculate the T_eq - Clausius-Clapeyron
+	T_eq = (reaction_enthalpy/R) / ((reaction_entropy/R) + std::log(p_r_g)); // unit of p in bar
+	//Alternative: Use T_s as T_eq and calculate p_eq - for Schaube kinetics
+	p_eq = exp((reaction_enthalpy/R)/T_s - (reaction_entropy/R));
+}
+
+
+double ReactionCaOH2::Ca_hydration()
+{
+	double dXdt;
+		// step 3, calculate dX/dt
+#ifdef SIMPLE_KINETICS
+	if ( T_s < T_eq ) // hydration - simple model
+#else
+	if ( p_r_g > p_eq ) // hydration - Schaube model
+#endif
+	{
+		//X_H = max(tol_l,X_H); //lower tolerance to avoid oscillations at onset of hydration reaction. Set here so that no residual reaction rate occurs at end of hydration.
+#ifdef SIMPLE_KINETICS // this is from P. Schmidt
+		dXdt = -1.0*(1.0-X_H) * (T_s - T_eq) / T_eq * 0.2 * conversion_rate::x_react;
+#else //this is from Schaube
+		if (X_H == tol_u || rho_s == rho_up)
+			dXdt = 0.0;
+		else if ( (T_eq-T_s) >= 50.0)
+			dXdt = 13945.0 * exp(-89486.0/R/T_s) * std::pow(p_r_g/p_eq - 1.0,0.83) * 3.0 * (X_D) * std::pow(-1.0*log(X_D),0.666);
+		else
+			dXdt = 1.0004e-34 * exp(5.3332e4/T_s) * std::pow(p_r_g, 6.0) * (X_D);
+#endif
+	}
+	else // dehydration
+	{
+		//X_D = max(tol_l,X_D); //lower tolerance to avoid oscillations at onset of dehydration reaction. Set here so that no residual reaction rate occurs at end of dehydration.
+#ifdef SIMPLE_KINETICS // this is from P. Schmidt
+		dXdt = -1.0* (1.0-X_D) * (T_s - T_eq) / T_eq * 0.05;
+#else
+		if (X_D == tol_u || rho_s == rho_low)
+			dXdt = 0.0;
+		else if (X_D < 0.2)
+			dXdt = -1.9425e12 * exp( -1.8788e5/R/T_s ) * std::pow(1.0-p_r_g/p_eq,3.0)*(X_H);
+		else
+			dXdt = -8.9588e9 * exp( -1.6262e5/R/T_s ) * std::pow(1.0-p_r_g/p_eq,3.0)*2.0 * std::pow(X_H, 0.5);
+#endif
+	}
+	return dXdt;
+}
+
+
+
+
+
+
+
+
+
+}

MaterialsLib/adsorption/reaction_CaOH2.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/ArrayRef.h"
+
+#include "reaction.h"
+#include "adsorption.h"
+
+
+namespace ProcessLib
+{
+template<typename>
+class TESFEMReactionAdaptorCaOH2;
+}
+
+namespace Ads
+{
+
+class ReactionCaOH2 final : public Reaction
+{
+public:
+    explicit ReactionCaOH2(BaseLib::ConfigTreeNew const& conf)
+        : ode_solver_config{conf.getConfSubtree("ode_solver_config")}
+    {
+        /*auto const param = conf.get_optional<double>("reaction_enthalpy");
+        if (param) {
+            _enthalpy = *param;
+        } else {
+            ERR("<reaction_enthalpy> not specified.");
+            std::abort();
+        }*/
+    }
+
+    double get_enthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+                        const double /*M_Ads*/) const override;
+
+    double get_reaction_rate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
+                             const double /*loading*/) const override;
+
+    const BaseLib::ConfigTreeNew& getOdeSolverConfig() const { return ode_solver_config; }
+
+
+    void eval(double /*t*/,
+              BaseLib::ArrayRef<const double, 1> const& y,
+              BaseLib::ArrayRef<double, 1>& dydx);
+    void update_param(double T_solid,
+                      double p_gas,
+                      double x_react,
+                      double rho_s_initial);
+
+private:
+    void calculate_qR();
+    void set_chemical_equilibrium();
+    double Ca_hydration();
+
+    static constexpr double R = Ads::GAS_CONST;  // [J/mol/K]
+    double rho_s;                 // solid phase density
+    double p_gas;                 // gas phase pressure in unit bar
+    double p_r_g;                 // pressure of H2O on gas phase;
+    double p_eq       = 1.0;      // equilibrium pressure; // unit in bar
+    double T_eq;                  // equilibrium temperature;
+    double T_s;                   // solid phase temperature;
+    double qR;                    // rate of solid density change;
+    double x_react;               // mass fraction of water in gas phase;
+    double X_D;                   // mass fraction of dehydration (CaO) in the solid phase;
+    double X_H;                   // mass fraction of hydration in the solid phase;
+    static constexpr double reaction_enthalpy = -1.12e+05; // in J/mol; negative for exothermic composition reaction
+    static constexpr double reaction_entropy  = -143.5;    // in J/mol/K
+    static constexpr double M_carrier = Ads::M_N2;         // inert component molar mass
+    static constexpr double M_react   = Ads::M_H2O;        // reactive component molar mass
+
+    static constexpr double tol_l   = 1e-4;
+    static constexpr double tol_u   = 1.0 - 1e-4;
+    static constexpr double tol_rho = 0.1;
+
+    const BaseLib::ConfigTreeNew ode_solver_config;
+
+    template<typename>
+    friend class ProcessLib::TESFEMReactionAdaptorCaOH2;
+
+public:
+    static constexpr double rho_low = 1656.0;              // lower density limit
+    static constexpr double rho_up = 2200.0;               // upper density limit
+};
+
+}

MaterialsLib/adsorption/reaction_inert.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "reaction.h"
+
+namespace Ads
+{
+
+class ReactionInert final : public Reaction
+{
+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
+    {
+        ERR("Method get_reaction_rate() should never be called directly");
+        std::abort();
+        return 0.0;
+    }
+};
+
+}

MaterialsLib/adsorption/reaction_sinusoidal.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <logog/include/logog.hpp>
+
+#include "reaction.h"
+#include "BaseLib/ConfigTreeNew.h"
+
+namespace Ads
+{
+
+class ReactionSinusoidal final : public Reaction
+{
+public:
+    explicit ReactionSinusoidal(BaseLib::ConfigTreeNew const& conf)
+        : _enthalpy(conf.getConfParam<double>("reaction_enthalpy"))
+    {
+    }
+
+    double get_enthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+                        const double /*M_Ads*/) const override
+    {
+        return _enthalpy;
+    }
+
+    double get_reaction_rate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
+                             const double /*loading*/) const override
+    {
+        ERR("Method get_reaction_rate() should never be called directly");
+        std::abort();
+        return 0.0;
+    }
+
+private:
+    double _enthalpy;
+};
+
+}