Skip to content

Commit

Permalink
[Ads] removed dead code, reformatted comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@chleh
chleh committed May 24, 2016
1 parent 6ed7c2b commit bb45bee
Show file tree
Hide file tree
Showing 16 changed files with 81 additions and 144 deletions.
44 changes: 17 additions & 27 deletions MaterialsLib/Adsorption/Adsorption.cpp
View file
Expand Up @@ -11,7 +11,7 @@
#include "Adsorption.h"

namespace {
const double k_rate = 6.0e-3; //to be specified
const double k_rate = 6.0e-3; // to be specified

template <typename T>
T square(const T& v)
Expand All @@ -26,23 +26,23 @@ 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
// 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
// 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
return pc * exp(exponent); // in Pa
}

//Evaporation enthalpy of water from Nunez
double AdsorptionReaction::getEvaporationEnthalpy(double T_Ads) //in kJ/kg
// Evaporation enthalpy of water from Nunez
double AdsorptionReaction::getEvaporationEnthalpy(double T_Ads) // in kJ/kg
{
T_Ads -= 273.15;
if (T_Ads <= 10.){
Expand All @@ -64,14 +64,14 @@ double AdsorptionReaction::getEvaporationEnthalpy(double T_Ads) //in kJ/kg
}


//evaluate specific heat capacity of adsorbate follwing Nunez
// 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)
return cp; // kJ/(kg*K)
}


Expand All @@ -97,15 +97,11 @@ double AdsorptionReaction::dMolarFraction(double xm, double M_this, double M_oth
double AdsorptionReaction::getReactionRate(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 = 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 0.0; // TODO [CL] for testing only

return k_rate * (C_eq - loading); //scaled with mass fraction
return k_rate * (C_eq - loading); // scaled with mass fraction
// this the rate in terms of loading!
}

Expand All @@ -132,14 +128,10 @@ void AdsorptionReaction::getDReactionRate(const double p_Ads, const double T_Ads
}


//Evaluate adsorbtion potential A
// 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
if (A < 0.0) {
// vapour partial pressure > saturation pressure
// A = 0.0; // TODO [CL] debug output
}
double A = GAS_CONST * T_Ads * log(getEquilibriumVapourPressure(T_Ads)/p_Ads) / (M_Ads*1.e3); // in kJ/kg = J/g
return A;
}

Expand All @@ -150,7 +142,7 @@ double AdsorptionReaction::getLoading(const double rho_curr, const double rho_dr
}


//Calculate sorption entropy
// Calculate sorption entropy
double AdsorptionReaction::getEntropy(const double T_Ads, const double A) const
{
const double epsilon = 1.0e-8;
Expand All @@ -167,7 +159,6 @@ double AdsorptionReaction::getEntropy(const double T_Ads, const double A) const
return 0.0;
}

// const double dAdlnW = 2.0*epsilon/(log(characteristicCurve(A+epsilon)) - log(characteristicCurve(A-epsilon)));
return dAdlnW * getAlphaT(T_Ads);
}

Expand All @@ -178,13 +169,12 @@ double AdsorptionReaction::getEnthalpy(const double p_Ads, const double T_Ads, c
// 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 0.0; // TODO [CL] for testing only
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
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);
Expand Down
20 changes: 2 additions & 18 deletions MaterialsLib/Adsorption/Adsorption.h
View file
Expand Up @@ -37,16 +37,12 @@ class AdsorptionReaction : public Reaction

static double getLoading(const double rho_curr, const double rho_dry);

// non-virtual members
double getEquilibriumLoading(const double p_Ads, const double T_Ads, const double M_Ads)
const override;

// virtual members:
virtual ~AdsorptionReaction() = default;

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;
const double M_Ads, const double loading) const override;
/**
* @brief get_d_reaction_rate
* @param p_Ads
Expand All @@ -66,28 +62,16 @@ class AdsorptionReaction : public Reaction
virtual double dCharacteristicCurve(const double A) const = 0;

private:
// non-virtual members
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) );

// 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 dCurvePolyfrac(const double* coeffs, const double x)
Expand Down
10 changes: 4 additions & 6 deletions MaterialsLib/Adsorption/Density100MPa.cpp
View file
Expand Up @@ -34,8 +34,7 @@ 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
// 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;
Expand All @@ -44,17 +43,16 @@ double Density100MPa::getAlphaT(const double T_Ads) const
return - drhodT / rho;
}


//Characteristic curve. Return W (A)
// Characteristic curve. Return W (A)
double Density100MPa::characteristicCurve(const double A) const
{
double W = curvePolyfrac(c, A); //cm^3/g
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
return W/1.e3; // m^3/kg
}

double Density100MPa::dCharacteristicCurve(const double A) const
Expand Down
7 changes: 2 additions & 5 deletions MaterialsLib/Adsorption/DensityConst.cpp
View file
Expand Up @@ -35,15 +35,12 @@ double DensityConst::getAdsorbateDensity(const double /*T_Ads*/) const
return rhoWaterHauer(150.0+273.15);
}


//Thermal expansivity model for water found in the works of Hauer
double DensityConst::getAlphaT(const double /*T_Ads*/) const
{
return 0.0;
}


//Characteristic curve. Return W (A)
// Characteristic curve. Return W (A)
double DensityConst::characteristicCurve(const double A) const
{
double W = curvePolyfrac(c, A); //cm^3/g
Expand All @@ -52,7 +49,7 @@ double DensityConst::characteristicCurve(const double A) const
W = 0.0; // TODO [CL] debug output
}

return W/1.e3; //m^3/kg
return W/1.e3; // m^3/kg
}

double DensityConst::dCharacteristicCurve(const double A) const
Expand Down
5 changes: 1 addition & 4 deletions MaterialsLib/Adsorption/DensityCook.cpp
View file
Expand Up @@ -34,15 +34,12 @@ double DensityCook::getAdsorbateDensity(const double T_Ads) const
return rhoWaterDean(T_Ads);
}


//Thermal expansivity model for water found in the works of Hauer
double DensityCook::getAlphaT(const double T_Ads) const
{
return alphaTWaterDean(T_Ads);
}


//Characteristic curve. Return W (A)
// Characteristic curve. Return W (A)
double DensityCook::characteristicCurve(const double A) const
{
double W = curvePolyfrac(c, A); //cm^3/g
Expand Down
1 change: 0 additions & 1 deletion MaterialsLib/Adsorption/DensityCook.h
View file
Expand Up @@ -23,7 +23,6 @@ class DensityCook : public AdsorptionReaction
double dCharacteristicCurve(const double A) const;
};


inline double rhoWaterDean(const double T_Ads)
{
const double Tcel = T_Ads - 273.15;
Expand Down
37 changes: 17 additions & 20 deletions MaterialsLib/Adsorption/DensityDubinin.cpp
View file
Expand Up @@ -38,56 +38,53 @@ double DensityDubinin::getAdsorbateDensity(const double T_Ads) const
if (T_Ads < Tb) {
return rhoWaterDean(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 Tc = 647.3; // K
const double pc = 221.2e5; // Pa
// boiling point density
const double rhob = rhoWaterDean(Tb);
//state values
// 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 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
// 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
// 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 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)
// Characteristic curve. Return W (A)
double DensityDubinin::characteristicCurve(const double A) const
{
double W = curvePolyfrac(c, A); //cm^3/g
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
return W/1.e3; // m^3/kg
}

double DensityDubinin::dCharacteristicCurve(const double A) const
Expand Down
10 changes: 4 additions & 6 deletions MaterialsLib/Adsorption/DensityHauer.cpp
View file
Expand Up @@ -34,8 +34,7 @@ double DensityHauer::getAdsorbateDensity(const double T_Ads) const
return rhoWaterHauer(T_Ads);
}


//Thermal expansivity model for water found in the works of Hauer
// Thermal expansivity model for water found in the works of Hauer
double DensityHauer::getAlphaT(const double T_Ads) const
{
// data like in python script
Expand All @@ -44,17 +43,16 @@ double DensityHauer::getAlphaT(const double T_Ads) const
return alpha0/(1. - alpha0 * (T_Ads-T0)); //in 1/K
}


//Characteristic curve. Return W (A)
// Characteristic curve. Return W (A)
double DensityHauer::characteristicCurve(const double A) const
{
double W = curvePolyfrac(c, A); //cm^3/g
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
return W/1.e3; // m^3/kg
}

double DensityHauer::dCharacteristicCurve(const double A) const
Expand Down
4 changes: 2 additions & 2 deletions MaterialsLib/Adsorption/DensityHauer.h
View file
Expand Up @@ -27,9 +27,9 @@ class DensityHauer : public AdsorptionReaction
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
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
return rho0 * (1. - alpha0 * (T_Ads-T0)); // in kg/m^3
}

}

0 comments on commit bb45bee

Please sign in to comment.