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

added support for modifying NASA polynomial coefficients

Browse source
This commit is contained in:
Dave Goodwin 2006-11-06 15:20:34 +00:00
parent d1029224fb
commit 704e084347
6 changed files with 85 additions and 42 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
Cantera/src/Constituents.h
View file

@ -194,27 +194,6 @@ namespace Cantera {
*/
doublereal size(int k) const { return m_speciesSize[k]; }
#ifdef INCL_DEPRECATED_METHODS
/**
* Return a SpeciesData structure containing species data.
*/
SpeciesData species(int k) const {
if (k < 0 || k >= nSpecies())
throw SpeciesRangeError("Constituents::charge",k,nSpecies());
SpeciesData s;
s.name = m_speciesNames[k];
s.phase = m_speciesPhase[k];
int offset = m_mm * k;
s.atoms.resize(m_mm);
for (int m = 0; m < m_mm; m++) {
s.atoms[m] = m_speciesComp[offset + m];
}
s.charge = m_speciesCharge[k];
s.molecularWeight = m_weight[k];
return s;
}
#endif
/**
* Prohibit addition of more species, and prepare for
* calculations with this set of elements and species.

33
Cantera/src/NasaPoly1.h
View file

@ -94,9 +94,9 @@ namespace Cantera {
return (SpeciesThermoInterpType *) np;
}
doublereal minTemp() const { return m_lowT;}
doublereal maxTemp() const { return m_highT;}
doublereal refPressure() const { return m_Pref; }
virtual doublereal minTemp() const { return m_lowT;}
virtual doublereal maxTemp() const { return m_highT;}
virtual doublereal refPressure() const { return m_Pref; }
virtual int reportType() const { return NASA1; }
/**
@ -116,7 +116,7 @@ namespace Cantera {
* tt[4] = 1.0/t;
* tt[5] = log(t);
*/
void updateProperties(const doublereal* tt,
virtual void updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
doublereal ct0 = m_coeff[2]; // a0
@ -148,7 +148,7 @@ namespace Cantera {
*
* (note: this is slow, but it is general)
*/
void updatePropertiesTemp(const doublereal temp,
virtual void updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
double tPoly[6];
@ -161,7 +161,7 @@ namespace Cantera {
updateProperties(tPoly, cp_R, h_RT, s_R);
}
void reportParameters(int &n, int &type,
virtual void reportParameters(int &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const {
@ -170,11 +170,28 @@ namespace Cantera {
tlow = m_lowT;
thigh = m_highT;
pref = m_Pref;
for (int i = 0; i < 7; i++) {
coeffs[i] = m_coeff[i];
coeffs[5] = m_coeff[0];
coeffs[6] = m_coeff[1];
for (int i = 2; i < 7; i++) {
coeffs[i-2] = m_coeff[i];
}
#ifdef WARN_ABOUT_CHANGES_FROM_VERSION_1_6
cout << "************************************************\n"
cout << "Warning: NasaPoly1::reportParameters now returns \n"
<< "the coefficient array in the same order as in\n"
<< "the input file. See file NasaPoly1.h" << endl;
cout << "************************************************\n"
#endif
}
virtual void modifyParameters(doublereal* coeffs) {
m_coeff[0] = coeffs[5];
m_coeff[1] = coeffs[6];
for (int i = 0; i < 5; i++) {
m_coeff[i+2] = coeffs[i];
}
}
protected:
doublereal m_lowT; // lowest valid temperature

57
Cantera/src/NasaThermo.h
View file

@ -74,6 +74,7 @@ namespace Cantera {
doublereal minTemp, doublereal maxTemp,
doublereal refPressure) {
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
@ -97,13 +98,21 @@ namespace Cantera {
vector_fp chigh(7);
copy(c + 8, c + 15, chigh.begin());
checkContinuity(name, tmid, clow, &chigh[0]);
m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh,
pref, &chigh[0]));
m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid,
pref, clow));
vector_fp clu(7), chu(7);
clu[5] = clow[0];
clu[6] = clow[1];
copy(clow+2, clow+7, clu.begin());
chu[5] = chigh[0];
chu[6] = chigh[1];
copy(chigh.begin()+2, chigh.begin()+7, chu.begin());
checkContinuity(name, tmid, &clu[0], &chu[0]);
if (tlow > m_tlow_max) m_tlow_max = tlow;
if (thigh < m_thigh_min) m_thigh_min = thigh;
m_tlow.push_back(tlow);
@ -244,6 +253,37 @@ namespace Cantera {
}
}
/**
* This utility function modifies the array of coefficients.
* The array is the same as that returned by reportParams, so
* a call can first be made to reportParams to populate the
* array, and then modifyParams can be called to alter
* selected values. For the NASA object, there are 15
* coefficients.
*/
virtual void modifyParams(int index, doublereal *c) {
int type = reportType(index);
if (type == NASA) {
int grp = m_group_map[index];
int pos = m_posInGroup_map[index];
vector<NasaPoly1> &mlg = m_low[grp-1];
vector<NasaPoly1> &mhg = m_high[grp-1];
NasaPoly1 *lowPoly = &(mlg[pos]);
NasaPoly1 *highPoly = &(mhg[pos]);
doublereal tmid = lowPoly->maxTemp();
if (c[0] != tmid) {
throw CanteraError(" ", "Tmid cannot be changed");
}
lowPoly->modifyParameters(c + 1);
highPoly->modifyParameters(c + 8);
checkContinuity(m_name[index], c[0], c + 1, c + 8);
} else {
throw CanteraError(" ", "confused");
}
}
protected:
vector<vector<NasaPoly1> > m_high;
@ -271,6 +311,7 @@ namespace Cantera {
* temperature polynomials for that species are storred.
*/
mutable map<int, int> m_posInGroup_map;
mutable map<int, string> m_name;
private:
@ -280,16 +321,16 @@ namespace Cantera {
/// for internal use by checkContinuity
doublereal 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;
return c[0] + 0.5*c[1]*t + OneThird*c[2]*t*t
+ 0.25*c[3]*t*t*t + 0.2*c[4]*t*t*t*t
+ c[5]/t;
}
/// for internal use by checkContinuity
doublereal 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];
return c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t
+ OneThird*c[3]*t*t*t + 0.25*c[4]*t*t*t*t
+ c[6];
}
};

3
Cantera/src/SpeciesThermo.h
View file

@ -159,6 +159,9 @@ namespace Cantera {
doublereal &minTemp,
doublereal &maxTemp,
doublereal &refPressure)=0;
virtual void modifyParams(int index, doublereal *c) {}
};
}

10
Cantera/src/SpeciesThermoFactory.cpp
View file

@ -160,11 +160,11 @@ namespace Cantera {
doublereal* chigh) {
// heat capacity
doublereal cplow = poly4(tmid, clow+2);
doublereal cphigh = poly4(tmid, chigh+2);
doublereal cplow = poly4(tmid, clow);
doublereal cphigh = poly4(tmid, chigh);
doublereal delta = cplow - cphigh;
if (fabs(delta/cplow) > 0.001) {
writelog("\n**** WARNING ****\nFor species "+name+
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in cp/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
@ -178,7 +178,7 @@ namespace Cantera {
doublereal hrthigh = enthalpy_RT(tmid, chigh);
delta = hrtlow - hrthigh;
if (fabs(delta/hrtlow) > 0.001) {
writelog("\n**** WARNING ****\nFor species "+name+
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in h/RT detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
@ -192,7 +192,7 @@ namespace Cantera {
doublereal srhigh = entropy_R(tmid, chigh);
delta = srlow - srhigh;
if (fabs(delta/srlow) > 0.001) {
writelog("\n**** WARNING ****\nFor species "+name+
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in s/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "

3
Cantera/src/SpeciesThermoInterpType.h
View file

@ -47,6 +47,9 @@ namespace Cantera {
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const = 0;
virtual void modifyParameters(doublereal* coeffs) {}
};
}