Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix #8599 - Incorrect psychrometric calculation causes error #8833

Merged
merged 14 commits into from Sep 6, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
21 changes: 11 additions & 10 deletions src/EnergyPlus/DataGlobalConstants.hh
Expand Up @@ -164,16 +164,17 @@ namespace DataGlobalConstants {
std::string::size_type constexpr MaxNameLength =
100; // Maximum Name Length in Characters -- should be the same as MaxAlphaArgLength in InputProcessor module
Real64 constexpr KelvinConv = 273.15; // Conversion factor for C to K and K to C
Real64 constexpr InitConvTemp = 5.05; // [deg C], standard init vol to mass flow conversion temp
Real64 constexpr AutoCalculate = -99999.0; // automatically calculate some fields.
Real64 constexpr CWInitConvTemp = 5.05; // [deg C], standard init chilled water vol to mass flow conversion temp
Real64 constexpr HWInitConvTemp = 60.0; // [deg C], standard init hot water vol to mass flow conversion temp
Real64 constexpr SteamInitConvTemp = 100.0; // [deg C], standard init steam vol to mass flow conversion temp
Real64 constexpr StefanBoltzmann = 5.6697E-8; // Stefan-Boltzmann constant in W/(m2*K4)
Real64 constexpr UniversalGasConst = 8314.462175; // Universal Gas Constant (J/mol*K)
Real64 constexpr convertJtoGJ = 1.0E-9; // Conversion factor for J to GJ
int constexpr MaxSpeedLevels = 10; // Maximum number of speed that supports
int constexpr ScheduleAlwaysOn = -1; // Value when passed to schedule routines gives back 1.0 (on)
Real64 constexpr TriplePointOfWaterTempKelvin = 273.16; // The triple point of water, in Kelvin
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New constant

Real64 constexpr InitConvTemp = 5.05; // [deg C], standard init vol to mass flow conversion temp
Real64 constexpr AutoCalculate = -99999.0; // automatically calculate some fields.
Real64 constexpr CWInitConvTemp = 5.05; // [deg C], standard init chilled water vol to mass flow conversion temp
Real64 constexpr HWInitConvTemp = 60.0; // [deg C], standard init hot water vol to mass flow conversion temp
Real64 constexpr SteamInitConvTemp = 100.0; // [deg C], standard init steam vol to mass flow conversion temp
Real64 constexpr StefanBoltzmann = 5.6697E-8; // Stefan-Boltzmann constant in W/(m2*K4)
Real64 constexpr UniversalGasConst = 8314.462175; // Universal Gas Constant (J/mol*K)
Real64 constexpr convertJtoGJ = 1.0E-9; // Conversion factor for J to GJ
int constexpr MaxSpeedLevels = 10; // Maximum number of speed that supports
int constexpr ScheduleAlwaysOn = -1; // Value when passed to schedule routines gives back 1.0 (on)
int constexpr MaxCTFTerms = 19; // Maximum number of CTF terms allowed to still allow stability //Note Duplicate of DataHeatBalance::MaxCTFTerms

ResourceType AssignResourceTypeNum(std::string const &ResourceTypeChar);
Expand Down
4 changes: 3 additions & 1 deletion src/EnergyPlus/General.cc
Expand Up @@ -971,7 +971,9 @@ void Iterate(Real64 &ResultX, // ResultX is the final Iteration result passed b

// New guess calculated from LINEAR FIT of most recent two points
DY = Y0 - Y1;
if (std::abs(DY) < small) DY = small;
if (std::abs(DY) < small) {
DY = small;
}
Comment on lines +974 to +976
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unrelated, but it's handy to put a debugger breakpoint... (We should have clang-format do that everywhere IMHO)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. It's standard practice but not enforced.

// new estimation

ResultX = (Y0 * X1 - Y1 * X0) / DY;
Expand Down
60 changes: 36 additions & 24 deletions src/EnergyPlus/Psychrometrics.cc
Expand Up @@ -55,6 +55,7 @@
#include <EnergyPlus/CommandLineInterface.hh>
#include <EnergyPlus/Data/EnergyPlusData.hh>
#include <EnergyPlus/DataEnvironment.hh>
#include <EnergyPlus/DataGlobalConstants.hh>
#include <EnergyPlus/General.hh>
#include <EnergyPlus/Psychrometrics.hh>
#include <EnergyPlus/UtilityRoutines.hh>
Expand Down Expand Up @@ -498,7 +499,9 @@ namespace Psychrometrics {
for (iter = 1; iter <= itmax; ++iter) {

// Assigning a value to WBT
if (WBT >= (tBoil - 0.09)) WBT = tBoil - 0.1;
if (WBT >= (tBoil - 0.09)) {
WBT = tBoil - 0.1;
}

// Determine the saturation pressure for wet bulb temperature
PSatstar =
Expand All @@ -523,10 +526,9 @@ namespace Psychrometrics {
WBT = ResultX;

// If converged, leave iteration loop.
if (icvg == 1) break;

// Error Trap for the Discontinuous nature of PsyPsatFnTemp function (Sat Press Curve) at ~0 Deg C.
if ((PSatstar > 611.000) && (PSatstar < 611.25) && (std::abs(error) <= 0.0001) && (iter > 4)) break;
if (icvg == 1) {
break;
}
Comment on lines +529 to +531
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the "Error Trap for the Discontinuous nature of PsyPsatFnTemp function (Sat Press Curve) at ~0 Deg C." since it's no longer discontinuous.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I vaguely remember this. It looks like the unit tests hits this so good there. Gotta check it though.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I remember this now, this is a check for where this used to have a problem. Now that the problem is corrected, this check is no longer needed.


} // End of Iteration Loop

Expand Down Expand Up @@ -682,6 +684,11 @@ namespace Psychrometrics {
// Compared to Table 3 values (August 2007) with average error of 0.00%, max .30%,
// min -.39%. (Spreadsheet available on request - Lawrie).

// Note: the ASHRAE Handbook of Fundamentals is being slightly inaccurate in its wording,
// and referring to Eq 5 applying to -100°C to 0°C and Eq 6 applying to 0°C to 200°C
// In fact, it is **not** 0°C, but the triple-point of water, which is 0.01°C. Evaluating the Eq 5 and 6 up to and from the triple-point
// is removing the discontinuity altogether.
Comment on lines +687 to +690
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New note

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you were actually the one that found this critical point. Nice!


// USE STATEMENTS:

// Return value
Expand Down Expand Up @@ -734,29 +741,34 @@ namespace Psychrometrics {

// If below -100C,set value of Pressure corresponding to Saturation Temperature of -100C.
if (Tkel < 173.15) {
Pascal = 0.0017;
Pascal = 0.001405102123874164;
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rraustad and I noted that the value below -100°C does not match the evaluation of the equation at -100°C. So align.


// If below freezing, calculate saturation pressure over ice.
} else if (Tkel < DataGlobalConstants::KelvinConv) { // Tkel >= 173.15
Real64 constexpr C1(-5674.5359); // Coefficient for TKel < KelvinConvK
Real64 constexpr C2(6.3925247); // Coefficient for TKel < KelvinConvK
Real64 constexpr C3(-0.9677843e-2); // Coefficient for TKel < KelvinConvK
Real64 constexpr C4(0.62215701e-6); // Coefficient for TKel < KelvinConvK
Real64 constexpr C5(0.20747825e-8); // Coefficient for TKel < KelvinConvK
Real64 constexpr C6(-0.9484024e-12); // Coefficient for TKel < KelvinConvK
Real64 constexpr C7(4.1635019); // Coefficient for TKel < KelvinConvK
} else if (Tkel < DataGlobalConstants::TriplePointOfWaterTempKelvin) { // Tkel >= 173.15, Tkel < 273.16 (0.01°C)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The actual fix

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

image

jmarrec marked this conversation as resolved.
Show resolved Hide resolved
Real64 constexpr C1(-5674.5359);
Real64 constexpr C2(6.3925247);
Real64 constexpr C3(-0.9677843e-2);
Real64 constexpr C4(0.62215701e-6);
Real64 constexpr C5(0.20747825e-8);
Real64 constexpr C6(-0.9484024e-12);
Real64 constexpr C7(4.1635019);
Pascal = std::exp(C1 / Tkel + C2 + Tkel * (C3 + Tkel * (C4 + Tkel * (C5 + C6 * Tkel))) + C7 * std::log(Tkel));

// If above freezing, calculate saturation pressure over liquid water.
} else if (Tkel <= 473.15) { // Tkel >= 173.15 // Tkel >= KelvinConv
} else if (Tkel <= 473.15) { // Tkel >= 173.15 // Tkel >= TriplePointOfWaterTempKelvin
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The else if above is Tkel < Real64 constexpr TriplePointOfWaterTempKelvin = 273.16; and this is <= 473.15 ... oh yeah, this is the top end so comment is wrong.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These comments are a little confusing but a second looks shows they say what they said before with updates to what was changed. I'd say this is OK.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the comments provide little value overall here, but there were pre-existing. and they should be correct.

#ifndef EP_IF97
Real64 constexpr C8(-5800.2206); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C9(1.3914993); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C10(-0.048640239); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C11(0.41764768e-4); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C12(-0.14452093e-7); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C13(6.5459673); // Coefficient for TKel >= KelvinConvK
Real64 constexpr C8(-5800.2206);
Real64 constexpr C9(1.3914993);
Real64 constexpr C10(-0.048640239);
Real64 constexpr C11(0.41764768e-4);
Real64 constexpr C12(-0.14452093e-7);
Real64 constexpr C13(6.5459673);
Pascal = std::exp(C8 / Tkel + C9 + Tkel * (C10 + Tkel * (C11 + Tkel * C12)) + C13 * std::log(Tkel));

// If above 200C, set value of Pressure corresponding to Saturation Temperature of 200C.
} else { // Tkel >= 173.15 // Tkel >= TriplePointOfWaterTempKelvin // Tkel > 473.15
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So shouldn't this be the top end? 473.15 ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, it's only a comment issue.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same for this comment, these are each true and true as a whole.

Pascal = 1555073.745636215;
}
#else
// Table 34 in IF97
Real64 constexpr N1(0.11670521452767e04);
Expand All @@ -776,12 +788,12 @@ namespace Psychrometrics {
Real64 const B = N3 * phi2 + N4 * phi + N5;
Real64 const C = N6 * phi2 + N7 * phi + N8;
Pascal = 1000000.0 * pow_4((2.0 * C) / (-B + std::sqrt((B * B) - 4.0 * A * C)));
#endif

// If above 200C, set value of Pressure corresponding to Saturation Temperature of 200C.
} else { // Tkel >= 173.15 // Tkel >= KelvinConv // Tkel > 473.15
Pascal = 1555000.0;
Pascal = 1554671.8682698254;
}

#endif
return Pascal;
}

Expand Down
6 changes: 3 additions & 3 deletions tst/EnergyPlus/unit/AirTerminalSingleDuctMixer.unit.cc
Expand Up @@ -7816,11 +7816,11 @@ TEST_F(EnergyPlusFixture, AirTerminalSingleDuctMixer_SimFCU_ATMInletSideTest)
EXPECT_NEAR(thisFanCoil.PLR, 0.78843, 0.00001);
// check mass flow rates
EXPECT_NEAR(PrimaryAirMassFlowRate, 0.2, 0.000001);
EXPECT_NEAR(SecondaryAirMassFlowRate, 0.369710, 0.000001);
EXPECT_NEAR(SecondaryAirMassFlowRate, 0.369714, 0.000001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.AirInNode).MassFlowRate, thisFan.InletAirMassFlowRate, 0.000001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.ATMixerPriNode).MassFlowRate, 0.2, 0.0001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.ATMixerSecNode).MassFlowRate, 0.369710, 0.000001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.ATMixerOutNode).MassFlowRate, 0.569710, 0.000001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.ATMixerSecNode).MassFlowRate, 0.369714, 0.000001);
EXPECT_NEAR(state->dataLoopNodes->Node(thisFanCoil.ATMixerOutNode).MassFlowRate, 0.569714, 0.000001);
}

TEST_F(EnergyPlusFixture, AirTerminalSingleDuctMixer_FCU_NightCycleTest)
Expand Down
21 changes: 21 additions & 0 deletions tst/EnergyPlus/unit/Psychrometrics.unit.cc
Expand Up @@ -639,3 +639,24 @@ TEST_F(EnergyPlusFixture, Psychrometrics_CSpline_Test)
// check error
EXPECT_LE(error, 1E-5);
}

TEST_F(EnergyPlusFixture, Psychrometrics_PsyTwbFnTdbWPb_Test_Discontinuity)
{
// Test for #8599
InitializePsychRoutines(*state);

state->dataGlobal->WarmupFlag = true;

// Test when wet bulb temperature is approaching zero. PsyPsatFnTemp used to have a discontinuity in Psat around 0.0°C that makes the calculation
// blow up. Before: PsyPsatFnTemp(-0.0001) = 611.1485382610978 PsyPsatFnTemp(+0.0001) = 611.2173076397495
// diff = 0.06876937865172295
Real64 TDB = 1.4333333333333331; // C
Real64 W = 0.0031902374172088472; // Kg.water/Kg.dryair
Real64 Pb = 101400.00000000001;

Real64 result = PsyTwbFnTdbWPb(*state, TDB, W, Pb);
Real64 expected_result = -0.1027; // expected result from psychrometrics chart
EXPECT_NEAR(result, expected_result, 0.001);

EXPECT_FALSE(has_err_output());
}
30 changes: 15 additions & 15 deletions tst/EnergyPlus/unit/UnitarySystem.unit.cc
Expand Up @@ -9562,22 +9562,22 @@ TEST_F(EnergyPlusFixture, UnitarySystemModel_MultispeedDXHeatingCoilOnly)
thisSys->sizeSystem(*state, FirstHVACIteration, AirLoopNum);
DXCoils::SizeDXCoil(*state, 1);

ASSERT_EQ(1, state->dataUnitarySystems->numUnitarySystems); // only 1 unitary system above so expect 1 as number of unitary system objects
EXPECT_EQ(1, state->dataUnitarySystems->numUnitarySystems); // only 1 unitary system above so expect 1 as number of unitary system objects

ASSERT_NEAR(thisSys->m_DesignHeatingCapacity, 1303.091, 0.001);
ASSERT_EQ(thisSys->m_DesignCoolingCapacity, 0.0);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(1), 325.773, 0.001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(2), 651.545, 0.001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(3), 977.318, 0.001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(4), 1303.091, 0.001);
ASSERT_NEAR(thisSys->m_HeatVolumeFlowRate[1], 0.0131, 0.0001);
ASSERT_NEAR(thisSys->m_HeatVolumeFlowRate[2], 0.0262, 0.0001);
ASSERT_NEAR(thisSys->m_HeatVolumeFlowRate[3], 0.0393, 0.0001);
ASSERT_NEAR(thisSys->m_HeatVolumeFlowRate[4], 0.0524, 0.0001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(1), 0.0131, 0.0001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(2), 0.0262, 0.0001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(3), 0.0393, 0.0001);
ASSERT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(4), 0.0524, 0.0001);
EXPECT_NEAR(thisSys->m_DesignHeatingCapacity, 1303.097, 0.001);
EXPECT_EQ(thisSys->m_DesignCoolingCapacity, 0.0);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(1), 325.774, 0.001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(2), 651.549, 0.001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(3), 977.323, 0.001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedTotCap(4), 1303.097, 0.001);
EXPECT_NEAR(thisSys->m_HeatVolumeFlowRate[1], 0.0131, 0.0001);
EXPECT_NEAR(thisSys->m_HeatVolumeFlowRate[2], 0.0262, 0.0001);
EXPECT_NEAR(thisSys->m_HeatVolumeFlowRate[3], 0.0393, 0.0001);
EXPECT_NEAR(thisSys->m_HeatVolumeFlowRate[4], 0.0524, 0.0001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(1), 0.0131, 0.0001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(2), 0.0262, 0.0001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(3), 0.0393, 0.0001);
EXPECT_NEAR(state->dataDXCoils->DXCoil(1).MSRatedAirVolFlowRate(4), 0.0524, 0.0001);
}

TEST_F(EnergyPlusFixture, UnitarySystemModel_MultiSpeedCoils_SingleMode)
Expand Down