From aceb896f624acab90d5899c0ca4447842053f4af Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 5 Mar 2019 14:34:40 -0500 Subject: [PATCH] [Thermo] Fix BinarySolutionTabulatedThermo updates when only T changes 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. --- .../thermo/BinarySolutionTabulatedThermo.h | 6 +++ src/thermo/BinarySolutionTabulatedThermo.cpp | 47 ++++++------------- .../BinarySolutionTabulatedThermo_Test.cpp | 9 +++- 3 files changed, 29 insertions(+), 33 deletions(-) diff --git a/include/cantera/thermo/BinarySolutionTabulatedThermo.h b/include/cantera/thermo/BinarySolutionTabulatedThermo.h index 11c2e9d7f..973164efa 100644 --- a/include/cantera/thermo/BinarySolutionTabulatedThermo.h +++ b/include/cantera/thermo/BinarySolutionTabulatedThermo.h @@ -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; diff --git a/src/thermo/BinarySolutionTabulatedThermo.cpp b/src/thermo/BinarySolutionTabulatedThermo.cpp index a90285795..166192a03 100644 --- a/src/thermo/BinarySolutionTabulatedThermo.cpp +++ b/src/thermo/BinarySolutionTabulatedThermo.cpp @@ -50,53 +50,36 @@ void BinarySolutionTabulatedThermo::compositionChanged() void BinarySolutionTabulatedThermo::_updateThermo() const { - double tnow = temperature(); double xnow = moleFraction(m_kk_tab); - std::pair 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; } } diff --git a/test/thermo/BinarySolutionTabulatedThermo_Test.cpp b/test/thermo/BinarySolutionTabulatedThermo_Test.cpp index fb9f41436..6145907f2 100644 --- a/test/thermo/BinarySolutionTabulatedThermo_Test.cpp +++ b/test/thermo/BinarySolutionTabulatedThermo_Test.cpp @@ -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); } }