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

[Thermo] Fix BinarySolutionTabulatedThermo updates when only T changes

Browse source 
In the case where temperature changes but the mole fractions are the same, we
still need to apply the enthalpy and entropy offsets to the tabulated species.
This commit is contained in:
Ray Speth 2019-03-05 14:34:40 -05:00
parent 47005a5008
commit aceb896f62
3 changed files with 29 additions and 33 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
include/cantera/thermo/BinarySolutionTabulatedThermo.h
View file

@ -158,6 +158,12 @@ protected:
//! Current tabulated species mole fraction
mutable double m_xlast;
//! Tabulated contribution to h0[m_kk_tab] at the current composition
mutable double m_h0_tab;
//! Tabulated contribution to s0[m_kk_tab] at the current composition
mutable double m_s0_tab;
//! Vector for storing tabulated thermo
vector_fp m_molefrac_tab;
vector_fp m_enthalpy_tab;

47
src/thermo/BinarySolutionTabulatedThermo.cpp
View file

@ -50,53 +50,36 @@ void BinarySolutionTabulatedThermo::compositionChanged()
void BinarySolutionTabulatedThermo::_updateThermo() const
{
double tnow = temperature();
double xnow = moleFraction(m_kk_tab);
std::pair<double,double> d;
double dS_corr = 0.0;
bool x_changed = (m_xlast != xnow);
if (m_xlast != xnow) {
d = interpolate(xnow);
if (xnow == 0)
{
dS_corr = -BigNumber;
} else if (xnow == 1)
{
dS_corr = BigNumber;
} else
{
dS_corr = GasConstant*std::log(xnow/(1.0-xnow)) +
if (x_changed) {
std::tie(m_h0_tab, m_s0_tab) = interpolate(xnow);
if (xnow == 0) {
m_s0_tab = -BigNumber;
} else if (xnow == 1) {
m_s0_tab = BigNumber;
} else {
m_s0_tab += GasConstant*std::log(xnow/(1.0-xnow)) +
GasConstant/Faraday*std::log(standardConcentration(1-m_kk_tab)
/standardConcentration(m_kk_tab));
}
m_xlast = xnow;
}
double tnow = temperature();
if (x_changed || m_tlast != tnow) {
// Update the thermodynamic functions of the reference state.
m_spthermo.update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
m_tlast = tnow;
double rrt = 1.0 / RT();
double rr = 1.0 / GasConstant;
m_h0_RT[m_kk_tab] += m_h0_tab * rrt;
m_s0_R[m_kk_tab] += m_s0_tab / GasConstant;
for (size_t k = 0; k < m_kk; k++) {
double deltaE = rrt * m_pe[k];
m_h0_RT[k] += deltaE;
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
}
m_h0_RT[m_kk_tab] += d.first*rrt;
m_s0_R[m_kk_tab] += (d.second + dS_corr)*rr;
m_g0_RT[m_kk_tab] = m_h0_RT[m_kk_tab] - m_s0_R[m_kk_tab];
m_tlast = tnow;
m_xlast = xnow;
} else if (m_tlast != tnow) {
// Update the thermodynamic functions of the reference state.
m_spthermo.update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
m_tlast = tnow;
double rrt = 1.0 / RT();
for (size_t k = 0; k < m_kk; k++) {
double deltaE = rrt * m_pe[k];
m_h0_RT[k] += deltaE;
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
}
m_tlast = tnow;
}
}

9
test/thermo/BinarySolutionTabulatedThermo_Test.cpp
View file

@ -52,6 +52,10 @@ TEST_F(BinarySolutionTabulatedThermo_Test,interp_h)
{
set_defect_X(xmin + i*dx);
EXPECT_NEAR(expected_result[i], test_phase->enthalpy_mole(), 1.e-6);
// enthalpy is temperature-independent in test data file (all species
// use constant cp model with cp = 0)
test_phase->setState_TP(310, 101325);
EXPECT_NEAR(expected_result[i], test_phase->enthalpy_mole(), 1.e-6);
}
}
@ -78,7 +82,10 @@ TEST_F(BinarySolutionTabulatedThermo_Test,interp_s)
for (int i = 0; i < 9; ++i)
{
set_defect_X(xmin + i*dx);
EXPECT_NEAR(expected_result[i], test_phase->entropy_mole(), 1.e-6);
// entropy is temperature-independent in test data file (all species use
// constant cp model with cp = 0)
test_phase->setState_TP(330.0, 101325);
EXPECT_NEAR(expected_result[i], test_phase->entropy_mole(), 1.e-6);
}
}