Commit 86c20a25 authored by Hang Zhou's avatar Hang Zhou

1. Add CCK model for char reactions; 2. Changing coal type...

1. Add CCK model for char reactions; 2. Changing coal type "Given_From_InputFile" to be "Custom_by_User".
parent 2b7e9342
Pipeline #1083 failed with stage
in 18 minutes and 9 seconds
......@@ -1143,6 +1143,68 @@ namespace ReactorEnsembleUtil{
}
};
/**
* @class AverageExpression
* @ingroup ReactorEnsembleUtil
* @brief Calculate the average of two expressions
* @author Hang
*
* The result is given as
* \f[
* Average = (Expression1 + Expression2) / 2
* \f]
*/
template<typename FieldT>
class AverageExpression : public Expr::Expression<FieldT> {
DECLARE_FIELDS( FieldT, expr1_, expr2_ );
AverageExpression( const Expr::Tag tag1,
const Expr::Tag tag2 )
: Expr::Expression<FieldT>(){
this->set_gpu_runnable( true );
expr1_ = this->template create_field_request<FieldT>( tag1 );
expr2_ = this->template create_field_request<FieldT>( tag2 );
}
public:
class Builder : public Expr::ExpressionBuilder {
const Expr::Tag tag1_, tag2_;
public:
/**
* @brief The mechanism for building a AverageExpression object
* @tparam FieldT
* @param averageTag average result of two expressions
* @param tag1 The first expression
* @param tag2 The second expression
*/
Builder( const Expr::Tag averageTag,
const Expr::Tag tag1,
const Expr::Tag tag2 ) :
Expr::ExpressionBuilder( averageTag ),
tag1_( tag1 ),
tag2_( tag2 ){}
~Builder(){}
Expr::ExpressionBase* build() const{
return new AverageExpression <FieldT>( tag1_, tag2_ );
}
};
void evaluate(){
using namespace SpatialOps;
this->value() <<= (expr1_->field_ref() + expr2_->field_ref()) / 2.0;
}
void sensitivity( const Expr::Tag& sensVarTag ){
FieldT& dfdv = this->sensitivity_result( sensVarTag );
dfdv <<= (expr1_->sens_field_ref( sensVarTag ) + expr2_->sens_field_ref( sensVarTag )) / 2.0;
}
};
}
......
......@@ -118,11 +118,13 @@ namespace Particles {
massEachPartTags = Expr::tag_list( particleNumArray, state, "mass_each_part_" );
partNumPerVolumeTags = Expr::tag_list( particleNumArray, state, "num_part_per_volume_" );
partNumFracTags = Expr::tag_list( particleNumArray, state, "num_frac_part_" );
partNumFracInitialTags = Expr::tag_list( particleNumArray, state, "num_frac_part_", "_initial" );
partTempTags = Expr::tag_list( particleNumArray, state, "temp_part_" );
partCpTags = Expr::tag_list( particleNumArray, state, "cp_part_" );
partMassFracTags = Expr::tag_list( particleNumArray, state, "mass_frac_part_" );
partRhoTags = Expr::tag_list( particleNumArray, state, "rho_part_" );
partSizeTags = Expr::tag_list( particleNumArray, state, "size_part_" );
partSizeInitialTags = Expr::tag_list( particleNumArray, state, "size_part_", "_initial" );
partConvecTags = Expr::tag_list( particleNumArray, state, "convection_part_" );
partGasRadiaTags = Expr::tag_list( particleNumArray, state, "radiation_part_gas_" );
partWallRadiaTags = Expr::tag_list( particleNumArray, state, "radiation_part_wall_" );
......@@ -271,6 +273,7 @@ namespace Particles {
particle_size(rootParser_["Particles"], ParticleSizeArray, nparsize_);
for( size_t i = 0;i < nparsize_;i++){
initRoots_.insert( initFactory_.register_expression( new ConstantT( tagsPart.partSizeTags[i], ParticleSizeArray[i] ) ) );
initRoots_.insert( initFactory_.register_expression( new ConstantT( tagsPart.partSizeInitialTags[i], ParticleSizeArray[i] ) ) );
}
for( size_t i = 0;i < nparsize_;i++){
......@@ -324,6 +327,7 @@ namespace Particles {
// viscosity is set to be constant in present code.
// execFactory.register_expression( new ViscosityT( tagsE.dynviscosityTag, tagsE. tempTag, tagsE.massTags ) );
execFactory.register_expression( new PlaceHolderT( tagsE.dynviscosityTag ));
integrator_->register_root_expression( new PlaceHolderT( tagsEPart.partRhoInitTag ));
if(reactorType_ == "PSRConstantVolume"){
execFactory.register_expression( new HeatCpT( tagsE.cpTag, tagsE.tempTag, tagsE.massTags ));
}
......@@ -332,6 +336,7 @@ namespace Particles {
execFactory.register_expression( new ConstantT( tagsE.scNumberTag, 1.0 ) ); // set sc to be constant
for( size_t i = 0;i < nparsize_;i++){
integrator_->register_root_expression( new PlaceHolderT( tagsEPart.partSizeInitialTags[i] ));
if(rootParser_["Particles"]["ParticleType"].as<std::string>() != "Solid"){
execFactory.register_expression( new ShNumberT( tagsEPart.shNumberTags[i], tagsE.scNumberTag, tagsEPart.reynoldsNumberTags[i] ) );
}
......@@ -472,8 +477,10 @@ namespace Particles {
else {
integrator_->copy_from_initial_condition_to_execution<FieldT>( tagsEPart.partSizeTags[i].name());
}
integrator_->copy_from_initial_condition_to_execution<FieldT>( tagsE.dynviscosityTag.name());
integrator_->copy_from_initial_condition_to_execution<FieldT>( tagsEPart.partSizeInitialTags[i].name());
}
integrator_->copy_from_initial_condition_to_execution<FieldT>( tagsE.dynviscosityTag.name());
integrator_->copy_from_initial_condition_to_execution<FieldT>( tagsEPart.partRhoInitTag.name());
}
for( size_t i = 0;i < nparsize_;i++ ){
if(reactorType_!="PFR"){
......@@ -491,8 +498,10 @@ namespace Particles {
else {
integrator_->lock_field<FieldT>( tagsEPart.partSizeTags[i] );
}
integrator_->lock_field<FieldT>( tagsEPart.partSizeInitialTags[i] );
}
integrator_->lock_field<FieldT>( tagsE.dynviscosityTag );
integrator_->lock_field<FieldT>( tagsEPart.partRhoInitTag );
if(coalimplement_){
CoalInterface_->initialize_coal();
......@@ -831,8 +840,10 @@ namespace Particles {
restartVars.push_back(tagsEPart.partSizeTags[i].name());
}
restartVars.push_back(tagsEPart.partMassFracInitialTags[i].name());
restartVars.push_back(tagsEPart.partSizeInitialTags[i].name());
}
restartVars.push_back(tagsE.dynviscosityTag.name());
restartVars.push_back(tagsEPart.partRhoInitTag.name());
if(coalimplement_){
CoalInterface_->restart_var_coal(restartVars);
}
......
......@@ -68,11 +68,13 @@ namespace Particles {
Expr::TagList massEachPartTags, ///< mass of each particle
partNumPerVolumeTags, ///< number of particles per volume
partNumFracTags, ///< number of particles per mass of gas in the reactor
partNumFracInitialTags, ///< initial number of particles per mass of gas in the reactor
partTempTags, ///< particle temperature
partCpTags, ///< particle heat capacity
partMassFracTags, ///< particle mass fraction: (kg particle)/(kg gas)
partRhoTags, ///< particle density
partSizeTags, ///< particle diameter
partSizeInitialTags, ///< initial particle diameter
partConvecTags, ///< convection term in particle temperature equation
partGasRadiaTags, ///< radiation term between particle and gas-phase in particle temperature equation
partWallRadiaTags, ///< radiation term between particle and wall in particle temperature equation
......
/*
* The MIT License
*
* Copyright (c) 2015-2017 The University of Utah
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal in the Software without restriction, including without limitation the
* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
* sell copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
* IN THE SOFTWARE.
*/
#ifndef AshFilm_Expr_h
#define AshFilm_Expr_h
#include <expression/Expression.h>
#include "CCKData.h"
namespace CCK{
/**
* \class AshFilm
*
* \brief Calculates ash film thickness and porosity and core diameter. These
* are calculated together because the values are interdependent.
*/
template< typename FieldT >
class AshFilm
: public Expr::Expression<FieldT>
{
DECLARE_FIELDS( FieldT, prtDiameter_, prtDensity_d_, initPrtDensity_d_, ashMassFrac_d_, initAshMassFrac_d_ )
const CCKData& cckData_;
const double thetaMin_, deltaMin_, rhoA_np_;
/*
* \param prtDiameter : particle diameter
* \param prtDensity_d : (mass of ash + char)/(particle volume)
* \param initPrtDensity_d: initial value of prtDensity_d
* \param ashMassFrac_d : (ash mass)/(ash mass + char mass)
*/
AshFilm( const Expr::Tag& prtDiameterTag,
const Expr::Tag& prtDensity_dTag,
const Expr::Tag& initPrtDensity_dTag,
const Expr::Tag& ashMassFrac_dTag,
const Expr::Tag& initAshMassFrac_dTag,
const CCKData& cckData )
: Expr::Expression<FieldT>(),
cckData_ ( cckData ),
thetaMin_(cckData_.get_min_ash_porosity()),
deltaMin_(cckData_.get_min_ash_thickness()),
rhoA_np_ (cckData_.get_nonporous_ash_density())
{
prtDiameter_ = this->template create_field_request<FieldT>( prtDiameterTag );
prtDensity_d_ = this->template create_field_request<FieldT>( prtDensity_dTag );
initPrtDensity_d_ = this->template create_field_request<FieldT>( initPrtDensity_dTag );
ashMassFrac_d_ = this->template create_field_request<FieldT>( ashMassFrac_dTag );
initAshMassFrac_d_= this->template create_field_request<FieldT>( initAshMassFrac_dTag);
}
public:
class Builder : public Expr::ExpressionBuilder
{
const Expr::Tag prtDiameterTag_, prtDensity_dTag_, initPrtDensity_dTag_, ashMassFrac_dTag_, initAshMassFrac_dTag_;
const CCKData& cckData_;
public:
/**
* @brief Build a AshFilm expression
* @param resultTags the tags for the value that this expression computes
*/
Builder( const Expr::TagList& resultTags,
const Expr::Tag& prtDiameterTag,
const Expr::Tag& prtDensity_dTag,
const Expr::Tag& initPrtDensity_dTag,
const Expr::Tag& ashMassFrac_dTag,
const Expr::Tag& initAshMassFrac_dTag,
const CCKData& cckData,
const int nghost = DEFAULT_NUMBER_OF_GHOSTS )
: ExpressionBuilder( resultTags, nghost ),
prtDiameterTag_ ( prtDiameterTag ),
prtDensity_dTag_ ( prtDensity_dTag ),
initPrtDensity_dTag_ ( initPrtDensity_dTag ),
ashMassFrac_dTag_ ( ashMassFrac_dTag ),
initAshMassFrac_dTag_( initAshMassFrac_dTag),
cckData_ ( cckData )
{}
~Builder(){}
Expr::ExpressionBase* build() const{
return new AshFilm<FieldT>( prtDiameterTag_, prtDensity_dTag_, initPrtDensity_dTag_, ashMassFrac_dTag_, initAshMassFrac_dTag_, cckData_ );
}
};
void evaluate()
{
using namespace SpatialOps;
typedef typename FieldT::const_iterator FieldIter;
typename Expr::Expression<FieldT>::ValVec& result = this->get_value_vec();
typename FieldT::iterator itheta = result[0]->begin(); // ash film porosity
typename FieldT::iterator idelta = result[1]->begin(); // ash film thickness
typename FieldT::iterator idC = result[2]->begin(); // carbonacious core diameter
FieldIter idP = prtDiameter_ ->field_ref().begin();
FieldIter irhoP_d = prtDensity_d_ ->field_ref().begin();
FieldIter irhoP0_d = initPrtDensity_d_ ->field_ref().begin();
FieldIter ixA_d = ashMassFrac_d_ ->field_ref().begin();
FieldIter ixA0_d = initAshMassFrac_d_ ->field_ref().begin();
const typename FieldT::const_iterator iEnd = prtDiameter_->field_ref().end();
for(; idP != iEnd; ++idP, ++irhoP_d, ++irhoP0_d, ++ixA_d, ++ixA0_d,
++itheta, ++idelta, ++idC ){
// Calculate trial value of particle core diameter and ash film thickness.
const double dC_t = pow(
(*ixA_d * *irhoP_d - rhoA_np_*( 1.0 - thetaMin_ ) )
/ ( *ixA0_d * *irhoP0_d - rhoA_np_*( 1.0 - thetaMin_ ) )
, 1.0/3.0 )* *idP;
const double delta_t = 0.5*(*idP - dC_t);
/* Calculate ash film porosity. If the ash film thickness is larger than its minimum
* value, the film porosity will be set to its minimum value and the core diameter and
* film thickness will be set to their trial values...
*/
if( delta_t >= deltaMin_ ){
*itheta = thetaMin_;
*idelta = delta_t;
*idC = dC_t;
}
/* ...otherwise, film thickness is set to its minimum value, film porosity is calculated
* using the minimum film thickness, and core diameter is calculated with the new value
* for film porosity.
*/
else{
*itheta = deltaMin_ > 0?
1.0 - ( *ixA_d * *irhoP_d - *ixA0_d * *irhoP0_d*pow(1 - 2*deltaMin_/(*idP), 3.0) )
/ ( rhoA_np_*( 1 - pow(1 - 2*deltaMin_/(*idP), 3.0) ) ):
thetaMin_;
*idelta = deltaMin_;
*idC = *idP - 2*deltaMin_;
}
}//for
}
void sensitivity( const Expr::Tag& var){
using namespace SpatialOps;
const Expr::TagList& targetTags = this->get_tags();
typename Expr::Expression<FieldT>::ValVec& result = this->get_value_vec();
std::vector<FieldT> dvaluedv;
for( size_t i = 0;i < targetTags.size();++i ){
dvaluedv.push_back( this->sensitivity_result( targetTags[i], var ));
}
const FieldT& dP = prtDiameter_ ->field_ref(); const FieldT& ddPdv = prtDiameter_ ->sens_field_ref(var);
const FieldT& rhoP_d = prtDensity_d_ ->field_ref(); const FieldT& drhoP_ddv = prtDensity_d_ ->sens_field_ref(var);
const FieldT& rhoP0_d = initPrtDensity_d_ ->field_ref(); const FieldT& drhoP0_ddv = initPrtDensity_d_ ->sens_field_ref(var);
const FieldT& xA_d = ashMassFrac_d_ ->field_ref(); const FieldT& dxA_ddv = ashMassFrac_d_ ->sens_field_ref(var);
const FieldT& xA0_d = initAshMassFrac_d_ ->field_ref(); const FieldT& dxA0_ddv = initAshMassFrac_d_ ->sens_field_ref(var);
SpatFldPtr<FieldT> dC_t = SpatialFieldStore::get<FieldT,FieldT>( dP );
SpatFldPtr<FieldT> delta_t = SpatialFieldStore::get<FieldT,FieldT>( dP );
SpatFldPtr<FieldT> ddC_tdv = SpatialFieldStore::get<FieldT,FieldT>( dP );
SpatFldPtr<FieldT> ddelta_tdv = SpatialFieldStore::get<FieldT,FieldT>( dP );
*dC_t <<= pow((xA_d * rhoP_d - rhoA_np_*( 1.0 - thetaMin_ ))/(xA0_d * rhoP0_d - rhoA_np_*( 1.0 - thetaMin_ ) ), 1.0/3.0 )* dP;
*delta_t <<= 0.5*(dP - *dC_t);
*ddC_tdv <<= 1.0/3.0 * pow((xA_d * rhoP_d - rhoA_np_*( 1.0 - thetaMin_ ))/(xA0_d * rhoP0_d - rhoA_np_*( 1.0 - thetaMin_ ) ), -2.0/3.0 )* dP
* ( (dxA_ddv*rhoP_d+xA_d*drhoP_ddv)/(xA0_d * rhoP0_d - rhoA_np_*( 1.0 - thetaMin_ ))
- (xA_d*rhoP_d-rhoA_np_*(1.0-thetaMin_))/square(xA0_d*rhoP0_d-rhoA_np_*(1.0-thetaMin_)) * (dxA0_ddv*rhoP0_d+xA0_d*drhoP0_ddv)
)
+ pow((xA_d * rhoP_d - rhoA_np_*( 1.0 - thetaMin_ ))/(xA0_d * rhoP0_d - rhoA_np_*( 1.0 - thetaMin_ ) ), 1.0/3.0 )* ddPdv;
*ddelta_tdv <<= 0.5*(ddPdv - *ddC_tdv);
dvaluedv[0] <<= cond(*delta_t>=deltaMin_ or deltaMin_<=0, 0)
(- (dxA_ddv*rhoP_d+xA_d*drhoP_ddv-(dxA0_ddv*rhoP0_d+xA0_d*drhoP0_ddv)*pow(1-2*deltaMin_/dP, 3.0)
- xA0_d*rhoP0_d*3.0*square(1-2*deltaMin_/dP)*2*deltaMin_/square(dP)*ddPdv)
/ (rhoA_np_*(1-pow(1-2*deltaMin_/dP, 3.0)))
- (xA_d*rhoP_d-xA0_d*rhoP0_d*pow(1-2*deltaMin_/dP, 3.0))/rhoA_np_/square((1-pow(1-2*deltaMin_/dP, 3.0)))
* 3.0*square(1-2*deltaMin_/dP) * 2*deltaMin_/square(dP)*ddPdv);
dvaluedv[1] <<= cond(*delta_t>=deltaMin_, *ddelta_tdv)(0);
dvaluedv[2] <<= cond(*delta_t>=deltaMin_, *ddC_tdv)(ddPdv);
}
};
}// namespace CCK
#endif // AshFilm_Expr_h
This diff is collapsed.
#include "CCKData.h"
#include <stdexcept>
#include <sstream>
#include <cmath>
namespace CCK{
CCKData::CCKData( const Coal::CoalType sel,
const YAML::Node& particleParser)
: coalType_(sel),
coalComp_(sel, particleParser),
charData_(sel, particleParser)
// note: neD and its standard deviation have units ln(kCal/mol)
{
/* Pre-exponential factors for each of the char surface reactions Evolve in time
* according to the following formula:
*
* A_i = A_i0*sqrt(f),
*
* where f is the fraction of active sites remaining in the char. The model implented
* assumes the initial value of f is a lognormal distribution wrt eD. f evolves in
* time according to the following ODE:
*
* df/dt = -f*aD*exp(-eD/RT).
*/
//Assemble vector of eD values:
size_t numel = 30;
eDVec_.resize( numel );
eDVec_[0] = 1;
eDVec_[numel-1] = 121;
double step = (eDVec_[numel-1] - eDVec_[0])/(numel - 1);
for(size_t i = 1; i<numel-1; ++i){
eDVec_[i] = eDVec_[0] + i*step;
}
/* The mode-of-burning parameter should be set to 0.2 for oxidation
* and 0.95 for gasification (without oxidation)
*/
modeOfBurningParam_ = 0.2; // for oxidation
tau_f_ = 12; // tau/f see section 6 of [1] for more details
ashDensity_ = 2650; // kg/m^3
c_ = coalComp_.get_C();
o_ = coalComp_.get_O();
fixedCarbon_ = coalComp_.get_fixed_c();
volatiles_ = coalComp_.get_vm();
xA_ = coalComp_.get_ash();
percentC_ = c_*100; // correlations use % basis
aFwd0_.assign(8,0.0);
aRev0_.assign(8,0.0);
eaFwd_.assign(8,0.0);
eaRev_.assign(8,0.0);
mean_neD_ = 2.8; // log-mean of thermal annealing activation energy
aD_ = exp(18.3); // thermal annealing frequency factor
stdDev_neD_ = 0.46; // std deviation of neD
deltaMin_ = 0.0; // minimum ash film thickness
thetaMin_ = 0.21; // minimum ash porosity
epsilon_ = 0.47; // char porosity
switch(coalType_){
case Coal::Highvale:
case Coal::Highvale_Char:
mean_neD_ = 2.0;
break;
case Coal::Eastern_Bituminous:
case Coal::Eastern_Bituminous_Char:
case Coal::Utah_Skyline:
case Coal::Utah_Skyline_Char_10atm:
case Coal::Utah_Skyline_Char_12_5atm:
case Coal::Utah_Skyline_Char_15atm:
case Coal::NorthDakota_Lignite:
case Coal::Gillette_Subbituminous:
case Coal::MontanaRosebud_Subbituminous:
case Coal::Illinois_Bituminous:
case Coal::Kentucky_Bituminous:
case Coal::Pittsburgh_Shaddix:
case Coal::Pittsburgh_Bituminous:
case Coal::Black_Thunder:
case Coal::Shenmu:
case Coal::Guizhou:
case Coal::Russian_Bituminous:
case Coal::Utah_Bituminous:
case Coal::Illinois_No_6:
case Coal::Illinois_No_6_Char:
case Coal::Custom_by_User:
break;
default:
std::ostringstream msg;
msg << __FILE__ << " : " << __LINE__ << std::endl
<< "Unsupported coal type" << std::endl
<< std::endl;
throw std::runtime_error( msg.str() );
}
switch(coalType_){
case Coal::Utah_Skyline:
case Coal::Utah_Skyline_Char_10atm:
case Coal::Utah_Skyline_Char_12_5atm:
case Coal::Utah_Skyline_Char_15atm:
// These data are for char generated from Utah Skyline coal.
// See Tables 3 & 5 from [1] for more details.
aFwd0_[6] = 9.97e+8; // s^-1
eaFwd_[6] = 166.0e+3; // J/mol
break;
case Coal::NorthDakota_Lignite:
case Coal::Gillette_Subbituminous:
case Coal::MontanaRosebud_Subbituminous:
case Coal::Illinois_Bituminous:
case Coal::Kentucky_Bituminous:
case Coal::Pittsburgh_Shaddix:
case Coal::Pittsburgh_Bituminous:
case Coal::Black_Thunder:
case Coal::Shenmu:
case Coal::Guizhou:
case Coal::Russian_Bituminous:
case Coal::Utah_Bituminous:
case Coal::Highvale:
case Coal::Highvale_Char:
case Coal::Eastern_Bituminous:
case Coal::Eastern_Bituminous_Char:
case Coal::Illinois_No_6:
case Coal::Illinois_No_6_Char:
case Coal::Custom_by_User:
aFwd0_[6] = exp(-0.3672*c_/o_ + 21.33);
eaFwd_[6] = 146.4e3;
break;
default:
std::ostringstream msg;
msg << __FILE__ << " : " << __LINE__ << std::endl
<< "Unsupported coal type" << std::endl
<< std::endl;
throw std::runtime_error( msg.str() );
}
/* The kinetic parameters below are from [3]. Afwd0[2] is dependent on
* pressure, so it is not set here.
*/
a0_over_a2_ = 1.0; // aFwd0_[0]/aFwd0_[2] (m^3/mol)
a1_over_a2_ = 0.05; // aFwd0_[1]/aFwd0_[2] (m^3/mol)
eaFwd_[0] = 25.0e+3;
eaFwd_[1] = 117.0e+3;
eaFwd_[2] = 133.8e+3;
/* Correlations for kinetic parameters below are given in [4], eqs. 6.33-6.45.
* 7 forward rate constants and 2 reverse rate constants (out of 8) are
* used (some defined above).
*/
aFwd0_[3] = (1.84e+3)*exp(-0.073*percentC_)*aFwd0_[6]/101325.0; // (6.33) (Pa^-1)
// aFwd0_[4] not set
aFwd0_[5] = (percentC_ <= 90.6)? 0.5*aFwd0_[6]/101325.0: // (6.35) (Pa^-1)
(0.021*percentC_ - 1.86)*aFwd0_[6]/101325.0; // (6.36) (Pa^-1)
aFwd0_[7] = 7.9e-5*aFwd0_[6]/101325.0; // (6.38) (Pa^-1)
aRev0_[3] = ( (3.57e-5)*percentC_ - 1.73e-3 )*aFwd0_[6]/101325.0;// (6.34) (Pa^-1)
aRev0_[5] = ( -(3.68e-8)*percentC_ + 3.2e-6 )*aFwd0_[6]/101325.0;// (6.37) (Pa^-1)
// activation energies in J/mol
eaFwd_[3] = eaFwd_[6] + 54.0e3; // (6.41)
// eaFwd_[4] not set
eaFwd_[5] = eaFwd_[6] - 16.0e+3; // (6.43)
eaFwd_[7] = eaFwd_[6] -26.0e+3; // (6.45)
eaRev_[3] = eaFwd_[6] - 53.9e+3; // (6.42)
eaRev_[5] = eaFwd_[6] - 156.0e+3; // (6.44)
}
//---------------------------------------------------------------------------
std::vector<double>
CCKData::forward_rate_constants(const double& pressure,
const double& pTemp,
const double& saFactor ) const{
/* calculation of aFwd0_[3] based on correlation given by eq. (10) on
* page 468 [3]. Correlations for 100kPa and "elevated" pressure are
* provided. Here The limit for elevated pressure will be defined to be
* 200kPa. For pressures between these two values aFwd0_[3] will be
* calculated with an exponential interpolation of the two separate
* correlations
*
* saFactor is the product of coefficients accounting for deactivation
* due to thermal annealing and conversion of char.
*/
std::vector<double> kFwd; // vector of forward rate constants
kFwd.assign(8,0.0);
const double eLow = 14.38 - 0.0764*percentC_; //"0.1 MPa only"
const double eHi = 12.22 - 0.0535*percentC_; //"elevated pressures"
const double pLow = 1.0e+5; //
const double pHi = 2.0e+5;
double expon;
if( pressure < pHi || pressure > pLow){
expon = (eHi - eLow)*(pressure - pLow )/(pHi - pLow) + eLow;
}
else if( pressure > pHi ) expon = eHi;
else expon = eLow;
const double aFwd0_2 = pow(expon,10.0);
const double RT = 8.3144621 * pTemp;
kFwd[0] = saFactor*(aFwd0_2*a0_over_a2_)*exp( -eaFwd_[0]/RT ); //1
kFwd[1] = saFactor*(aFwd0_2*a1_over_a2_)*exp( -eaFwd_[1]/RT ); //2
kFwd[2] = saFactor*aFwd0_2*exp( -eaFwd_[2]/RT ); //3
kFwd[3] = saFactor*aFwd0_[3]*exp( -eaFwd_[3]/RT ); //4
// kFwd[4] (rxn 5) is included in the gamma parameter
kFwd[5] = saFactor*aFwd0_[5]*exp( -eaFwd_[5]/RT ); //6
kFwd[6] = saFactor*aFwd0_[6]*exp( -eaFwd_[6]/RT ); //7
kFwd[7] = saFactor*aFwd0_[7]*exp( -eaFwd_[7]/RT ); //8
return kFwd;
}
// -----------------------------------------------------------------------------
std::vector<double>
CCKData::reverse_rate_constants( const double& pTemp,
const double& saFactor) const{
/* saFactor is the product of coefficients accounting for deactivation
* due to thermal annealing and conversion of char.
*/
std::vector<double> kRev; // vector reverse rate constants
kRev.assign(8,0.0);
const double RT = 8.3144621 * pTemp;
kRev[3] = saFactor*aRev0_[3]*exp( -eaRev_[3]/RT ); //4
kRev[5] = saFactor*aRev0_[5]*exp( -eaRev_[5]/RT ); //6
return kRev;
}
// -----------------------------------------------------------------------------
const double calc_gamma( const double& pTemp){
// calculates gamma from eq. 43 in [4]
double T;
if(pTemp < 1573 || pTemp > 1073 ) T = pTemp;
else if(pTemp >= 1573 ) T = 1573;
else T = 1073;
return (6.92e+2)*exp( (-5.0e+4)/(8.31446211*T) );
}
// -----------------------------------------------------------------------------
}// namespace CCK
/*
[1] A. Lewis et. al. Pulverized Steam Gasification Rates of Three Bituminous
Coal Chars in an Entrained-Flow Reactor at Pressurized Conditions. Energy
and Fuels. 2015, 29, 1479-1493. http://pubs.acs.org/doi/abs/10.1021/ef502608y
[2] R. Shurtz. Effects of Pressure on the Properties of Coal Char Under
Gasification Conditions at High Initial Heating Rates. (2011). All Theses
and Dissertations. Paper 2877. http://scholarsarchive.byu.edu/etd/2877/
[3] S. Niksa et. al. Coal conversion submodels for design applications at elevated