465 lines
15 KiB
C++
465 lines
15 KiB
C++
/*!
|
|
* @file NasaThermo.cpp Implementation of class Cantera::NasaThermo
|
|
*/
|
|
#include "NasaThermo.h"
|
|
|
|
#include "cantera/base/utilities.h"
|
|
#include "cantera/numerics/DenseMatrix.h"
|
|
#include "cantera/numerics/ctlapack.h"
|
|
|
|
namespace Cantera
|
|
{
|
|
|
|
NasaThermo::NasaThermo() :
|
|
ID(NASA),
|
|
m_tlow_max(0.0),
|
|
m_thigh_min(1.e30),
|
|
m_p0(-1.0),
|
|
m_ngroups(0) {
|
|
m_t.resize(6);
|
|
}
|
|
|
|
NasaThermo::NasaThermo(const NasaThermo& right) :
|
|
ID(NASA),
|
|
m_tlow_max(0.0),
|
|
m_thigh_min(1.e30),
|
|
m_p0(-1.0),
|
|
m_ngroups(0) {
|
|
*this = operator=(right);
|
|
}
|
|
|
|
NasaThermo& NasaThermo::operator=(const NasaThermo& right)
|
|
{
|
|
/*
|
|
* Check for self assignment.
|
|
*/
|
|
if (this == &right) {
|
|
return *this;
|
|
}
|
|
|
|
m_high = right.m_high;
|
|
m_low = right.m_low;
|
|
m_index = right.m_index;
|
|
m_tmid = right.m_tmid;
|
|
m_tlow_max = right.m_tlow_max;
|
|
m_thigh_min = right.m_thigh_min;
|
|
m_tlow = right.m_tlow;
|
|
m_thigh = right.m_thigh;
|
|
m_p0 = right.m_p0;
|
|
m_ngroups = right.m_ngroups;
|
|
m_t = right.m_t;
|
|
m_group_map = right.m_group_map;
|
|
m_posInGroup_map = right.m_posInGroup_map;
|
|
m_name = right.m_name;
|
|
|
|
return *this;
|
|
}
|
|
|
|
void NasaThermo::install(const std::string& name, size_t index, int type,
|
|
const doublereal* c,
|
|
doublereal min_temp, doublereal max_temp,
|
|
doublereal ref_pressure)
|
|
{
|
|
m_name[index] = name;
|
|
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
|
|
std::vector<NasaPoly1> v;
|
|
m_high.push_back(v);
|
|
m_low.push_back(v);
|
|
m_tmid.push_back(c[0]);
|
|
m_index[imid] = igrp = static_cast<int>(m_high.size());
|
|
m_ngroups++;
|
|
}
|
|
|
|
m_group_map[index] = igrp;
|
|
m_posInGroup_map[index] = (int) m_low[igrp-1].size();
|
|
|
|
doublereal tlow = min_temp;
|
|
doublereal tmid = c[0];
|
|
doublereal thigh = max_temp;
|
|
|
|
vector_fp chigh(c+8, c+15);
|
|
vector_fp clow(c+1, c+8);
|
|
|
|
checkContinuity(name, tmid, &clow[0], &chigh[0]);
|
|
|
|
|
|
m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh,
|
|
ref_pressure, &chigh[0]));
|
|
m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid,
|
|
ref_pressure, &clow[0]));
|
|
|
|
m_tlow_max = std::max(tlow, m_tlow_max);
|
|
m_thigh_min = std::min(thigh, m_thigh_min);
|
|
if (m_tlow.size() < index + 1) {
|
|
m_tlow.resize(index + 1, tlow);
|
|
m_thigh.resize(index + 1, thigh);
|
|
}
|
|
m_tlow[index] = tlow;
|
|
m_thigh[index] = thigh;
|
|
if (m_p0 < 0.0) {
|
|
m_p0 = ref_pressure;
|
|
} else if (fabs(m_p0 - ref_pressure) > 0.1) {
|
|
std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, "
|
|
+ fp2str(ref_pressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
|
|
writelog(logmsg);
|
|
logmsg = " This is now a fatal error\n";
|
|
writelog(logmsg);
|
|
throw CanteraError("install()", "species have different reference pressures");
|
|
}
|
|
m_p0 = ref_pressure;
|
|
}
|
|
|
|
void NasaThermo::update_one(size_t k, doublereal t, doublereal* cp_R,
|
|
doublereal* h_RT, doublereal* s_R) const
|
|
{
|
|
m_t[0] = t;
|
|
m_t[1] = t*t;
|
|
m_t[2] = m_t[1]*t;
|
|
m_t[3] = m_t[2]*t;
|
|
m_t[4] = 1.0/t;
|
|
m_t[5] = log(t);
|
|
|
|
size_t grp = m_group_map[k];
|
|
size_t pos = m_posInGroup_map[k];
|
|
const std::vector<NasaPoly1> &mlg = m_low[grp-1];
|
|
const NasaPoly1* nlow = &(mlg[pos]);
|
|
|
|
doublereal tmid = nlow->maxTemp();
|
|
if (t < tmid) {
|
|
nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R);
|
|
} else {
|
|
const std::vector<NasaPoly1> &mhg = m_high[grp-1];
|
|
const NasaPoly1* nhigh = &(mhg[pos]);
|
|
nhigh->updateProperties(&m_t[0], cp_R, h_RT, s_R);
|
|
}
|
|
}
|
|
|
|
void NasaThermo::update(doublereal t, doublereal* cp_R,
|
|
doublereal* h_RT, doublereal* s_R) const
|
|
{
|
|
int i;
|
|
|
|
// load functions of temperature into m_t vector
|
|
m_t[0] = t;
|
|
m_t[1] = t*t;
|
|
m_t[2] = m_t[1]*t;
|
|
m_t[3] = m_t[2]*t;
|
|
m_t[4] = 1.0/t;
|
|
m_t[5] = log(t);
|
|
|
|
// iterate over the groups
|
|
std::vector<NasaPoly1>::const_iterator _begin, _end;
|
|
for (i = 0; i != m_ngroups; i++) {
|
|
if (t > m_tmid[i]) {
|
|
_begin = m_high[i].begin();
|
|
_end = m_high[i].end();
|
|
} else {
|
|
_begin = m_low[i].begin();
|
|
_end = m_low[i].end();
|
|
}
|
|
for (; _begin != _end; ++_begin) {
|
|
_begin->updateProperties(&m_t[0], cp_R, h_RT, s_R);
|
|
}
|
|
}
|
|
}
|
|
|
|
void NasaThermo::reportParams(size_t index, int& type,
|
|
doublereal* const c,
|
|
doublereal& minTemp,
|
|
doublereal& maxTemp,
|
|
doublereal& refPressure) const
|
|
{
|
|
type = reportType(index);
|
|
if (type == NASA) {
|
|
size_t grp = m_group_map[index];
|
|
size_t pos = m_posInGroup_map[index];
|
|
const std::vector<NasaPoly1> &mlg = m_low[grp-1];
|
|
const std::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;
|
|
size_t n;
|
|
double ttemp;
|
|
lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
|
|
c + 1);
|
|
if (n != index) {
|
|
throw CanteraError(" ", "confused");
|
|
}
|
|
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");
|
|
}
|
|
}
|
|
|
|
doublereal NasaThermo::reportOneHf298(const size_t k) const
|
|
{
|
|
size_t grp = m_group_map[k];
|
|
size_t pos = m_posInGroup_map[k];
|
|
const std::vector<NasaPoly1> &mlg = m_low[grp-1];
|
|
const NasaPoly1* nlow = &(mlg[pos]);
|
|
doublereal tmid = nlow->maxTemp();
|
|
double h;
|
|
if (298.15 <= tmid) {
|
|
h = nlow->reportHf298(0);
|
|
} else {
|
|
const std::vector<NasaPoly1> &mhg = m_high[grp-1];
|
|
const NasaPoly1* nhigh = &(mhg[pos]);
|
|
h = nhigh->reportHf298(0);
|
|
}
|
|
return h;
|
|
}
|
|
|
|
void NasaThermo::modifyOneHf298(const size_t k, const doublereal Hf298New)
|
|
{
|
|
size_t grp = m_group_map[k];
|
|
size_t pos = m_posInGroup_map[k];
|
|
std::vector<NasaPoly1> &mlg = m_low[grp-1];
|
|
NasaPoly1* nlow = &(mlg[pos]);
|
|
std::vector<NasaPoly1> &mhg = m_high[grp-1];
|
|
NasaPoly1* nhigh = &(mhg[pos]);
|
|
doublereal tmid = nlow->maxTemp();
|
|
|
|
double hnow = reportOneHf298(k);
|
|
double delH = Hf298New - hnow;
|
|
if (298.15 <= tmid) {
|
|
nlow->modifyOneHf298(k, Hf298New);
|
|
double h = nhigh->reportHf298(0);
|
|
double hnew = h + delH;
|
|
nhigh->modifyOneHf298(k, hnew);
|
|
} else {
|
|
nhigh->modifyOneHf298(k, Hf298New);
|
|
double h = nlow->reportHf298(0);
|
|
double hnew = h + delH;
|
|
nlow->modifyOneHf298(k, hnew);
|
|
}
|
|
}
|
|
|
|
doublereal NasaThermo::cp_R(double t, const doublereal* c)
|
|
{
|
|
return poly4(t, c+2);
|
|
}
|
|
|
|
doublereal NasaThermo::enthalpy_RT(double t, const doublereal* c) {
|
|
return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t
|
|
+ 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t
|
|
+ c[0]/t;
|
|
}
|
|
|
|
doublereal NasaThermo::entropy_R(double t, const doublereal* c) {
|
|
return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t
|
|
+ OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t
|
|
+ c[1];
|
|
}
|
|
|
|
doublereal NasaThermo::checkContinuity(const std::string& name, double tmid,
|
|
doublereal* clow, doublereal* chigh)
|
|
{
|
|
// heat capacity
|
|
doublereal cplow = cp_R(tmid, clow);
|
|
doublereal cphigh = cp_R(tmid, chigh);
|
|
doublereal delta = cplow - cphigh;
|
|
doublereal maxError = std::abs(delta);
|
|
if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) {
|
|
writelog("\n\n**** WARNING ****\nFor species "+name+
|
|
", discontinuity in cp/R detected at Tmid = "
|
|
+fp2str(tmid)+"\n");
|
|
writelog("\tValue computed using low-temperature polynomial: "
|
|
+fp2str(cplow)+".\n");
|
|
writelog("\tValue computed using high-temperature polynomial: "
|
|
+fp2str(cphigh)+".\n");
|
|
}
|
|
|
|
// enthalpy
|
|
doublereal hrtlow = enthalpy_RT(tmid, clow);
|
|
doublereal hrthigh = enthalpy_RT(tmid, chigh);
|
|
delta = hrtlow - hrthigh;
|
|
maxError = std::max(std::abs(delta), maxError);
|
|
if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) {
|
|
writelog("\n\n**** WARNING ****\nFor species "+name+
|
|
", discontinuity in h/RT detected at Tmid = "
|
|
+fp2str(tmid)+"\n");
|
|
writelog("\tValue computed using low-temperature polynomial: "
|
|
+fp2str(hrtlow)+".\n");
|
|
writelog("\tValue computed using high-temperature polynomial: "
|
|
+fp2str(hrthigh)+".\n");
|
|
}
|
|
|
|
// entropy
|
|
doublereal srlow = entropy_R(tmid, clow);
|
|
doublereal srhigh = entropy_R(tmid, chigh);
|
|
delta = srlow - srhigh;
|
|
maxError = std::max(std::abs(delta), maxError);
|
|
if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) {
|
|
writelog("\n\n**** WARNING ****\nFor species "+name+
|
|
", discontinuity in s/R detected at Tmid = "
|
|
+fp2str(tmid)+"\n");
|
|
writelog("\tValue computed using low-temperature polynomial: "
|
|
+fp2str(srlow)+".\n");
|
|
writelog("\tValue computed using high-temperature polynomial: "
|
|
+fp2str(srhigh)+".\n");
|
|
}
|
|
|
|
return maxError;
|
|
}
|
|
|
|
void NasaThermo::fixDiscontinuities(doublereal Tlow, doublereal Tmid,
|
|
doublereal Thigh, doublereal* clow,
|
|
doublereal* chigh)
|
|
{
|
|
// The thermodynamic parameters can be written in terms nondimensionalized
|
|
// coefficients A[i] and the nondimensional temperature t = T/Tmid as:
|
|
//
|
|
// C_low(t) = A[0] + A[i] * t**i
|
|
// H_low(t) = A[0] + A[i] / (i+1) * t**i + A[5] / t
|
|
// S_low(t) = A[0]*ln(t) + A[i] / i * t**i + A[6]
|
|
//
|
|
// where the implicit sum is over the range 1 <= i <= 4 and the
|
|
// nondimensional coefficients are related to the dimensional coefficients
|
|
// a[i] by:
|
|
//
|
|
// A[0] = a[0]
|
|
// A[i] = Tmid**i * a[i], 1 <= i <= 4
|
|
// A[5] = a[5] / Tmid
|
|
// A[6] = a[6] + a[0] * ln(Tmid)
|
|
//
|
|
// and corresponding relationships hold for the high-temperature
|
|
// polynomial coefficients B[i]. This nondimensionalization is necessary
|
|
// in order for the resulting matrix to be well-conditioned.
|
|
//
|
|
// The requirement that C_low(1) = C_high(1) is satisfied by:
|
|
//
|
|
// B[0] = A[0] + (A[i] - B[i])
|
|
// C_high(t) = A[0] + (A[i] + B[i] * t**i - 1)
|
|
//
|
|
// The requirement that H_low(1) = H_high(1) is satisfied by:
|
|
//
|
|
// B[5] = A[5] + (i / (i+1) * (B[i] - A[i]))
|
|
// H_high(t) = A[0] + A[5] / t + (1 - i / (i+1) / t) * A[i] +
|
|
// (t**i / (i+1) - 1 + i / (i+1) / t) * B[i]
|
|
//
|
|
// The requirement that S_low(1) = S_high(1) is satisfied by:
|
|
//
|
|
// B[6] = A[6] + (A[i] - B[i]) / i
|
|
// S_high(t) = A[0] * ln(t) + A[6] + (ln(t) + 1 / i) * A[i] +
|
|
// (-ln(t) + t**i / i - 1 / i) * B[i]
|
|
|
|
// Formulate a linear least squares problem for the nondimensionalized
|
|
// coefficients. In the system of equations M*x = b:
|
|
// - each row of M consists of the factors in one of the above equations
|
|
// for C_low, H_high, etc. evaluated at some temperature between Tlow
|
|
// and Thigh
|
|
// - x is a vector of the 11 independent coefficients (A[0] through A[6]
|
|
// and B[1] through B[4])
|
|
// - B is a vector of the corresponding value of C, H, or S computed using
|
|
// the original polynomial.
|
|
|
|
const size_t nTemps = 12;
|
|
const size_t nCols = 11; // number of independent coefficients
|
|
const size_t nRows = 3*nTemps; // Evaluate C, H, and S at each temperature
|
|
DenseMatrix M(nRows, nCols, 0.0);
|
|
vector_fp b(nRows);
|
|
doublereal sqrtDeltaT = sqrt(Thigh) - sqrt(Tlow);
|
|
vector_fp tpow(5);
|
|
for (size_t j = 0; j < nTemps; j++) {
|
|
double T = pow(sqrt(Tlow) + sqrtDeltaT * j / (nTemps - 1.0), 2);
|
|
double t = T / Tmid; // non-dimensionalized temperature
|
|
double logt = std::log(t);
|
|
size_t n = 3 * j; // row index
|
|
for (int i = 1; i <= 4; i++) {
|
|
tpow[i] = pow(t, i);
|
|
}
|
|
|
|
// row n: Cp/R
|
|
// row n+1: H/RT
|
|
// row n+2: S/R
|
|
// columns 0 through 6 are for the low-T coefficients
|
|
// columns 7 through 10 are for the independent high-T coefficients
|
|
M(n, 0) = 1.0;
|
|
M(n+1,0) = 1.0;
|
|
M(n+2,0) = logt;
|
|
M(n+1,5) = 1.0 / t;
|
|
M(n+2,6) = 1.0;
|
|
if (t <= 1.0) {
|
|
for (int i = 1; i <= 4; i++) {
|
|
M(n,i) = tpow[i];
|
|
M(n+1,i) = tpow[i] / (i+1);
|
|
M(n+2,i) = tpow[i] / i;
|
|
}
|
|
b[n] = cp_R(T, clow);
|
|
b[n+1] = enthalpy_RT(T, clow);
|
|
b[n+2] = entropy_R(T, clow);
|
|
} else {
|
|
for (int i = 1; i <= 4; i++) {
|
|
M(n,i) = 1.0;
|
|
M(n,i+6) = tpow[i] - 1.0;
|
|
M(n+1,i) = 1 - i / ((i + 1.0) * t);
|
|
M(n+1,i+6) = -1 + tpow[i] / (i+1) + i / ((i+1) * t);
|
|
M(n+2,i) = logt + 1.0 / i;
|
|
M(n+2,i+6) = -logt + (tpow[i] - 1.0) / i;
|
|
}
|
|
b[n] = cp_R(T, chigh);
|
|
b[n+1] = enthalpy_RT(T, chigh);
|
|
b[n+2] = entropy_R(T, chigh);
|
|
}
|
|
}
|
|
|
|
// Solve the least squares problem
|
|
vector_fp sigma(nRows);
|
|
size_t rank;
|
|
int info;
|
|
vector_fp work(1);
|
|
int lwork = -1;
|
|
// First get the desired size of the work array
|
|
ct_dgelss(nRows, nCols, 1, &M(0,0), nRows, &b[0], nRows,
|
|
&sigma[0], -1, rank, &work[0], lwork, info);
|
|
work.resize(static_cast<size_t>(work[0]));
|
|
lwork = static_cast<int>(work[0]);
|
|
ct_dgelss(nRows, nCols, 1, &M(0,0), nRows, &b[0], nRows,
|
|
&sigma[0], -1, rank, &work[0], lwork, info);
|
|
|
|
AssertTrace(info == 0);
|
|
AssertTrace(rank == nCols);
|
|
AssertTrace(sigma[0] / sigma[10] < 1e20); // condition number
|
|
|
|
// Compute the full set of nondimensionalized coefficients
|
|
// (dgelss returns the solution of M*x = b in b).
|
|
|
|
// Note that clow and chigh store the coefficients in the order:
|
|
// clow = [a[5], a[6], a[0], a[1], a[2], a[3], a[4]]
|
|
clow[2] = chigh[2] = b[0];
|
|
clow[0] = chigh[0] = b[5];
|
|
clow[1] = chigh[1] = b[6];
|
|
for (int i = 1; i <= 4; i++) {
|
|
clow[2+i] = b[i];
|
|
chigh[2+i] = b[6+i];
|
|
chigh[2] += clow[2+i] - chigh[2+i];
|
|
chigh[0] += i / (i + 1.0) * (chigh[2+i] - clow[2+i]);
|
|
chigh[1] += (clow[2+i] - chigh[2+i]) / i;
|
|
}
|
|
|
|
// redimensionalize
|
|
for (int i = 1; i <= 4; i++) {
|
|
clow[2+i] /= pow(Tmid, i);
|
|
chigh[2+i] /= pow(Tmid, i);
|
|
}
|
|
clow[0] *= Tmid;
|
|
chigh[0] *= Tmid;
|
|
clow[1] -= clow[2] * std::log(Tmid);
|
|
chigh[1] -= chigh[2] * std::log(Tmid);
|
|
}
|
|
|
|
}
|