From c2c1d28b7f8ba22a8f78f26e8bbd0cc404283816 Mon Sep 17 00:00:00 2001
From: Matthew Boyd <30417543+Matthew-Boyd@users.noreply.github.com>
Date: Thu, 20 Apr 2023 07:11:15 -0600
Subject: [PATCH] Fix cold tank issues in the solar water heating model (#1011)

* Fix unphysical code node tank temperature

* Deliver timestep averaged temperature when tank is drained

* Fix hot tank node temperature after draining

* Update swh test
---
 ssc/cmod_swh.cpp                | 27 ++++++++++++++++-----------
 test/ssc_test/cmod_swh_test.cpp |  2 +-
 2 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/ssc/cmod_swh.cpp b/ssc/cmod_swh.cpp
index 81a25783f..0309fb752 100644
--- a/ssc/cmod_swh.cpp
+++ b/ssc/cmod_swh.cpp
@@ -672,7 +672,8 @@ class cm_swh : public compute_module
 						V_hot = V_hot_prev + ts_sec*mdot_total/rho_water;
 						V_cold = V_tank - V_hot;
 						T_hot = (T_hot_prev*V_hot_prev + ts_sec*(mdot_total/rho_water)*(T_cold_prev + dT_collector))/V_hot;
-						T_cold = (V_tank/V_cold)*T_tank - (V_hot/V_cold)*T_hot;
+                        T_cold = (V_tank/V_cold)*T_tank - (V_hot/V_cold)*T_hot;                                 // weighted average to enforce T_tank based on T_hot
+                        if (T_cold < std::min(T_mains_use, T_room)) T_cold = std::min(T_mains_use, T_room);     // above relation breaks-down at small V_cold causing unphysical T_cold
 						T_top = T_hot;
 						T_bot = T_cold;
 						T_deliv = T_top;
@@ -715,9 +716,11 @@ class cm_swh : public compute_module
 						V_hot = 0;
                     }
 
+                    double T_hot_drained = std::numeric_limits<double>::quiet_NaN();
 					if (V_hot == 0)	// cold water drawn into the bottom of the tank in previous timesteps has completely flushed hot water from the tank
 					{
-						T_hot = T_hot_prev;
+                        double time_to_drain_sec = V_hot_prev * rho_water / mdot_mix;
+                        T_hot_drained = (T_hot_prev * time_to_drain_sec + T_cold * (ts_sec - time_to_drain_sec)) / ts_sec;
 					}
 					else
 					{
@@ -726,7 +729,6 @@ class cm_swh : public compute_module
 						double m_hot = V_hot_prev*rho_water;
 						T_hot = ((T_hot_prev * Cp_water * m_hot) + (ts_sec*U_tank*A_hot * T_room))/((m_hot*Cp_water) + (ts_sec*U_tank*A_hot)); // IMPLICIT NON-STEADY (Euler)
 					}
-					hotLoss = U_tank * A_hot * (T_hot - T_room);
 
 					// Cold node calculations
 					V_cold = V_tank-V_hot;
@@ -742,20 +744,23 @@ class cm_swh : public compute_module
 						T_cold = ((T_cold_prev*m_cold*Cp_water) + (ts_sec*U_tank*A_cold*T_room) + (ts_sec*mdot_mix*Cp_water*T_mains_use))
 							/((m_cold*Cp_water) + (ts_sec*A_cold*U_tank) + (mdot_mix*ts_sec*Cp_water) ); // IMPLICIT NON-STEADY
 					}
-					coldLoss = U_tank*A_cold*(T_cold - T_room);
 
-					Q_tankloss = hotLoss + coldLoss;
-					T_tank = (V_hot / V_tank) * T_hot + (V_cold / V_tank) * T_cold;
-					T_top = T_tank + 0.33*dT_collector;
-					T_bot = T_tank - 0.67*dT_collector;
-					// T_top = T_hot
-					// T_bot = T_cold
                     if (V_hot > 0) {
                         T_deliv = T_hot;
                     }
                     else {
-                        T_deliv = T_cold;
+                        T_deliv = T_hot_drained;
+                        T_hot = T_cold;     // hot water completely flushed from tank
                     }
+					T_tank = (V_hot / V_tank) * T_hot + (V_cold / V_tank) * T_cold;
+					T_top = T_tank + 0.33*dT_collector;
+					T_bot = T_tank - 0.67*dT_collector;
+					// T_top = T_hot
+					// T_bot = T_cold
+
+					hotLoss = U_tank * A_hot * (T_hot - T_room);
+					coldLoss = U_tank*A_cold*(T_cold - T_room);
+					Q_tankloss = hotLoss + coldLoss;
 				}
 
 				// calculate pumping losses (pump size is user entered) -
diff --git a/test/ssc_test/cmod_swh_test.cpp b/test/ssc_test/cmod_swh_test.cpp
index 96d868d1f..def0d0c76 100644
--- a/test/ssc_test/cmod_swh_test.cpp
+++ b/test/ssc_test/cmod_swh_test.cpp
@@ -43,7 +43,7 @@ TEST_F(CM_SWH, ResidentialDefault_cmod_swh) {
 
     ssc_number_t annual_energy;
     ssc_data_get_number(data, "annual_energy", &annual_energy);
-    EXPECT_NEAR(annual_energy, 2362.2, 0.1); 
+    EXPECT_NEAR(annual_energy, 2362.5, 0.1); 
 
 }