Revision d8b826d344caa66abbcc3d86852ea99e6a14a329 authored by Jörg Buchwald on 08 March 2021, 20:38:20 UTC, committed by Dmitry Yu. Naumov on 14 March 2021, 16:08:33 UTC
1 parent 883db3d
ReactionCaOH2.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, 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 "ReactionCaOH2.h"
#include <cassert>
#include "Adsorption.h"
#include "BaseLib/Error.h"
#include "MaterialLib/PhysicalConstant.h"
namespace Adsorption
{
const double ReactionCaOH2::_reaction_enthalpy = -1.12e+05;
const double ReactionCaOH2::_reaction_entropy = -143.5;
const double ReactionCaOH2::_M_carrier =
MaterialLib::PhysicalConstant::MolarMass::N2;
const double ReactionCaOH2::_M_react =
MaterialLib::PhysicalConstant::MolarMass::Water;
const double ReactionCaOH2::_tol_l = 1e-4;
const double ReactionCaOH2::_tol_u = 1.0 - 1e-4;
const double ReactionCaOH2::_tol_rho = 0.1;
const double ReactionCaOH2::rho_low = 1656.0;
const double ReactionCaOH2::rho_up = 2200.0;
double ReactionCaOH2::getEnthalpy(const double /*p_Ads*/,
const double /*T_Ads*/,
const double /*M_Ads*/) const
{
return - _reaction_enthalpy/_M_react;
}
double ReactionCaOH2::getReactionRate(const double /*p_Ads*/,
const double /*T_Ads*/,
const double /*M_Ads*/,
const double /*loading*/) const
{
OGS_FATAL("get_reaction_rate do not call directly");
}
double ReactionCaOH2::getReactionRate(double const solid_density)
{
_rho_s = solid_density;
calculateQR();
return _qR;
}
void ReactionCaOH2::updateParam(
double T_solid,
double p_gas,
double x_react,
double rho_s_initial)
{
_T_s = T_solid;
_p_gas = p_gas / 1e5; // convert Pa to bar
_x_react = x_react;
_rho_s = rho_s_initial;
}
void ReactionCaOH2::calculateQR()
{
// Convert mass fraction into mole fraction
const double mol_frac_react = AdsorptionReaction::getMolarFraction(_x_react, _M_react, _M_carrier);
_p_r_g = std::max(mol_frac_react * _p_gas, 1.0e-3); // avoid illdefined log
setChemicalEquilibrium();
const double dXdt = CaHydration();
_qR = (rho_up - rho_low) * dXdt;
}
// determine equilibrium temperature and pressure according to van't Hoff
void ReactionCaOH2::setChemicalEquilibrium()
{
const double R = MaterialLib::PhysicalConstant::IdealGasConstant;
_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 = std::exp((_reaction_enthalpy/R)/_T_s - (_reaction_entropy/R));
}
double ReactionCaOH2::CaHydration()
{
const double R = MaterialLib::PhysicalConstant::IdealGasConstant;
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;
}
} // namespace Adsorption
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...