General commit for a reworking of the Species reference state thermo
calculation. This is a reclarification of the reference state thermo calculations for individual species, and an expansion to handle liquid phase thermo needs. There is now a virtual base class for the calculation of reference state thermo functions for individual species. It is called SpeciesThermoInterpType. There is also a class which allows for a complete general calculation of the reference state species thermo for a phase, GeneralSpeciesThermo. Some of this new functionality may be relegated to ifdef blocks in the future to limit the amount of code for users who don't need the new functionality.
This commit is contained in:
Harry Moffat
parent
ece819061c
commit
5de3c69245
18 changed files with 1797 additions and 82 deletions
132
Cantera/src/ConstCpPoly.cpp
Normal file
132
Cantera/src/ConstCpPoly.cpp
Normal file
|
|
@ -0,0 +1,132 @@
|
|||
/**
|
||||
* @file ConstCpPoly.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
|
||||
|
||||
#include "ConstCpPoly.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
ConstCpPoly::ConstCpPoly()
|
||||
: m_t0(0.0),
|
||||
m_cp0_R(0.0),
|
||||
m_h0_R(0.0),
|
||||
m_s0_R(0.0),
|
||||
m_logt0(0.0),
|
||||
m_lowT(0.0),
|
||||
m_highT(0.0),
|
||||
m_Pref(0.0),
|
||||
m_index(0) {
|
||||
}
|
||||
|
||||
ConstCpPoly::ConstCpPoly(int n, doublereal tlow, doublereal thigh,
|
||||
doublereal pref,
|
||||
const doublereal* coeffs) :
|
||||
m_lowT (tlow),
|
||||
m_highT (thigh),
|
||||
m_Pref (pref),
|
||||
m_index (n) {
|
||||
m_t0 = coeffs[0];
|
||||
m_h0_R = coeffs[1] / GasConstant;
|
||||
m_s0_R = coeffs[2] / GasConstant;
|
||||
m_cp0_R = coeffs[3] / GasConstant;
|
||||
m_logt0 = log(m_t0);
|
||||
}
|
||||
|
||||
ConstCpPoly::ConstCpPoly(const ConstCpPoly& b) :
|
||||
m_t0 (b.m_t0),
|
||||
m_cp0_R (b.m_cp0_R),
|
||||
m_h0_R (b.m_h0_R),
|
||||
m_s0_R (b.m_s0_R),
|
||||
m_logt0 (b.m_logt0),
|
||||
m_lowT (b.m_lowT),
|
||||
m_highT (b.m_highT),
|
||||
m_Pref (b.m_Pref),
|
||||
m_index (b.m_index)
|
||||
{
|
||||
}
|
||||
|
||||
ConstCpPoly& ConstCpPoly::operator=(const ConstCpPoly& b) {
|
||||
if (&b != this) {
|
||||
m_t0 = b.m_t0;
|
||||
m_cp0_R = b.m_cp0_R;
|
||||
m_h0_R = b.m_h0_R;
|
||||
m_s0_R = b.m_s0_R;
|
||||
m_logt0 = b.m_logt0;
|
||||
m_lowT = b.m_lowT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
ConstCpPoly::~ConstCpPoly(){}
|
||||
|
||||
SpeciesThermoInterpType *
|
||||
ConstCpPoly::duplMyselfAsSpeciesThermoInterpType() const {
|
||||
ConstCpPoly* newCCP = new ConstCpPoly(*this);
|
||||
return (SpeciesThermoInterpType*) newCCP;
|
||||
}
|
||||
|
||||
doublereal ConstCpPoly::minTemp() const {
|
||||
return m_lowT;
|
||||
}
|
||||
doublereal ConstCpPoly::maxTemp() const {
|
||||
return m_highT;
|
||||
}
|
||||
doublereal ConstCpPoly::refPressure() const {
|
||||
return m_Pref;
|
||||
}
|
||||
|
||||
void ConstCpPoly::updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
double t = *tt;
|
||||
doublereal logt = log(t);
|
||||
doublereal rt = 1.0/t;
|
||||
cp_R[m_index] = m_cp0_R;
|
||||
h_RT[m_index] = rt*(m_h0_R + (t - m_t0) * m_cp0_R);
|
||||
s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0);
|
||||
}
|
||||
|
||||
void ConstCpPoly::updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
doublereal logt = log(temp);
|
||||
doublereal rt = 1.0/temp;
|
||||
cp_R[m_index] = m_cp0_R;
|
||||
h_RT[m_index] = rt*(m_h0_R + (temp - m_t0) * m_cp0_R);
|
||||
s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0);
|
||||
}
|
||||
|
||||
void ConstCpPoly::reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = CONSTANT_CP;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
coeffs[0] = m_t0;
|
||||
coeffs[1] = m_h0_R * GasConstant;
|
||||
coeffs[2] = m_s0_R * GasConstant;
|
||||
coeffs[3] = m_cp0_R * GasConstant;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
66
Cantera/src/ConstCpPoly.h
Normal file
66
Cantera/src/ConstCpPoly.h
Normal file
|
|
@ -0,0 +1,66 @@
|
|||
/**
|
||||
* @file ConstCpPoly.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
|
||||
#ifndef CT_CONSTCPPOLY_H
|
||||
#define CT_CONSTCPPOLY_H
|
||||
|
||||
#include "SpeciesThermoInterpType.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
|
||||
class ConstCpPoly: public SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
ConstCpPoly();
|
||||
ConstCpPoly(int n, doublereal tlow, doublereal thigh,
|
||||
doublereal pref,
|
||||
const doublereal* coeffs);
|
||||
ConstCpPoly(const ConstCpPoly&);
|
||||
ConstCpPoly& operator=(const ConstCpPoly&);
|
||||
virtual ~ConstCpPoly();
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const;
|
||||
|
||||
doublereal minTemp() const;
|
||||
doublereal maxTemp() const;
|
||||
doublereal refPressure() const;
|
||||
virtual int reportType() const { return CONSTANT_CP; }
|
||||
|
||||
void updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const;
|
||||
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const;
|
||||
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const;
|
||||
protected:
|
||||
doublereal m_t0;
|
||||
doublereal m_cp0_R;
|
||||
doublereal m_h0_R;
|
||||
doublereal m_s0_R;
|
||||
doublereal m_logt0;
|
||||
doublereal m_lowT, m_highT, m_Pref;
|
||||
int m_index;
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
228
Cantera/src/GeneralSpeciesThermo.cpp
Normal file
228
Cantera/src/GeneralSpeciesThermo.cpp
Normal file
|
|
@ -0,0 +1,228 @@
|
|||
/**
|
||||
* @file GeneralSpeciesThermo.cpp
|
||||
*
|
||||
*/
|
||||
|
||||
// Copyright 2001-2004 California Institute of Technology
|
||||
|
||||
#include "GeneralSpeciesThermo.h"
|
||||
#include "NasaPoly1.h"
|
||||
#include "NasaPoly2.h"
|
||||
#include "ShomatePoly.h"
|
||||
#include "ConstCpPoly.h"
|
||||
#include "Mu0Poly.h"
|
||||
#include "SpeciesThermoFactory.h"
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
|
||||
/*
|
||||
* Constructors
|
||||
*/
|
||||
GeneralSpeciesThermo::GeneralSpeciesThermo() :
|
||||
SpeciesThermo(),
|
||||
m_tlow_max(0.0),
|
||||
m_thigh_min(1.0E30),
|
||||
m_p0(OneAtm),
|
||||
m_kk(0)
|
||||
{
|
||||
m_tlow_max = 0.0;
|
||||
m_thigh_min = 1.0E30;
|
||||
}
|
||||
|
||||
GeneralSpeciesThermo::
|
||||
GeneralSpeciesThermo(const GeneralSpeciesThermo &b) :
|
||||
m_tlow_max(b.m_tlow_max),
|
||||
m_thigh_min(b.m_thigh_min),
|
||||
m_kk(b.m_kk) {
|
||||
m_sp = b.m_sp;
|
||||
}
|
||||
|
||||
const GeneralSpeciesThermo&
|
||||
GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo &b) {
|
||||
if (&b != this) {
|
||||
m_tlow_max = b.m_tlow_max;
|
||||
m_thigh_min = b.m_thigh_min;
|
||||
m_kk = b.m_kk;
|
||||
m_sp = b.m_sp;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
GeneralSpeciesThermo::~GeneralSpeciesThermo() {
|
||||
for (int k = 0; k < m_kk; k++) {
|
||||
SpeciesThermoInterpType *sp = m_sp[k];
|
||||
if (sp) {
|
||||
delete (sp);
|
||||
m_sp[k] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
SpeciesThermo *
|
||||
GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const {
|
||||
GeneralSpeciesThermo *gsth = new GeneralSpeciesThermo(*this);
|
||||
return (SpeciesThermo *) gsth;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Install parameterization for a species.
|
||||
* @param index Species index
|
||||
* @param type ignored, since only NASA type is supported
|
||||
* @param c coefficients. These are
|
||||
* - c[0] midpoint temperature
|
||||
* - c[1] - c[7] coefficients for low T range
|
||||
* - c[8] - c[14] coefficients for high T range
|
||||
*/
|
||||
void GeneralSpeciesThermo::install(string name,
|
||||
int index,
|
||||
int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp,
|
||||
doublereal maxTemp,
|
||||
doublereal refPressure) {
|
||||
/*
|
||||
* Resize the arrays if necessary, filling the empty
|
||||
* slots with the zero pointer.
|
||||
*/
|
||||
if (index > m_kk - 1) {
|
||||
m_sp.resize(index+1, 0);
|
||||
m_kk = index+1;
|
||||
}
|
||||
|
||||
/*
|
||||
* Create the necessary object
|
||||
*/
|
||||
switch (type) {
|
||||
case NASA1:
|
||||
m_sp[index] = new NasaPoly1(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
case SHOMATE1:
|
||||
m_sp[index] = new ShomatePoly(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
case CONSTANT_CP:
|
||||
case SIMPLE:
|
||||
m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
case MU0_INTERP:
|
||||
m_sp[index] = new Mu0Poly(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
case SHOMATE2:
|
||||
m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
case NASA2:
|
||||
m_sp[index] = new NasaPoly2(index, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
break;
|
||||
default:
|
||||
throw UnknownSpeciesThermoModel(
|
||||
"GeneralSpeciesThermo::install",
|
||||
"unknown species type", int2str(type));
|
||||
break;
|
||||
}
|
||||
m_tlow_max = max(minTemp, m_tlow_max);
|
||||
m_thigh_min = min(maxTemp, m_thigh_min);
|
||||
}
|
||||
|
||||
/**
|
||||
* Update the properties for all species;
|
||||
*/
|
||||
void GeneralSpeciesThermo::
|
||||
update_one(int k, doublereal t, doublereal* cp_R,
|
||||
doublereal* h_RT, doublereal* s_R) const {
|
||||
SpeciesThermoInterpType * sp_ptr = m_sp[k];
|
||||
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Update the properties for all species;
|
||||
*/
|
||||
void GeneralSpeciesThermo::
|
||||
update(doublereal t, doublereal* cp_R,
|
||||
doublereal* h_RT, doublereal* s_R) const {
|
||||
vector<SpeciesThermoInterpType *>::const_iterator _begin, _end;
|
||||
_begin = m_sp.begin();
|
||||
_end = m_sp.end();
|
||||
SpeciesThermoInterpType * sp_ptr;
|
||||
for (; _begin != _end; ++_begin) {
|
||||
sp_ptr = *(_begin);
|
||||
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* This utility function reports the type of parameterization
|
||||
* used for the species, index.
|
||||
*/
|
||||
int GeneralSpeciesThermo::reportType(int index) const {
|
||||
SpeciesThermoInterpType *sp = m_sp[index];
|
||||
return sp->reportType();
|
||||
}
|
||||
|
||||
/**
|
||||
* This utility function reports back the type of
|
||||
* parameterization and all of the parameters for the
|
||||
* species, index.
|
||||
* For the NASA object, there are 15 coefficients.
|
||||
*/
|
||||
void GeneralSpeciesThermo::
|
||||
reportParams(int index, int &type,
|
||||
doublereal * const c,
|
||||
doublereal &minTemp,
|
||||
doublereal &maxTemp,
|
||||
doublereal &refPressure) {
|
||||
SpeciesThermoInterpType *sp = m_sp[index];
|
||||
int n;
|
||||
sp->reportParameters(n, type, minTemp, maxTemp,
|
||||
refPressure, c);
|
||||
if (n != index) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the lowest temperature at which the thermodynamic
|
||||
* parameterization is valid. If no argument is supplied, the
|
||||
* value is the one for which all species parameterizations
|
||||
* are valid. Otherwise, if an integer argument is given, the
|
||||
* value applies only to the species with that index.
|
||||
*/
|
||||
doublereal GeneralSpeciesThermo::minTemp(int k) const {
|
||||
if (k < 0)
|
||||
return m_tlow_max;
|
||||
else {
|
||||
SpeciesThermoInterpType *sp = m_sp[k];
|
||||
return sp->minTemp();
|
||||
}
|
||||
}
|
||||
|
||||
doublereal GeneralSpeciesThermo::maxTemp(int k) const {
|
||||
if (k < 0) {
|
||||
return m_thigh_min;
|
||||
} else {
|
||||
SpeciesThermoInterpType *sp = m_sp[k];
|
||||
return sp->maxTemp();
|
||||
}
|
||||
}
|
||||
|
||||
doublereal GeneralSpeciesThermo::refPressure(int k) const {
|
||||
if (k < 0) {
|
||||
return m_p0;
|
||||
} else {
|
||||
SpeciesThermoInterpType *sp = m_sp[k];
|
||||
return sp->refPressure();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
123
Cantera/src/GeneralSpeciesThermo.h
Normal file
123
Cantera/src/GeneralSpeciesThermo.h
Normal file
|
|
@ -0,0 +1,123 @@
|
|||
/**
|
||||
* @file GeneralSpeciesThermo.h
|
||||
*/
|
||||
|
||||
/*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
|
||||
#ifndef CT_GENERALSPECIESTHERMO_H
|
||||
#define CT_GENERALSPECIESTHERMO_H
|
||||
#include <string>
|
||||
#include "ct_defs.h"
|
||||
#include "SpeciesThermoMgr.h"
|
||||
#include "NasaPoly1.h"
|
||||
#include "speciesThermoTypes.h"
|
||||
#include "polyfit.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* A species thermodynamic property manager for a phase.
|
||||
* This is a general manager that can handle a wide variety
|
||||
* of species thermodynamic polynomials for individual species.
|
||||
* It is slow, however.
|
||||
*
|
||||
*
|
||||
*/
|
||||
class GeneralSpeciesThermo : public SpeciesThermo {
|
||||
|
||||
public:
|
||||
|
||||
GeneralSpeciesThermo();
|
||||
GeneralSpeciesThermo(const GeneralSpeciesThermo &);
|
||||
const GeneralSpeciesThermo & operator=(const GeneralSpeciesThermo &);
|
||||
virtual ~GeneralSpeciesThermo();
|
||||
virtual SpeciesThermo *duplMyselfAsSpeciesThermo() const ;
|
||||
|
||||
/**
|
||||
* Install parameterization for a species.
|
||||
* @param index Species index
|
||||
* @param type ignored, since only NASA type is supported
|
||||
* @param c coefficients. These are
|
||||
* - c[0] midpoint temperature
|
||||
* - c[1] - c[7] coefficients for low T range
|
||||
* - c[8] - c[14] coefficients for high T range
|
||||
*/
|
||||
virtual void install(string name, int index, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp,
|
||||
doublereal refPressure);
|
||||
/**
|
||||
* update the properties for only one species.
|
||||
*/
|
||||
virtual void update_one(int k, doublereal t, doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const;
|
||||
|
||||
virtual void update(doublereal t, doublereal* cp_R,
|
||||
doublereal* h_RT, doublereal* s_R) const;
|
||||
|
||||
/**
|
||||
* Return the lowest temperature at which the thermodynamic
|
||||
* parameterization is valid. If no argument is supplied, the
|
||||
* value is the one for which all species parameterizations
|
||||
* are valid. Otherwise, if an integer argument is given, the
|
||||
* value applies only to the species with that index.
|
||||
*/
|
||||
virtual doublereal minTemp(int k=-1) const;
|
||||
|
||||
virtual doublereal maxTemp(int k=-1) const;
|
||||
|
||||
virtual doublereal refPressure(int k = -1) const;
|
||||
|
||||
/**
|
||||
* This utility function reports the type of parameterization
|
||||
* used for the species, index.
|
||||
*/
|
||||
virtual int reportType(int index) const;
|
||||
|
||||
/**
|
||||
* This utility function reports back the type of
|
||||
* parameterization and all of the parameters for the
|
||||
* species, index.
|
||||
* For the NASA object, there are 15 coefficients.
|
||||
*/
|
||||
virtual void reportParams(int index, int &type,
|
||||
doublereal * const c,
|
||||
doublereal &minTemp,
|
||||
doublereal &maxTemp,
|
||||
doublereal &refPressure);
|
||||
|
||||
protected:
|
||||
|
||||
/**
|
||||
* This is the main unknown in the object. It is
|
||||
* a list of pointers to type SpeciesThermoInterpType.
|
||||
* Note, this object owns the objects, so they are deleted
|
||||
* in the destructor of this object.
|
||||
*/
|
||||
vector<SpeciesThermoInterpType *> m_sp;
|
||||
doublereal m_tlow_max;
|
||||
doublereal m_thigh_min;
|
||||
doublereal m_p0;
|
||||
|
||||
/**
|
||||
* Internal variable indicating the length of the
|
||||
* number of species in the phase.
|
||||
*/
|
||||
int m_kk;
|
||||
|
||||
private:
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -39,13 +39,16 @@ BASE = $(BASE_OBJ)
|
|||
|
||||
# thermodynamic properties
|
||||
THERMO_OBJ = ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o \
|
||||
SpeciesThermoFactory.o \
|
||||
SpeciesThermoFactory.o ConstCpPoly.o \
|
||||
Mu0Poly.o GeneralSpeciesThermo.o \
|
||||
ThermoFactory.o @phase_object_files@
|
||||
|
||||
THERMO_H = ThermoPhase.h IdealGasPhase.h ConstDensityThermo.h \
|
||||
SpeciesThermoFactory.h \
|
||||
ThermoFactory.h NasaPoly1.h NasaThermo.h ShomateThermo.h \
|
||||
ShomatePoly.h SimpleThermo.h SpeciesThermoMgr.h \
|
||||
SpeciesThermoFactory.h ThermoFactory.h \
|
||||
NasaPoly1.h NasaPoly2.h NasaThermo.h \
|
||||
ShomateThermo.h ShomatePoly.h ConstCpPoly.h\
|
||||
SimpleThermo.h SpeciesThermoMgr.h \
|
||||
GeneralSpeciesThermo.h Mu0Poly.h \
|
||||
speciesThermoTypes.h SpeciesThermo.h SurfPhase.h \
|
||||
EdgePhase.h polyfit.h Func1.h \
|
||||
FuncEval.h StoichManager.h @phase_header_files@
|
||||
|
|
|
|||
409
Cantera/src/Mu0Poly.cpp
Normal file
409
Cantera/src/Mu0Poly.cpp
Normal file
|
|
@ -0,0 +1,409 @@
|
|||
/**
|
||||
* @file Mu0Poly.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
|
||||
#include "Mu0Poly.h"
|
||||
#include "ctexceptions.h"
|
||||
#include "speciesThermoTypes.h"
|
||||
#include "SpeciesThermo.h"
|
||||
#include "xml.h"
|
||||
#include "ctml.h"
|
||||
|
||||
using namespace ctml;
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* The Mu0Poly class implements a linear interpolation
|
||||
* of the standard state chemical potential of one
|
||||
* species at a single reference pressure.
|
||||
* The chemical potential is input as a series of (T, mu0)
|
||||
* values. The first temperature is assumed to be equal
|
||||
* to 298.15 K; however, this may be relaxed in the future.
|
||||
* This information, and an assumption of a constant
|
||||
* heat capacity within each interval is enough to
|
||||
* calculate all thermodynamic functions.
|
||||
*
|
||||
* The basic equation for going from point 1 to point 2
|
||||
* are as follows for T, T1 <= T <= T2
|
||||
*
|
||||
* mu1 = H1 - T1 * S1
|
||||
*
|
||||
* mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1)
|
||||
*
|
||||
* S2 = S1 + Cp1(ln(T2/T1))
|
||||
*
|
||||
* H2 = H1 + Cp1(T2 - T1)
|
||||
*
|
||||
* In the future, a better assumption about the heat
|
||||
* capacity may be employed, so that it can be continuous.
|
||||
*
|
||||
* Notes about temperature interpolation for T < T1 and T > Tn
|
||||
* These are achieved by assuming a constant heat capacity
|
||||
* equal to the value in the closest temperature interval.
|
||||
* No error is thrown.
|
||||
*/
|
||||
Mu0Poly::Mu0Poly() : m_numIntervals(0),
|
||||
m_H298(0.0),
|
||||
m_lowT(0.0),
|
||||
m_highT(0.0),
|
||||
m_Pref(0.0),
|
||||
m_index(0) {
|
||||
}
|
||||
|
||||
/**
|
||||
* Mu0Poly():
|
||||
*
|
||||
* In the constructor, we calculate and store the
|
||||
* piecewise linear approximation to the thermodynamic
|
||||
* functions.
|
||||
*
|
||||
* coeffs[0] = number of points (integer)
|
||||
* 1 = H298(J/kmol)
|
||||
* 2 = T1 (Kelvin)
|
||||
* 3 = mu1 (J/kmol)
|
||||
* 4 = T2 (Kelvin)
|
||||
* 5 = mu2 (J/kmol)
|
||||
* 6 = T3 (Kelvin)
|
||||
* 7 = mu3 (J/kmol)
|
||||
* ........
|
||||
*/
|
||||
Mu0Poly::Mu0Poly(int n, doublereal tlow, doublereal thigh,
|
||||
doublereal pref,
|
||||
const doublereal* coeffs) :
|
||||
m_numIntervals(0),
|
||||
m_H298(0.0),
|
||||
m_lowT (tlow),
|
||||
m_highT (thigh),
|
||||
m_Pref (pref),
|
||||
m_index (n) {
|
||||
|
||||
int i, iindex;
|
||||
double T1, T2;
|
||||
int nPoints = (int) coeffs[0];
|
||||
if (nPoints < 2) {
|
||||
throw CanteraError("Mu0Poly",
|
||||
"nPoints must be >= 2");
|
||||
}
|
||||
m_numIntervals = nPoints - 1;
|
||||
m_H298 = coeffs[1] / GasConstant;
|
||||
int iT298 = 0;
|
||||
/*
|
||||
* Resize according to the number of points
|
||||
*/
|
||||
m_t0_int.resize(nPoints);
|
||||
m_h0_R_int.resize(nPoints);
|
||||
m_s0_R_int.resize(nPoints);
|
||||
m_cp0_R_int.resize(nPoints);
|
||||
m_mu0_R_int.resize(nPoints);
|
||||
/*
|
||||
* Calculate the T298 interval and make sure that
|
||||
* the temperatures are strictly monotonic.
|
||||
* Also distribute the data into the internal arrays.
|
||||
*/
|
||||
bool ifound = false;
|
||||
for (i = 0, iindex = 2; i < nPoints; i++) {
|
||||
T1 = coeffs[iindex];
|
||||
m_t0_int[i] = T1;
|
||||
m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant;
|
||||
if (T1 == 298.15) {
|
||||
iT298 = i;
|
||||
ifound = true;
|
||||
}
|
||||
if (i < nPoints - 1) {
|
||||
T2 = coeffs[iindex+2];
|
||||
if (T2 <= T1) {
|
||||
throw CanteraError("Mu0Poly",
|
||||
"Temperatures are not monotonic increasing");
|
||||
}
|
||||
}
|
||||
iindex += 2;
|
||||
}
|
||||
if (!ifound) {
|
||||
throw CanteraError("Mu0Poly",
|
||||
"One temperature has to be 298.15");
|
||||
}
|
||||
|
||||
/*
|
||||
* Starting from the interval with T298, we go up
|
||||
*/
|
||||
doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT;
|
||||
T1 = m_t0_int[iT298];
|
||||
doublereal mu1 = m_mu0_R_int[iT298];
|
||||
m_h0_R_int[iT298] = m_H298;
|
||||
m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1;
|
||||
for (i = iT298; i < m_numIntervals; i++) {
|
||||
T1 = m_t0_int[i];
|
||||
s1 = m_s0_R_int[i];
|
||||
h1 = m_h0_R_int[i];
|
||||
mu1 = m_mu0_R_int[i];
|
||||
T2 = m_t0_int[i+1];
|
||||
mu2 = m_mu0_R_int[i+1];
|
||||
deltaMu = mu2 - mu1;
|
||||
deltaT = T2 - T1;
|
||||
cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
|
||||
h2 = h1 + cpi * deltaT;
|
||||
s2 = s1 + cpi * log(T2/T1);
|
||||
m_cp0_R_int[i] = cpi;
|
||||
m_h0_R_int[i+1] = h2;
|
||||
m_s0_R_int[i+1] = s2;
|
||||
m_cp0_R_int[i+1] = cpi;
|
||||
}
|
||||
|
||||
/*
|
||||
* Starting from the interval with T298, we go down
|
||||
*/
|
||||
if (iT298 > 0) {
|
||||
T2 = m_t0_int[iT298];
|
||||
mu2 = m_mu0_R_int[iT298];
|
||||
m_h0_R_int[iT298] = m_H298;
|
||||
m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
|
||||
for (i = iT298 - 1; i >= 0; i--) {
|
||||
T1 = m_t0_int[i];
|
||||
mu1 = m_mu0_R_int[i];
|
||||
T2 = m_t0_int[i+1];
|
||||
mu2 = m_mu0_R_int[i+1];
|
||||
s2 = m_s0_R_int[i+1];
|
||||
h2 = m_h0_R_int[i+1];
|
||||
deltaMu = mu2 - mu1;
|
||||
deltaT = T2 - T1;
|
||||
cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
|
||||
h1 = h2 - cpi * deltaT;
|
||||
s1 = s2 - cpi * log(T2/T1);
|
||||
m_cp0_R_int[i] = cpi;
|
||||
m_h0_R_int[i] = h1;
|
||||
m_s0_R_int[i] = s1;
|
||||
if (i == (m_numIntervals-1)) {
|
||||
m_cp0_R_int[i+1] = cpi;
|
||||
}
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_HKM_NOT
|
||||
printf(" Temp mu0(J/kmol) cp0(J/kmol/K) "
|
||||
" h0(J/kmol) s0(J/kmol/K) \n");
|
||||
for (i = 0; i < nPoints; i++) {
|
||||
printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n",
|
||||
m_t0_int[i], m_mu0_R_int[i] * GasConstant,
|
||||
m_cp0_R_int[i]* GasConstant,
|
||||
m_h0_R_int[i]* GasConstant,
|
||||
m_s0_R_int[i]* GasConstant);
|
||||
fflush(stdout);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
Mu0Poly::Mu0Poly(const Mu0Poly &b)
|
||||
: m_numIntervals (b.m_numIntervals),
|
||||
m_H298 (b.m_H298),
|
||||
m_t0_int (b.m_t0_int),
|
||||
m_mu0_R_int (b.m_mu0_R_int),
|
||||
m_h0_R_int (b.m_h0_R_int),
|
||||
m_s0_R_int (b.m_s0_R_int),
|
||||
m_cp0_R_int (b.m_cp0_R_int),
|
||||
m_lowT (b.m_lowT),
|
||||
m_highT (b.m_highT),
|
||||
m_Pref (b.m_Pref),
|
||||
m_index (b.m_index) {
|
||||
}
|
||||
|
||||
Mu0Poly& Mu0Poly::operator=(const Mu0Poly& b) {
|
||||
if (&b != this) {
|
||||
m_numIntervals = b.m_numIntervals;
|
||||
m_H298 = b.m_H298;
|
||||
m_t0_int = b.m_t0_int;
|
||||
m_mu0_R_int = b.m_mu0_R_int;
|
||||
m_h0_R_int = b.m_h0_R_int;
|
||||
m_s0_R_int = b.m_s0_R_int;
|
||||
m_cp0_R_int = b.m_cp0_R_int;
|
||||
m_lowT = b.m_lowT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
/**
|
||||
* Destructor:
|
||||
*/
|
||||
Mu0Poly::~Mu0Poly(){
|
||||
}
|
||||
|
||||
SpeciesThermoInterpType *
|
||||
Mu0Poly::duplMyselfAsSpeciesThermoInterpType() const {
|
||||
Mu0Poly* mp = new Mu0Poly(*this);
|
||||
return (SpeciesThermoInterpType *) mp;
|
||||
}
|
||||
|
||||
doublereal Mu0Poly::minTemp() const { return m_lowT;}
|
||||
doublereal Mu0Poly::maxTemp() const { return m_highT;}
|
||||
doublereal Mu0Poly::refPressure() const { return m_Pref; }
|
||||
|
||||
/**
|
||||
* updateProperties is the main workhorse program.
|
||||
* Given a temperature (*tt), it calculates the thermodynamic
|
||||
* functions H/RT, S_R, and cp_R, and returns the answer.
|
||||
*
|
||||
* Note, it returns an answer by inserting the values into the
|
||||
* index position, m_index in vectors of H/RT, S_R, and cp_R.
|
||||
*
|
||||
*
|
||||
* Input
|
||||
* -------
|
||||
* *tt = Temperature (Kelvin)
|
||||
*
|
||||
*/
|
||||
void Mu0Poly::
|
||||
updateProperties(const doublereal* tt, doublereal* cp_R,
|
||||
doublereal* h_RT, doublereal* s_R) const {
|
||||
int j = m_numIntervals;
|
||||
double T = *tt;
|
||||
for (int i = 0; i < m_numIntervals; i++) {
|
||||
double T2 = m_t0_int[i+1];
|
||||
if (T <=T2) {
|
||||
j = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
double T1 = m_t0_int[j];
|
||||
double cp_Rj = m_cp0_R_int[j];
|
||||
|
||||
doublereal rt = 1.0/T;
|
||||
cp_R[m_index] = cp_Rj;
|
||||
h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj);
|
||||
s_R[m_index] = m_s0_R_int[j] + cp_Rj * (log(T/T1));
|
||||
}
|
||||
|
||||
void Mu0Poly::
|
||||
updatePropertiesTemp(const doublereal T,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
updateProperties(&T, cp_R, h_RT, s_R);
|
||||
}
|
||||
|
||||
/**
|
||||
* report all of the parameters that make up this
|
||||
* interpolation.
|
||||
*/
|
||||
void Mu0Poly::reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = MU0_INTERP;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
coeffs[0] = m_numIntervals+1;
|
||||
coeffs[1] = m_H298 * GasConstant;
|
||||
int j = 2;
|
||||
for (int i = 0; i < m_numIntervals+1; i++) {
|
||||
coeffs[j] = m_t0_int[i];
|
||||
coeffs[j+1] = m_mu0_R_int[i] * GasConstant;
|
||||
j += 2;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Install a Mu0 polynomial thermodynamic reference state property
|
||||
* parameterization for species k into a SpeciesThermo instance,
|
||||
* getting the information from an XML database.
|
||||
*/
|
||||
void installMu0ThermoFromXML(string speciesName,
|
||||
SpeciesThermo& sp, int k,
|
||||
const XML_Node* Mu0Node_ptr) {
|
||||
|
||||
doublereal tmin, tmax;
|
||||
bool dimensionlessMu0Values = false;
|
||||
const XML_Node& Mu0Node = *Mu0Node_ptr;
|
||||
|
||||
tmin = fpValue(Mu0Node["Tmin"]);
|
||||
tmax = fpValue(Mu0Node["Tmax"]);
|
||||
doublereal pref = fpValue(Mu0Node["Pref"]);
|
||||
|
||||
doublereal h298 = 0.0;
|
||||
if (Mu0Node.hasChild("H298")) {
|
||||
h298 = getFloat(Mu0Node, "H298", "actEnergy");
|
||||
}
|
||||
|
||||
int numPoints = 1;
|
||||
if (Mu0Node.hasChild("numPoints")) {
|
||||
numPoints = getInteger(Mu0Node, "numPoints");
|
||||
}
|
||||
|
||||
vector_fp cValues(numPoints);
|
||||
const XML_Node *valNode_ptr =
|
||||
getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Values");
|
||||
if (!valNode_ptr) {
|
||||
throw CanteraError("installMu0ThermoFromXML",
|
||||
"missing required while processing "
|
||||
+ speciesName);
|
||||
}
|
||||
getFloatArray(*valNode_ptr, cValues, true, "actEnergy");
|
||||
/*
|
||||
* Check to see whether the Mu0's were input in a dimensionless
|
||||
* form. If they were, then the assumed temperature needs to be
|
||||
* adjusted from the assumed T = 273.15
|
||||
*/
|
||||
string uuu = (*valNode_ptr)["units"];
|
||||
if (uuu == "Dimensionless") {
|
||||
dimensionlessMu0Values = true;
|
||||
}
|
||||
int ns = cValues.size();
|
||||
if (ns != numPoints) {
|
||||
throw CanteraError("installMu0ThermoFromXML",
|
||||
"numPoints inconsistent while processing "
|
||||
+ speciesName);
|
||||
}
|
||||
|
||||
vector_fp cTemperatures(numPoints);
|
||||
const XML_Node *tempNode_ptr =
|
||||
getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Temperatures");
|
||||
if (!tempNode_ptr) {
|
||||
throw CanteraError("installMu0ThermoFromXML",
|
||||
"missing required while processing + "
|
||||
+ speciesName);
|
||||
}
|
||||
getFloatArray(*tempNode_ptr, cTemperatures, false);
|
||||
ns = cTemperatures.size();
|
||||
if (ns != numPoints) {
|
||||
throw CanteraError("installMu0ThermoFromXML",
|
||||
"numPoints inconsistent while processing "
|
||||
+ speciesName);
|
||||
}
|
||||
|
||||
/*
|
||||
* Fix up dimensionless Mu0 values if input
|
||||
*/
|
||||
if (dimensionlessMu0Values) {
|
||||
for (int i = 0; i < numPoints; i++) {
|
||||
cValues[i] *= cTemperatures[i] / 273.15;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
vector_fp c(2 + 2 * numPoints);
|
||||
|
||||
c[0] = numPoints;
|
||||
c[1] = h298;
|
||||
for (int i = 0; i < numPoints; i++) {
|
||||
c[2+i*2] = cTemperatures[i];
|
||||
c[2+i*2+1] = cValues[i];
|
||||
}
|
||||
|
||||
sp.install(speciesName, k, MU0_INTERP, c.begin(), tmin, tmax, pref);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
136
Cantera/src/Mu0Poly.h
Normal file
136
Cantera/src/Mu0Poly.h
Normal file
|
|
@ -0,0 +1,136 @@
|
|||
/**
|
||||
* @file Mu0Poly.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
|
||||
|
||||
#ifndef CT_MU0POLY_H
|
||||
#define CT_MU0POLY_H
|
||||
|
||||
#include "SpeciesThermoInterpType.h"
|
||||
|
||||
namespace Cantera {
|
||||
class SpeciesThermo;
|
||||
class XML_Node;
|
||||
/**
|
||||
* The Mu0Poly class implements a linear interpolation
|
||||
* of the standard state chemical potential of one
|
||||
* species at a single reference pressure.
|
||||
* The chemical potential is input as a series of (T, mu0)
|
||||
* values. The first temperature is assumed to be equal
|
||||
* to 298.15 K; however, this may be relaxed in the future.
|
||||
* This information, and an assumption of a constant
|
||||
* heat capacity within each interval is enough to
|
||||
* calculate all thermodynamic functions.
|
||||
*
|
||||
* The basic equation for going from point 1 to point 2
|
||||
* are as follows for T, T1 <= T <= T2
|
||||
*
|
||||
* mu1 = H1 - T1 * S1
|
||||
*
|
||||
* mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1)
|
||||
*
|
||||
* S2 = S1 + Cp1(ln(T2/T1))
|
||||
*
|
||||
* H2 = H1 + Cp1(T2 - T1)
|
||||
*
|
||||
* In the future, a better assumption about the heat
|
||||
* capacity may be employed, so that it can be continuous.
|
||||
*
|
||||
* Notes about temperature interpolation for T < T1 and T > Tn
|
||||
* These are achieved by assuming a constant heat capacity
|
||||
* equal to the value in the closest temperature interval.
|
||||
* No error is thrown.
|
||||
*/
|
||||
class Mu0Poly: public SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
Mu0Poly();
|
||||
|
||||
Mu0Poly(int n, doublereal tlow, doublereal thigh,
|
||||
doublereal pref, const doublereal* coeffs);
|
||||
Mu0Poly(const Mu0Poly &);
|
||||
Mu0Poly& operator=(const Mu0Poly&);
|
||||
virtual ~Mu0Poly();
|
||||
SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const;
|
||||
|
||||
doublereal minTemp() const;
|
||||
doublereal maxTemp() const;
|
||||
doublereal refPressure() const;
|
||||
virtual int reportType() const { return MU0_INTERP; }
|
||||
|
||||
/**
|
||||
* Update all of the properties, using the polynomial
|
||||
* tPoly[]
|
||||
*
|
||||
* tPoly[0] = temp (Kelvin)
|
||||
*/
|
||||
void updateProperties(const doublereal* tPoly,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const ;
|
||||
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const ;
|
||||
|
||||
/**
|
||||
* report all of the parameters that make up this
|
||||
* interpolation.
|
||||
*/
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const;
|
||||
protected:
|
||||
/**
|
||||
* Number of intervals in the interpolating linear
|
||||
* approximation. Number of points is one more than the
|
||||
* number of intervals.
|
||||
*/
|
||||
int m_numIntervals;
|
||||
/**
|
||||
* Value of the enthalpy at T = 298.15.
|
||||
* This value is tied to the Heat of formation of
|
||||
* the species at 298.15.
|
||||
*/
|
||||
doublereal m_H298;
|
||||
/**
|
||||
* Points at which the standard state chemical potential
|
||||
* are given.
|
||||
*/
|
||||
vector_fp m_t0_int;
|
||||
|
||||
/*
|
||||
* Mu0's are primary input data. They aren't strictly
|
||||
* needed, but are kept here for convenience.
|
||||
*/
|
||||
vector_fp m_mu0_R_int;
|
||||
vector_fp m_h0_R_int;
|
||||
vector_fp m_s0_R_int;
|
||||
vector_fp m_cp0_R_int;
|
||||
doublereal m_lowT, m_highT, m_Pref;
|
||||
int m_index;
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
void installMu0ThermoFromXML(string speciesName,
|
||||
SpeciesThermo& sp, int k,
|
||||
const XML_Node* Mu0Node_ptr);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -13,6 +13,8 @@
|
|||
#ifndef CT_NASAPOLY1_H
|
||||
#define CT_NASAPOLY1_H
|
||||
|
||||
#include "SpeciesThermoInterpType.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
|
|
@ -41,7 +43,7 @@ namespace Cantera {
|
|||
* This class is designed specifically for use by class NasaThermo.
|
||||
* @ingroup spthermo
|
||||
*/
|
||||
class NasaPoly1 {
|
||||
class NasaPoly1 : public SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
|
|
@ -59,20 +61,59 @@ namespace Cantera {
|
|||
copy(coeffs, coeffs + 7, m_coeff.begin());
|
||||
}
|
||||
|
||||
~NasaPoly1(){}
|
||||
NasaPoly1(const NasaPoly1& b) :
|
||||
m_lowT (b.m_lowT),
|
||||
m_highT (b.m_highT),
|
||||
m_Pref (b.m_Pref),
|
||||
m_index (b.m_index),
|
||||
m_coeff (array_fp(7)) {
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 7,
|
||||
m_coeff.begin());
|
||||
}
|
||||
|
||||
NasaPoly1& operator=(const NasaPoly1& b) {
|
||||
if (&b != this) {
|
||||
m_lowT = b.m_lowT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 7,
|
||||
m_coeff.begin());
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
virtual ~NasaPoly1(){}
|
||||
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const {
|
||||
NasaPoly1* np = new NasaPoly1(*this);
|
||||
return (SpeciesThermoInterpType *) np;
|
||||
}
|
||||
|
||||
doublereal minTemp() const { return m_lowT;}
|
||||
doublereal maxTemp() const { return m_highT;}
|
||||
doublereal refPressure() const { return m_Pref; }
|
||||
virtual int reportType() const { return NASA1; }
|
||||
|
||||
/**
|
||||
* Update the properties for this species. This method is called
|
||||
* with a pointer to an array containing the functions of
|
||||
* temperature needed by this
|
||||
* temperature needed by this
|
||||
* parameterization, and three pointers to arrays where the
|
||||
* computed property values
|
||||
* computed property values
|
||||
* should be written. This method updates only one value in
|
||||
* each array.
|
||||
*
|
||||
* Temperature Polynomial:
|
||||
* tt[0] = t;
|
||||
* tt[1] = t*t;
|
||||
* tt[2] = m_t[1]*t;
|
||||
* tt[3] = m_t[2]*t;
|
||||
* tt[4] = 1.0/t;
|
||||
* tt[5] = log(t);
|
||||
*/
|
||||
void updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
|
||||
|
|
@ -97,10 +138,32 @@ namespace Cantera {
|
|||
s_R[m_index] = s;
|
||||
}
|
||||
|
||||
void reportParameters(int &n, doublereal &tlow, doublereal &thigh,
|
||||
/**
|
||||
* updatePropertiesTemp():
|
||||
* This formulation creates its own temperature
|
||||
* polynomial. Then, it calls updateProperties();
|
||||
*
|
||||
* (note: this is slow, but it is general)
|
||||
*/
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
double tPoly[6];
|
||||
tPoly[0] = temp;
|
||||
tPoly[1] = temp * temp;
|
||||
tPoly[2] = tPoly[1] * temp;
|
||||
tPoly[3] = tPoly[2] * temp;
|
||||
tPoly[4] = 1.0 / temp;
|
||||
tPoly[5] = log(temp);
|
||||
updateProperties(tPoly, cp_R, h_RT, s_R);
|
||||
}
|
||||
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = NASA1;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
|
|
@ -122,11 +185,5 @@ namespace Cantera {
|
|||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
217
Cantera/src/NasaPoly2.h
Normal file
217
Cantera/src/NasaPoly2.h
Normal file
|
|
@ -0,0 +1,217 @@
|
|||
/**
|
||||
* @file NasaPoly1.h
|
||||
*/
|
||||
|
||||
/* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
|
||||
#ifndef CT_NASAPOLY2_H
|
||||
#define CT_NASAPOLY2_H
|
||||
|
||||
#include "SpeciesThermoInterpType.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
*
|
||||
*
|
||||
* The NASA polynomial parameterization for one temperature range.
|
||||
* This parameterization expresses the heat capacity as a
|
||||
* fourth-order polynomial. Note that this is the form used in the
|
||||
* 1971 NASA equilibrium program and by the Chemkin software
|
||||
* package, but differs from the form used in the more recent NASA
|
||||
* equilibrium program.
|
||||
*
|
||||
* Seven coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
|
||||
* \f$ c_p^0(T)\f$, \f$ h^0(T)\f$, and \f$ s^0(T) \f$ as
|
||||
* polynomials in \f$ T \f$ :
|
||||
* \f[
|
||||
* \frac{c_p(T)}{R} = a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4
|
||||
* \f]
|
||||
* \f[
|
||||
* \frac{h^0(T)}{RT} = a_0 + \frac{a_1}{2} T + \frac{a_2}{3} T^2
|
||||
* + \frac{a_3}{4} T^3 + \frac{a_4}{5} T^4 + \frac{a_5}{T}.
|
||||
* \f]
|
||||
* \f[
|
||||
* \frac{s^0(T)}{R} = a_0\ln T + a_1 T + \frac{a_2}{2} T^2
|
||||
+ \frac{a_3}{3} T^3 + \frac{a_4}{4} T^4 + a_6.
|
||||
* \f]
|
||||
*
|
||||
* This class is designed specifically for use by class
|
||||
* GeneralSpeciesThermo.
|
||||
* @ingroup spthermo
|
||||
*/
|
||||
class NasaPoly2 : public SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
NasaPoly2()
|
||||
: m_lowT(0.0),
|
||||
m_midT(0.0),
|
||||
m_highT (0.0),
|
||||
m_Pref(0.0),
|
||||
mnp_low(0),
|
||||
mnp_high(0),
|
||||
m_index(0),
|
||||
m_coeff(array_fp(15)) {
|
||||
}
|
||||
|
||||
NasaPoly2(int n, doublereal tlow, doublereal thigh, doublereal pref,
|
||||
const doublereal* coeffs) :
|
||||
m_lowT(tlow),
|
||||
m_highT(thigh),
|
||||
m_Pref(pref),
|
||||
mnp_low(0),
|
||||
mnp_high(0),
|
||||
m_index(n),
|
||||
m_coeff(array_fp(15)) {
|
||||
|
||||
copy(coeffs, coeffs + 15, m_coeff.begin());
|
||||
m_midT = coeffs[0];
|
||||
mnp_low = new NasaPoly1(m_index, m_lowT, m_midT,
|
||||
m_Pref, m_coeff.begin()+1);
|
||||
mnp_high = new NasaPoly1(m_index, m_midT, m_highT,
|
||||
m_Pref, m_coeff.begin()+8);
|
||||
}
|
||||
|
||||
NasaPoly2(const NasaPoly2& b) :
|
||||
m_lowT(b.m_lowT),
|
||||
m_midT(b.m_midT),
|
||||
m_highT(b.m_highT),
|
||||
m_Pref(b.m_Pref),
|
||||
mnp_low(0),
|
||||
mnp_high(0),
|
||||
m_index(b.m_index),
|
||||
m_coeff(array_fp(15)) {
|
||||
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 15,
|
||||
m_coeff.begin());
|
||||
mnp_low = new NasaPoly1(m_index, m_lowT, m_midT,
|
||||
m_Pref, m_coeff.begin()+1);
|
||||
mnp_high = new NasaPoly1(m_index, m_midT, m_highT,
|
||||
m_Pref, m_coeff.begin()+8);
|
||||
}
|
||||
|
||||
NasaPoly2& operator=(const NasaPoly2& b) {
|
||||
if (&b != this) {
|
||||
m_lowT = b.m_lowT;
|
||||
m_midT = b.m_midT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 15,
|
||||
m_coeff.begin());
|
||||
if (mnp_low) delete mnp_low;
|
||||
if (mnp_high) delete mnp_high;
|
||||
mnp_low = new NasaPoly1(m_index, m_lowT, m_midT,
|
||||
m_Pref, m_coeff.begin()+1);
|
||||
mnp_high = new NasaPoly1(m_index, m_midT, m_highT,
|
||||
m_Pref, m_coeff.begin()+8);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
virtual ~NasaPoly2(){
|
||||
delete mnp_low;
|
||||
delete mnp_high;
|
||||
}
|
||||
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const {
|
||||
NasaPoly2* np = new NasaPoly2(*this);
|
||||
return (SpeciesThermoInterpType *) np;
|
||||
}
|
||||
|
||||
doublereal minTemp() const { return m_lowT;}
|
||||
doublereal maxTemp() const { return m_highT;}
|
||||
doublereal refPressure() const { return m_Pref; }
|
||||
virtual int reportType() const { return NASA2; }
|
||||
|
||||
/**
|
||||
* Update the properties for this species. This method is called
|
||||
* with a pointer to an array containing the functions of
|
||||
* temperature needed by this
|
||||
* parameterization, and three pointers to arrays where the
|
||||
* computed property values
|
||||
* should be written. This method updates only one value in
|
||||
* each array.
|
||||
*
|
||||
* Temperature Polynomial:
|
||||
* tt[0] = t;
|
||||
* tt[1] = t*t;
|
||||
* tt[2] = m_t[1]*t;
|
||||
* tt[3] = m_t[2]*t;
|
||||
* tt[4] = 1.0/t;
|
||||
* tt[5] = log(t);
|
||||
*/
|
||||
void updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
|
||||
|
||||
double T = tt[0];
|
||||
if (T <= m_midT) {
|
||||
mnp_low->updateProperties(tt, cp_R, h_RT, s_R);
|
||||
} else {
|
||||
mnp_high->updateProperties(tt, cp_R, h_RT, s_R);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* updatePropertiesTemp():
|
||||
* This formulation creates its own temperature
|
||||
* polynomial. Then, it calls updateProperties();
|
||||
*
|
||||
* (note: this is slow, but it is general)
|
||||
*/
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
if (temp <= m_midT) {
|
||||
mnp_low->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
|
||||
} else {
|
||||
mnp_high->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
|
||||
}
|
||||
}
|
||||
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = NASA2;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
for (int i = 0; i < 15; i++) {
|
||||
coeffs[i] = m_coeff[i];
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
|
||||
doublereal m_lowT; // lowest valid temperature
|
||||
doublereal m_midT;
|
||||
doublereal m_highT; // highest valid temperature
|
||||
doublereal m_Pref; // standard-state pressure
|
||||
NasaPoly1 *mnp_low;
|
||||
NasaPoly1 *mnp_high;
|
||||
int m_index; // species index
|
||||
array_fp m_coeff; // array of polynomial coefficients
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -65,7 +65,8 @@ namespace Cantera {
|
|||
* - c[1] - c[7] coefficients for low T range
|
||||
* - c[8] - c[14] coefficients for high T range
|
||||
*/
|
||||
virtual void install(string name, int index, int type, const doublereal* c,
|
||||
virtual void install(string name, int index, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp,
|
||||
doublereal refPressure) {
|
||||
|
||||
|
|
@ -184,7 +185,9 @@ namespace Cantera {
|
|||
return m_thigh[k];
|
||||
}
|
||||
|
||||
virtual doublereal refPressure() const {return m_p0;}
|
||||
virtual doublereal refPressure(int k = -1) const {
|
||||
return m_p0;
|
||||
}
|
||||
|
||||
/**
|
||||
* This utility function reports the type of parameterization
|
||||
|
|
@ -211,21 +214,27 @@ namespace Cantera {
|
|||
const vector<NasaPoly1> &mhg = m_high[grp-1];
|
||||
const NasaPoly1 *lowPoly = &(mlg[pos]);
|
||||
const NasaPoly1 *highPoly = &(mhg[pos]);
|
||||
|
||||
int itype = NASA;
|
||||
doublereal tmid = lowPoly->maxTemp();
|
||||
c[0] = tmid;
|
||||
int n;
|
||||
double ttemp;
|
||||
lowPoly->reportParameters(n, minTemp, ttemp, refPressure,
|
||||
lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
|
||||
c + 1);
|
||||
if (n != index) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
highPoly->reportParameters(n, ttemp, maxTemp, refPressure,
|
||||
if (itype != NASA1) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure,
|
||||
c + 8);
|
||||
if (n != index) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
if (itype != NASA1) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
} else {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
|
|
|
|||
|
|
@ -12,6 +12,8 @@
|
|||
#ifndef CT_SHOMATEPOLY1_H
|
||||
#define CT_SHOMATEPOLY1_H
|
||||
|
||||
#include "SpeciesThermoInterpType.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
|
|
@ -32,7 +34,7 @@ namespace Cantera {
|
|||
* \f]
|
||||
*/
|
||||
|
||||
class ShomatePoly {
|
||||
class ShomatePoly : public SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
|
|
@ -50,23 +52,59 @@ namespace Cantera {
|
|||
copy(coeffs, coeffs + 7, m_coeff.begin());
|
||||
}
|
||||
|
||||
~ShomatePoly(){}
|
||||
ShomatePoly(const ShomatePoly& b) :
|
||||
m_lowT (b.m_lowT),
|
||||
m_highT (b.m_highT),
|
||||
m_Pref (b.m_Pref),
|
||||
m_coeff (array_fp(7)),
|
||||
m_index (b.m_index) {
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 7,
|
||||
m_coeff.begin());
|
||||
}
|
||||
|
||||
ShomatePoly& operator=(const ShomatePoly& b) {
|
||||
if (&b != this) {
|
||||
m_lowT = b.m_lowT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 7,
|
||||
m_coeff.begin());
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
virtual ~ShomatePoly(){}
|
||||
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const {
|
||||
ShomatePoly* sp = new ShomatePoly(*this);
|
||||
return (SpeciesThermoInterpType *) sp;
|
||||
}
|
||||
|
||||
doublereal minTemp() const { return m_lowT;}
|
||||
doublereal maxTemp() const { return m_highT;}
|
||||
doublereal refPressure() const { return m_Pref; }
|
||||
virtual int reportType() const { return SHOMATE; }
|
||||
|
||||
/**
|
||||
* t is T/1000.
|
||||
*
|
||||
* tt[0] t
|
||||
* tt[1] t*t
|
||||
* tt[2] t*t*t
|
||||
* tt[3] t^4
|
||||
* tt[4] ln t
|
||||
* This formulation calculates the thermo functions
|
||||
* given the native formulation of the temperature
|
||||
* polynomial
|
||||
*
|
||||
* tt is T/1000.
|
||||
* m_t[0] = tt;
|
||||
* m_t[1] = tt*tt;
|
||||
* m_t[2] = m_t[1]*tt;
|
||||
* m_t[3] = 1.0/m_t[1];
|
||||
* m_t[4] = log(tt);
|
||||
* m_t[5] = 1.0/GasConstant;
|
||||
* m_t[6] = 1.0/(GasConstant * T);
|
||||
*/
|
||||
void updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
|
||||
doublereal A = m_coeff[0];
|
||||
doublereal Bt = m_coeff[1]*tt[0];
|
||||
|
|
@ -80,18 +118,49 @@ namespace Cantera {
|
|||
cp = A + Bt + Ct2 + Dt3 + Etm2;
|
||||
h = tt[0]*(A + 0.5*Bt + OneThird*Ct2 + 0.25*Dt3 - Etm2) + F;
|
||||
s = A*tt[4] + Bt + 0.5*Ct2 + OneThird*Dt3 - 0.5*Etm2 + G;
|
||||
h *= 1.e6;
|
||||
|
||||
/*
|
||||
* Shomate polynomials parameterizes assuming units of
|
||||
* J/(gmol*K) for cp_r and s_R and kJ/(gmol) for h.
|
||||
* However, Cantera assumes default MKS units of
|
||||
* J/(kmol*K). This requires us to multiply cp and s
|
||||
* by 1.e3 and h by 1.e6, before we then nondimensionlize
|
||||
* the results by dividing by (GasConstant * T),
|
||||
* where GasConstant has units of J/(kmol * K).
|
||||
*/
|
||||
cp_R[m_index] = 1.e3 * cp * tt[5];
|
||||
h_RT[m_index] = h * tt[6];
|
||||
s_R[m_index] = 1.e3 * s * tt[5];
|
||||
h_RT[m_index] = 1.e6 * h * tt[6];
|
||||
s_R[m_index] = 1.e3 * s * tt[5];
|
||||
}
|
||||
|
||||
/**
|
||||
* updatePropertiesTemp():
|
||||
* This formulation creates its own temperature
|
||||
* polynomial. Then, it calls updateProperties();
|
||||
* -> general, but slow.
|
||||
*/
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
double tPoly[7];
|
||||
doublereal tt = 1.e-3*temp;
|
||||
tPoly[0] = tt;
|
||||
tPoly[1] = tt * tt;
|
||||
tPoly[2] = tPoly[1] * tt;
|
||||
tPoly[3] = 1.0/tPoly[1];
|
||||
tPoly[4] = log(tt);
|
||||
tPoly[5] = 1.0/GasConstant;
|
||||
tPoly[6] = 1.0/(GasConstant * temp);
|
||||
updateProperties(tPoly, cp_R, h_RT, s_R);
|
||||
}
|
||||
|
||||
void reportParameters(int &n, doublereal &tlow, doublereal &thigh,
|
||||
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = SHOMATE;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
|
|
@ -109,12 +178,159 @@ namespace Cantera {
|
|||
|
||||
};
|
||||
|
||||
|
||||
|
||||
class ShomatePoly2 : public SpeciesThermoInterpType {
|
||||
public:
|
||||
|
||||
ShomatePoly2()
|
||||
: m_lowT(0.0),
|
||||
m_midT(0.0),
|
||||
m_highT (0.0),
|
||||
m_Pref(0.0),
|
||||
msp_low(0),
|
||||
msp_high(0),
|
||||
m_index(0) {
|
||||
m_coeff.resize(15);
|
||||
}
|
||||
|
||||
ShomatePoly2(int n, doublereal tlow, doublereal thigh, doublereal pref,
|
||||
const doublereal* coeffs) :
|
||||
m_lowT (tlow),
|
||||
m_midT(0.0),
|
||||
m_highT (thigh),
|
||||
m_Pref (pref),
|
||||
msp_low(0),
|
||||
msp_high(0),
|
||||
m_index (n) {
|
||||
m_coeff.resize(15);
|
||||
copy(coeffs, coeffs + 15, m_coeff.begin());
|
||||
m_midT = coeffs[0];
|
||||
msp_low = new ShomatePoly(n, tlow, m_midT, pref, coeffs+1);
|
||||
msp_high = new ShomatePoly(n, m_midT, thigh, pref, coeffs+8);
|
||||
}
|
||||
|
||||
ShomatePoly2(const ShomatePoly2& b) :
|
||||
m_lowT (b.m_lowT),
|
||||
m_midT (b.m_midT),
|
||||
m_highT (b.m_highT),
|
||||
m_Pref (b.m_Pref),
|
||||
msp_low(0),
|
||||
msp_high(0),
|
||||
m_coeff (array_fp(15)),
|
||||
m_index (b.m_index) {
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 15,
|
||||
m_coeff.begin());
|
||||
msp_low = new ShomatePoly(m_index, m_lowT, m_midT,
|
||||
m_Pref, m_coeff.begin()+1);
|
||||
msp_high = new ShomatePoly(m_index, m_midT, m_highT,
|
||||
m_Pref, m_coeff.begin()+8);
|
||||
}
|
||||
|
||||
ShomatePoly2& operator=(const ShomatePoly2& b) {
|
||||
if (&b != this) {
|
||||
m_lowT = b.m_lowT;
|
||||
m_midT = b.m_midT;
|
||||
m_highT = b.m_highT;
|
||||
m_Pref = b.m_Pref;
|
||||
m_index = b.m_index;
|
||||
copy(b.m_coeff.begin(),
|
||||
b.m_coeff.begin() + 15,
|
||||
m_coeff.begin());
|
||||
if (msp_low) delete msp_low;
|
||||
if (msp_high) delete msp_high;
|
||||
msp_low = new ShomatePoly(m_index, m_lowT, m_midT,
|
||||
m_Pref, m_coeff.begin()+1);
|
||||
msp_high = new ShomatePoly(m_index, m_midT, m_highT,
|
||||
m_Pref, m_coeff.begin()+8);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
virtual ~ShomatePoly2(){
|
||||
delete msp_low;
|
||||
delete msp_high;
|
||||
}
|
||||
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const {
|
||||
ShomatePoly2* sp = new ShomatePoly2(*this);
|
||||
return (SpeciesThermoInterpType *) sp;
|
||||
}
|
||||
|
||||
doublereal minTemp() const { return m_lowT;}
|
||||
doublereal maxTemp() const { return m_highT;}
|
||||
doublereal refPressure() const { return m_Pref; }
|
||||
virtual int reportType() const { return SHOMATE2; }
|
||||
|
||||
/**
|
||||
* This formulation calculates the thermo functions
|
||||
* given the native formulation of the temperature
|
||||
* polynomial
|
||||
*
|
||||
* tt is T/1000.
|
||||
* m_t[0] = tt;
|
||||
* m_t[1] = tt*tt;
|
||||
* m_t[2] = m_t[1]*tt;
|
||||
* m_t[3] = 1.0/m_t[1];
|
||||
* m_t[4] = log(tt);
|
||||
* m_t[5] = 1.0/GasConstant;
|
||||
* m_t[6] = 1.0/(GasConstant * T);
|
||||
*/
|
||||
void updateProperties(const doublereal* tt,
|
||||
doublereal* cp_R, doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
double T = 1000 * tt[0];
|
||||
if (T <= m_midT) {
|
||||
msp_low->updateProperties(tt, cp_R, h_RT, s_R);
|
||||
} else {
|
||||
msp_high->updateProperties(tt, cp_R, h_RT, s_R);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* updatePropertiesTemp():
|
||||
* This formulation creates its own temperature
|
||||
* polynomial. Then, it calls updateProperties();
|
||||
* -> general, but slow.
|
||||
*/
|
||||
void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const {
|
||||
if (temp <= m_midT) {
|
||||
msp_low->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
|
||||
} else {
|
||||
msp_high->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
|
||||
}
|
||||
}
|
||||
|
||||
void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const {
|
||||
n = m_index;
|
||||
type = SHOMATE2;
|
||||
tlow = m_lowT;
|
||||
thigh = m_highT;
|
||||
pref = m_Pref;
|
||||
for (int i = 0; i < 15; i++) {
|
||||
coeffs[i] = m_coeff[i];
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
|
||||
doublereal m_lowT, m_midT;
|
||||
doublereal m_highT;
|
||||
doublereal m_Pref;
|
||||
ShomatePoly *msp_low;
|
||||
ShomatePoly *msp_high;
|
||||
array_fp m_coeff;
|
||||
int m_index;
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -10,12 +10,20 @@
|
|||
*
|
||||
* \f[ S = A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G \f]
|
||||
*
|
||||
* - Cp = heat capacity (J/mol*K)
|
||||
* - H = standard enthalpy (kJ/mol)
|
||||
* - \f$ \Delta_f H_298.15 \f$ = enthalpy of formation at 298.15 K (kJ/mol)
|
||||
* - S = standard entropy (J/mol*K)
|
||||
* - Cp = heat capacity (J/gmol*K)
|
||||
* - H = standard enthalpy (kJ/gmol)
|
||||
* - \f$ \Delta_f H_298.15 \f$ = enthalpy of formation at 298.15 K (kJ/gmol)
|
||||
* - S = standard entropy (J/gmol*K)
|
||||
* - t = temperature (K) / 1000.
|
||||
*
|
||||
* Note, the polynomial data (i.e., A, ... , G) is entered in dimensional form.
|
||||
* This is in contrast to the NASA database polynomials which are entered in
|
||||
* nondimensional form (i.e., NASA parameterizes C_p/R, while Shomate
|
||||
* parameterizes C_p assuming units of J/gmol*K - and kJ/gmol*K for H).
|
||||
* Note, also that the H - H_298.15 equation has units of kJ/gmol, because of
|
||||
* the implicit integration of (t = T 1000), which provides a
|
||||
* multiplier of 1000 to the Enthalpy equation.
|
||||
*
|
||||
*/
|
||||
|
||||
#ifndef CT_SHOMATETHERMO_H
|
||||
|
|
@ -54,8 +62,10 @@ namespace Cantera {
|
|||
* in the same units as used in the NIST Chemistry WebBook.
|
||||
*
|
||||
*/
|
||||
virtual void install(string name, int index, int type, const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
|
||||
virtual void install(string name, int index, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp,
|
||||
doublereal refPressure) {
|
||||
int imid = int(c[0]); // midpoint temp converted to integer
|
||||
int igrp = m_index[imid]; // has this value been seen before?
|
||||
if (igrp == 0) { // if not, prepare new group
|
||||
|
|
@ -139,7 +149,7 @@ namespace Cantera {
|
|||
_end = m_low[i].end();
|
||||
}
|
||||
for (; _begin != _end; ++_begin) {
|
||||
_begin->updateProperties(m_t.begin(), cp_R, h_RT, s_R);
|
||||
_begin->updateProperties(m_t.begin(), cp_R, h_RT, s_R);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -158,7 +168,9 @@ namespace Cantera {
|
|||
return m_thigh[k];
|
||||
}
|
||||
|
||||
virtual doublereal refPressure() const {return m_p0;}
|
||||
virtual doublereal refPressure(int k=-1) const {
|
||||
return m_p0;
|
||||
}
|
||||
|
||||
virtual int reportType(int index) const { return SHOMATE; }
|
||||
|
||||
|
|
@ -174,9 +186,10 @@ namespace Cantera {
|
|||
doublereal &maxTemp,
|
||||
doublereal &refPressure) {
|
||||
type = reportType(index);
|
||||
if (type == NASA) {
|
||||
if (type == SHOMATE) {
|
||||
int grp = m_group_map[index];
|
||||
int pos = m_posInGroup_map[index];
|
||||
int itype = SHOMATE;
|
||||
const vector<ShomatePoly> &mlg = m_low[grp-1];
|
||||
const vector<ShomatePoly> &mhg = m_high[grp-1];
|
||||
const ShomatePoly *lowPoly = &(mlg[pos]);
|
||||
|
|
@ -185,16 +198,22 @@ namespace Cantera {
|
|||
c[0] = tmid;
|
||||
int n;
|
||||
double ttemp;
|
||||
lowPoly->reportParameters(n, minTemp, ttemp, refPressure,
|
||||
lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
|
||||
c + 1);
|
||||
if (n != index) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
highPoly->reportParameters(n, ttemp, maxTemp, refPressure,
|
||||
c + 8);
|
||||
if (itype != SHOMATE && itype != SHOMATE1) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
highPoly->reportParameters(n, itype, ttemp, maxTemp,
|
||||
refPressure, c + 8);
|
||||
if (n != index) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
if (itype != SHOMATE && itype != SHOMATE1) {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
} else {
|
||||
throw CanteraError(" ", "confused");
|
||||
}
|
||||
|
|
|
|||
|
|
@ -27,7 +27,8 @@ namespace Cantera {
|
|||
|
||||
virtual ~SimpleThermo() {}
|
||||
|
||||
virtual void install(string name, int index, int type, const doublereal* c,
|
||||
virtual void install(string name, int index, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
|
||||
m_logt0.push_back(log(c[0]));
|
||||
m_t0.push_back(c[0]);
|
||||
|
|
@ -80,7 +81,7 @@ namespace Cantera {
|
|||
return m_thigh[k];
|
||||
}
|
||||
|
||||
virtual doublereal refPressure() const {return m_p0;}
|
||||
virtual doublereal refPressure(int k=-1) const {return m_p0;}
|
||||
|
||||
virtual int reportType(int index) const { return SIMPLE; }
|
||||
|
||||
|
|
|
|||
|
|
@ -64,14 +64,18 @@ namespace Cantera {
|
|||
* parameterization.
|
||||
* @see speciesThermoTypes.h
|
||||
*/
|
||||
virtual void install(string name, int index, int type, const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp, doublereal refPressure)=0;
|
||||
virtual void install(string name, int index, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp,
|
||||
doublereal maxTemp,
|
||||
doublereal refPressure)=0;
|
||||
|
||||
/**
|
||||
* Compute the standard-state properties for all species.
|
||||
* Given temperature T in K, this method updates the values of
|
||||
* the non-dimensional heat capacity at constant pressure,
|
||||
* enthalpy, and entropy.
|
||||
* enthalpy, and entropy, at the reference pressure Pref
|
||||
* of each of the standard states.
|
||||
*/
|
||||
virtual void update(doublereal T,
|
||||
doublereal* cp_R,
|
||||
|
|
@ -79,7 +83,7 @@ namespace Cantera {
|
|||
doublereal* s_R) const=0;
|
||||
|
||||
/**
|
||||
* Like update, but only updates the species k.
|
||||
* Like update(), but only updates the single species k.
|
||||
*/
|
||||
virtual void update_one(int k, doublereal T,
|
||||
doublereal* cp_R,
|
||||
|
|
@ -107,10 +111,17 @@ namespace Cantera {
|
|||
virtual doublereal maxTemp(int k=-1) const =0;
|
||||
|
||||
/**
|
||||
* The standard-state pressure. All parameterizations must be
|
||||
* for the same standard-state pressure.
|
||||
* The reference-state pressure for species k.
|
||||
*
|
||||
* returns the reference state pressure in Pascals for
|
||||
* species k. If k is left out of the argument list,
|
||||
* it returns the reference state pressure for the first
|
||||
* species.
|
||||
* Note that some SpeciesThermo implementations, such
|
||||
* as those for ideal gases, require that all species
|
||||
* in the same phase have the same reference state pressures.
|
||||
*/
|
||||
virtual doublereal refPressure() const =0;
|
||||
virtual doublereal refPressure(int k=-1) const =0;
|
||||
|
||||
/**
|
||||
* This utility function reports the type of parameterization
|
||||
|
|
|
|||
|
|
@ -16,6 +16,8 @@
|
|||
#include "ShomateThermo.h"
|
||||
//#include "PolyThermoMgr.h"
|
||||
#include "SimpleThermo.h"
|
||||
#include "GeneralSpeciesThermo.h"
|
||||
#include "Mu0Poly.h"
|
||||
|
||||
#include "SpeciesThermoMgr.h"
|
||||
#include "speciesThermoTypes.h"
|
||||
|
|
@ -31,7 +33,8 @@ namespace Cantera {
|
|||
|
||||
|
||||
static void getSpeciesThermoTypes(XML_Node* node,
|
||||
int& has_nasa, int& has_shomate, int& has_simple) {
|
||||
int& has_nasa, int& has_shomate, int& has_simple,
|
||||
int &has_other) {
|
||||
const XML_Node& sparray = *node;
|
||||
vector<XML_Node*> sp;
|
||||
sparray.getChildren("species",sp);
|
||||
|
|
@ -44,10 +47,11 @@ namespace Cantera {
|
|||
if (th.hasChild("Shomate")) has_shomate = 1;
|
||||
if (th.hasChild("const_cp")) has_simple = 1;
|
||||
if (th.hasChild("poly")) {
|
||||
if (th.child("poly")["order"] == "1") has_simple = 1;
|
||||
else throw CanteraError("newSpeciesThermo",
|
||||
"poly with order > 1 not yet supported");
|
||||
if (th.child("poly")["order"] == "1") has_simple = 1;
|
||||
else throw CanteraError("newSpeciesThermo",
|
||||
"poly with order > 1 not yet supported");
|
||||
}
|
||||
if (th.hasChild("Mu0")) has_other = 1;
|
||||
} else {
|
||||
throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
|
||||
spNode->attrib("name"), "missing");
|
||||
|
|
@ -61,8 +65,15 @@ namespace Cantera {
|
|||
* specified in a CTML phase specification.
|
||||
*/
|
||||
SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(XML_Node* node) {
|
||||
int inasa = 0, ishomate = 0, isimple = 0;
|
||||
getSpeciesThermoTypes(node, inasa, ishomate, isimple);
|
||||
int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
|
||||
try {
|
||||
getSpeciesThermoTypes(node, inasa, ishomate, isimple, iother);
|
||||
} catch (UnknownSpeciesThermoModel) {
|
||||
iother = 1;
|
||||
}
|
||||
if (iother) {
|
||||
return new GeneralSpeciesThermo();
|
||||
}
|
||||
return newSpeciesThermo(NASA*inasa
|
||||
+ SHOMATE*ishomate + SIMPLE*isimple);
|
||||
}
|
||||
|
|
@ -70,10 +81,17 @@ namespace Cantera {
|
|||
SpeciesThermo* SpeciesThermoFactory::
|
||||
newSpeciesThermo(vector<XML_Node*> nodes) {
|
||||
int n = static_cast<int>(nodes.size());
|
||||
int inasa = 0, ishomate = 0, isimple = 0;
|
||||
int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
|
||||
for (int j = 0; j < n; j++) {
|
||||
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple);
|
||||
try {
|
||||
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
|
||||
} catch (UnknownSpeciesThermoModel) {
|
||||
iother = 1;
|
||||
}
|
||||
}
|
||||
if (iother) {
|
||||
return new GeneralSpeciesThermo();
|
||||
}
|
||||
return newSpeciesThermo(NASA*inasa
|
||||
+ SHOMATE*ishomate + SIMPLE*isimple);
|
||||
}
|
||||
|
|
@ -82,14 +100,18 @@ namespace Cantera {
|
|||
SpeciesThermo* SpeciesThermoFactory::
|
||||
newSpeciesThermoOpt(vector<XML_Node*> nodes) {
|
||||
int n = static_cast<int>(nodes.size());
|
||||
int inasa = 0, ishomate = 0, isimple = 0;
|
||||
int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
|
||||
for (int j = 0; j < n; j++) {
|
||||
try {
|
||||
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple);
|
||||
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
|
||||
} catch (UnknownSpeciesThermoModel) {
|
||||
|
||||
iother = 1;
|
||||
popError();
|
||||
}
|
||||
}
|
||||
if (iother) {
|
||||
return new GeneralSpeciesThermo();
|
||||
}
|
||||
return newSpeciesThermo(NASA*inasa
|
||||
+ SHOMATE*ishomate + SIMPLE*isimple);
|
||||
}
|
||||
|
|
@ -98,7 +120,6 @@ namespace Cantera {
|
|||
SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(int type) {
|
||||
|
||||
switch (type) {
|
||||
|
||||
case NASA:
|
||||
return new NasaThermo;
|
||||
case SHOMATE:
|
||||
|
|
@ -308,8 +329,10 @@ namespace Cantera {
|
|||
getFloatArray(f0.child("floatArray"), c0, false);
|
||||
if (dualRange)
|
||||
getFloatArray(f1ptr->child("floatArray"), c1, false);
|
||||
else
|
||||
c1.resize(7,0.0);
|
||||
else {
|
||||
c1.resize(7,0.0);
|
||||
copy(c0.begin(), c0.begin()+7, c1.begin());
|
||||
}
|
||||
}
|
||||
else if (fabs(tmax1 - tmin0) < 0.01) {
|
||||
tmin = tmin1;
|
||||
|
|
@ -385,6 +408,9 @@ namespace Cantera {
|
|||
else if (f->name() == "NASA") {
|
||||
installNasaThermoFromXML(s["name"], spthermo, k, f, 0);
|
||||
}
|
||||
else if (f->name() == "Mu0") {
|
||||
installMu0ThermoFromXML(s["name"], spthermo, k, f);
|
||||
}
|
||||
else {
|
||||
throw UnknownSpeciesThermoModel("installSpecies",
|
||||
s["name"], f->name());
|
||||
|
|
|
|||
49
Cantera/src/SpeciesThermoInterpType.h
Normal file
49
Cantera/src/SpeciesThermoInterpType.h
Normal file
|
|
@ -0,0 +1,49 @@
|
|||
/**
|
||||
* @file SpeciesThermoInterpType.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
#include "speciesThermoTypes.h"
|
||||
#ifndef CT_SPECIESTHERMOINTERPTYPE_H
|
||||
#define CT_SPECIESTHERMOINTERPTYPE_H
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
|
||||
class SpeciesThermoInterpType {
|
||||
|
||||
public:
|
||||
|
||||
virtual SpeciesThermoInterpType *
|
||||
duplMyselfAsSpeciesThermoInterpType() const = 0;
|
||||
|
||||
virtual doublereal minTemp() const = 0;
|
||||
virtual doublereal maxTemp() const = 0;
|
||||
virtual doublereal refPressure() const = 0;
|
||||
virtual int reportType() const = 0;
|
||||
|
||||
virtual void updateProperties(const doublereal* tempPoly,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const = 0;
|
||||
|
||||
virtual void updatePropertiesTemp(const doublereal temp,
|
||||
doublereal* cp_R,
|
||||
doublereal* h_RT,
|
||||
doublereal* s_R) const = 0;
|
||||
|
||||
virtual void reportParameters(int &n, int &type,
|
||||
doublereal &tlow, doublereal &thigh,
|
||||
doublereal &pref,
|
||||
doublereal* const coeffs) const = 0;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -104,14 +104,19 @@ namespace Cantera {
|
|||
SpeciesThermoDuo() {}
|
||||
virtual ~SpeciesThermoDuo(){}
|
||||
|
||||
virtual void install(string name, int sp, int type, const doublereal* c,
|
||||
doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
|
||||
virtual void install(string name, int sp, int type,
|
||||
const doublereal* c,
|
||||
doublereal minTemp,
|
||||
doublereal maxTemp,
|
||||
doublereal refPressure) {
|
||||
m_p0 = refPressure;
|
||||
if (type == m_thermo1.ID) {
|
||||
m_thermo1.install(name, sp, 0, c, minTemp, maxTemp, refPressure);
|
||||
m_thermo1.install(name, sp, 0, c, minTemp, maxTemp,
|
||||
refPressure);
|
||||
speciesToType[sp] = m_thermo1.ID;
|
||||
} else if (type == m_thermo2.ID) {
|
||||
m_thermo2.install(name, sp, 0, c, minTemp, maxTemp, refPressure);
|
||||
m_thermo2.install(name, sp, 0, c, minTemp, maxTemp,
|
||||
refPressure);
|
||||
speciesToType[sp] = m_thermo2.ID;
|
||||
} else {
|
||||
throw UnknownSpeciesThermo("SpeciesThermoDuo:install",type);
|
||||
|
|
@ -136,7 +141,7 @@ namespace Cantera {
|
|||
return (tm1 < tm2 ? tm1 : tm2);
|
||||
}
|
||||
|
||||
virtual doublereal refPressure() const {
|
||||
virtual doublereal refPressure(int k = -1) const {
|
||||
return m_p0;
|
||||
}
|
||||
|
||||
|
|
@ -158,7 +163,7 @@ namespace Cantera {
|
|||
if (ctype == m_thermo1.ID) {
|
||||
m_thermo1.reportParams(index, type, c, minTemp, maxTemp,
|
||||
refPressure);
|
||||
} else if (ctype == m_thermo1.ID) {
|
||||
} else if (ctype == m_thermo2.ID) {
|
||||
m_thermo2.reportParams(index, type, c, minTemp, maxTemp,
|
||||
refPressure);
|
||||
} else {
|
||||
|
|
@ -225,7 +230,7 @@ namespace Cantera {
|
|||
return m_thermo[k].maxTemp();
|
||||
}
|
||||
|
||||
virtual doublereal refPressure() const {
|
||||
virtual doublereal refPressure(int k = -1) const {
|
||||
return m_pref;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -20,15 +20,23 @@
|
|||
|
||||
// NASA Polynomials
|
||||
#define NASA 4
|
||||
#define NASA2 4
|
||||
|
||||
// Shomate Polynomials used in NIST database
|
||||
#define SHOMATE 8
|
||||
#define SHOMATE2 8
|
||||
|
||||
// Tiger Polynomials
|
||||
#define TIGER 16
|
||||
|
||||
#define SIMPLE 32
|
||||
|
||||
#define MU0_INTERP 64
|
||||
|
||||
#define SHOMATE1 128
|
||||
|
||||
#define NASA1 256
|
||||
|
||||
#include "ct_defs.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue