research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Thermo] Adjust NASA polynomials to eliminate discontinuities

Browse source 
We expect the high- and low-temperature NASA polynomials for each species to be
continuous in a number of places, despite many frequently used thermo databases
having data that is inconsistent.

In addition to warning about discontinuities, We now modify the provided
coefficients to ensure continuity at the midpoint temperature. This resolves
numerical issues, e.g. with the algorithm used for setting the state of a
reactor network.
This commit is contained in:
Ray Speth 2013-04-18 22:07:00 +00:00
parent 548f0b1c56
commit ae96e48b4f
2 changed files with 36 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

17
src/thermo/NasaThermo.h
View file

@ -162,7 +162,7 @@ public:
vector_fp chigh(c+8, c+15);
vector_fp clow(c+1, c+8);
checkContinuity(name, tmid, &clow[0], &chigh[0]);
ensureContinuity(name, tmid, &clow[0], &chigh[0]);
m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh,
refPressure, &chigh[0]));
@ -460,17 +460,22 @@ protected:
mutable std::map<size_t, std::string> m_name;
private:
//! see SpeciesThermoFactory.cpp for the definition
//! Adjust polynomials to be continuous at the midpoint temperature.
/*!
* Check to see if the provided coefficients are nearly continuous. Adjust
* the values to get more precise contintinuity to avoid convergence
* issues with algorithms that expect these quantities to be continuous.
* See SpeciesThermoFactory.cpp for the definition.
*
* @param name string name of species
* @param tmid Mid temperature, between the two temperature regions
* @param clow coefficients for lower temperature region
* @param chigh coefficients for higher temperature region
*/
void checkContinuity(const std::string& name, double tmid,
const doublereal* clow, doublereal* chigh);
void ensureContinuity(const std::string& name, double tmid,
doublereal* clow, doublereal* chigh);
//! for internal use by checkContinuity
//! for internal use by ensureContinuity
/*!
* @param t temperature
* @param c coefficient array
@ -481,7 +486,7 @@ private:
+ c[0]/t;
}
//! for internal use by checkContinuity
//! for internal use by ensureContinuity
/*!
* @param t temperature
* @param c coefficient array

27
src/thermo/SpeciesThermoFactory.cpp
View file

@ -198,8 +198,8 @@ SpeciesThermo* SpeciesThermoFactory::newSpeciesThermoManager(std::string& stype)
return (SpeciesThermo*) 0;
}
void NasaThermo::checkContinuity(const std::string& name, double tmid,
const doublereal* clow, doublereal* chigh)
void NasaThermo::ensureContinuity(const std::string& name, double tmid,
doublereal* clow, doublereal* chigh)
{
// heat capacity
doublereal cplow = poly4(tmid, clow + 2);
@ -215,6 +215,13 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
+fp2str(cphigh)+".\n");
}
// Adjust coefficients to eliminate any discontinuity
chigh[2] += 0.5 * delta;
clow[2] -= 0.5 * delta;
AssertThrowMsg(std::abs(poly4(tmid, clow+2) - poly4(tmid, chigh+2)) < 1e-12,
"NasaThermo::ensureContinuity", "Cp/R does not match");
// enthalpy
doublereal hrtlow = enthalpy_RT(tmid, clow);
doublereal hrthigh = enthalpy_RT(tmid, chigh);
@ -229,6 +236,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
+fp2str(hrthigh)+".\n");
}
// Adjust coefficients to eliminate any discontinuity
chigh[0] += 0.5 * delta * tmid;
clow[0] -= 0.5 * delta * tmid;
AssertThrowMsg(std::abs(enthalpy_RT(tmid, clow) -
enthalpy_RT(tmid, chigh)) < 1e-12,
"NasaThermo::ensureContinuity", "H/RT does not match");
// entropy
doublereal srlow = entropy_R(tmid, clow);
doublereal srhigh = entropy_R(tmid, chigh);
@ -242,6 +257,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(srhigh)+".\n");
}
// Adjust coefficients to eliminate any discontinuity
chigh[1] += 0.5 * delta;
clow[1] -= 0.5 * delta;
AssertThrowMsg(std::abs(entropy_R(tmid, clow) -
entropy_R(tmid, chigh)) < 1e-12,
"NasaThermo::ensureContinuity", "S/R does not match");
}