Merge remote-tracking branch 'origin/develop' into extendedHIspeedFix

design/FY2024/chiller_heater_part_load_fix.md

@@ -0,0 +1,106 @@
+DEFECT: Fix for Chiller Heater Always Assuming Evaporator is at Full Load
+================
+
+**Rick Strand, University of Illinois at Urbana-Champaign**
+
+ - Original Date: July 22, 2024
+ - Revision Date: July 30, 2024
+
+
+## Justification for New Feature ##
+
+The current heater mode portion of the chiller heater model in PlantCentralGSHP.cc is written with the built-in assumption that the evaporator is running at full load.  When the condenser load comes back at less than full load, the evaporator load is never adjusted and reports load and electric power at full load.  This is not correct and resulted in a defect being logged (Defect 10065).  This document is a plan for a potential solution to this problem.
+
+## E-mail and Conference Call Conclusions ##
+
+July 24, 2024: Discussed this on the technicalities call.  Decision was made to not implement an iteration strategy but to simply make an approximation of the PLR from the condenser load and then multiple full load evaporator load, compressor power, and false loading by that PLR.  Not ideal, but given all of the suspected problems in this model, it was decided to not invest too heavily in this now and turn this into a potential development topic in the future.
+
+## Overview ##
+
+The heat pump model for the ChillerHeater model is contained in the PlantCentralGSHP.cc file.  The model allows for various modes of operation: off (0), cooling only (1), heating only (2), heat recovery (3), cooling dominant (4), and heating dominant (5).  Off mode is obvious--no loads, nothing happening.  When in cooling or heating only mode, the heat rejection is lost/sent to the other thermal environment.  When in heat recovery mode, heat rejection at the condenser is used for heating purposes in the HVAC system.  Cooling and heating dominant modes have heat recovery but there is more heat recovery than is needed so the excess is rejected to whatever the outside environment is.  The cooling controlled modes (1, 3, and 4) are simulated using the CalcChillerModel in CentralPlantGSHP.cc.  The heating controlled modes (2 and 5) are simulated using the CalcChillerHeaterModel in CentralPlantGSHP.cc.  The cooling controlled modes seem to be working without any known issues.  The heating modes run the condenser to the correct load for heating.  However, a user noticed that there was an issue with the reported evaporator load and power consumption which always seemed to be relatively constant at a high level whenever the condenser was needed for a heating load.  This was traced back to the assumptions in the heater portion of the chiller heater model.
+
+The simlation flow in the heater portion of the model is summarized as follows.  Once it is identified that the chiller is in one of the heating modes, the chiller is assumed to run at full capacity to get an evaporator load, a compressor power consumption, and any false load for when the chiller is below the minimum part load ratio.  A condenser load at full power is calculated from this information.  This condenser load is then adjusted to fit the actual heating need.  From this condenser load, flows or temperatures are adjusted for this load as needed on the condenser side.  The simulation then moves on from there.
+
+The problem here is that the evaporator load and the compressor power are still at full load and are never adjusted when the condenser load gets reduced because the heating load does not require full load.  This is the source of the error--evaporator load and compressor power never change in heating mode regardless of the actual part load ratio based on the condenser load.  PLR simply stays at near 100%.  This is not correct and leads to over-estimation of both the evaporator load and the compressor power consumption.
+
+## Original Approach ##
+
+Note: before the actual fix takes place, it was decided to make a code improvement pass through the current chiller heater model.  This has already taken place and has been merged into develop.  The point was to make the code re-usable within the chiller heater model but it also realized some improvements in the cooling mode subroutine as well.  The changes took several different code sections and turned them into smaller subroutines.  The heating mode code is now much easier to follow, reducing the size of the routine by a factor of more than 3 (based on printouts of the routine before and after restructuring).  The real benefit will be seen when the problem is fixed as the algorithm should stay fairly compact and easy to follow (hopefully).
+
+The approach to this problem is to for the most part leave the initial pass at full load in tact.  The code still needs to know what the condenser is capable of from a heating standpoint.  If the required heating load (HeatingLoadToMeet) is larger than the condenser load at full power, then there is no change to be made.  However, when the condenser load is reduced because the HeatingLoadToMeet is less than the condenser load at full power, there needs to be additional code to handle this case.
+
+When the QCondenser is initially calculated based on full load evaporator conditions, a new variable (qCondOrig or something like that) tracking the original value will be set equal to this value.  After QCondenser is compared to HeatingLoadToMeet and then potentially modified by the new restructured subroutine adjustChillerHeaterFlowTemp, QCondenser will be compared to qCondOrig.  If QCondenser is lower than qCondOrig, then this means that we are no longer at full load and need to go into an adjustment procedure.
+
+Adjustment Procedure: The adjustment process to obtain a revised/appropriate evaporator load and compressor power will take place in the code after the initial pass to adjust the chiller heater flow rate and temperature.  The code in CalcChillerHeaterModel, after restructuring, looks like this:
+
+                if (CurrentMode == 2 || this->SimulHtgDominant) {
+                    if (CondMassFlowRate > DataBranchAirLoopPlant::MassFlowTolerance && CondDeltaTemp > 0.0) {
+                        this->adjustChillerHeaterFlowTemp(state, QCondenser, CondMassFlowRate, CondOutletTemp, CondInletTemp, CondDeltaTemp);
+                    } else {
+                        QCondenser = 0.0;
+                        CondOutletTemp = CondInletTemp;
+                    }
+                }
+
+New code will be added right after the call to adjustChillerHeaterFlowTemp.  There will need to be an iteration loop as there is no guarantee that one pass will obtain the correct evaporator conditions.  The iteration loop will end after 20 attempts or when the latest condenser load is sufficiently similar to the condenser load from the previous iteration.
+
+Before the iteration loop, set iteration number to zero and old condenser load to full load (qCondOrig) and new condenser load to existing condenser load.
+
+During the iteration loop (while iteration is less than 20 and difference between old and new QCondenser is outside of tolerance), the following functions will be done:
+
+Step 0: Update last iteration condenser load and iteration number
+
+Step 1: Calculate PLR based on comparison of last iteration condenser load to full condenser load, limit to MaxPartLoadRatio.
+
+Step 2: Calculate QEvaporator and EvapOutletTemp based on new PartLoadRatio.  Could potentially be turned into a new subroutine.
+
+Step 3: Call checkEvapOutletTemp. This in all likelihood won't do anything because it just makes sure that the minimum evaporator temperature limits are not violated, but it still needs to be run just in case.
+
+Step 4: Call calcPLRAndCyclingRatio.  This shouldn't change PLR because it was just used to calculate a new QEvaporator.  It could, however, result in some false loading if PLR drops below the minimum.
+
+Step 5: Recalculate CHPower, the compressor power using the same equation used from previously in this routine.  Could potentially be turned into a new function (one line).
+
+Step 6: Recalculate ActualCOP.  Could Could potentially be turned into a new function.
+
+Step 7: Calculate the new QCondenser for this iteration.  Reuse existing code that limits this based on minPLR, HeatingLoadToMeet.
+
+Step 8: Call adjustChillerHeaterFlowTemp to adjust flow rate and temperature if necessary.
+
+At this point, a new QCondenser has been calculated so the iteration cycle is done.  No additional code is needed after the iteration cycle as it should just be able to pick up where it left off as it currently does.
+
+## Modified Approach ##
+
+During the technicalities call, it was suggested that rather than iterating, we should just approximate the PLR from the condenser load and then multiply evaporator load and compressor power by this PLR.  The false load was also factored in this way though in this case it was probably zero at full load anyway.  Other problems in the algorithm were also fixed along the way.  No guarantees that this model is now 100% bug free but it should be improved.
+
+## Testing/Validation/Data Sources ##
+
+Testing will be done using the existing user input file that shows the problem.  Comparisons will be made between develop and the new version to establish that the results have changed after the fix has been implemented and that the new output makes sense.
+
+Unit testing: as this is being handled as part of a bug fix, at least two unit tests will be generated.  As there are several new subroutines, the unit tests will likely center on these new subroutines.
+
+## Input Output Reference Documentation ##
+
+No changes needed--this is an algorithm fix that does not require input changes.
+
+## Input Description ##
+
+No changes needed--this is an algorithm fix that does not require input changes.
+
+## Outputs Description ##
+
+No changes needed--this is an algorithm fix that does not require input changes or new output.
+
+## Engineering Reference ##
+
+Currently, the Engineering Reference has a section for the chiller heater model (ChillerHeaterPerformance\:Electric:EIR).  The subsection entitled "Heating-only mode and Simultaneous cooling-heating mode" essentially outlines the calculation process in heating only (2) and heating dominant (5) modes.  Additional text will be added to the end of this section to describe the work implemented as part of this fix that will outline the steps in a similar fashion to what is shown above in the Approach section of this document.
+
+## Example File and Transition Changes ##
+
+No transition changes are needed since there is no change to the input.  A change to an existing input file that has a chiller heater equipment may be needed to show differences in output, but it likely will without any changes to the .idf.  So, there are no changes anticipated to existing example files either.
+
+## References ##
+
+Current code in PlantCentralGSHP.cc
+
+
+

doc/engineering-reference/src/special-modules-reporting/resilience-metrics.tex

@@ -66,17 +66,69 @@ \subsubsection{Heat Index}\label{heat-index}
 \end{tabular}
 \end{table}
 
-Before version 24.2, the computation of the heat index is a refinement of a result obtained by
+The computation of the heat index is a refinement of a result obtained by
 multiple regression analysis carried out by Lans P. Rothfusz and described in a
-1990 National Weather Service (NWS) Technical Attachment (SR 90-23) [4-5].
+1990 National Weather Service (NWS) Technical Attachment (SR 90-23) [4-5]. The
+calculation is based on degree Fahrenheit.
 
-Starting from version 24.2, the heat index calculation adopts the extended heat
-index method developed by Lu \& Romps [17]. The previous heat index gives
-unrealistic results for very hot and humid or very cold and dry conditions. The
-extended index extends the domain of the heat index calculation to all
-combinations of temperature and relative humidity and gives a more realistic heat
-index for the extreme conditions. The implementation in EnergyPlus is based on
-the released Python code by Lu and Romps [18].
+The regression equation of Rothfusz is
+\begin{equation}  \label{eq:rm-1}
+HI = c_1 + c_2T + c_3R + c_4TR + c_5T^2 + c_6R^2 + c_7T^2R + c_8TR^2 + c_9T^2R^2
+\end{equation}
+
+where
+
+HI = heat index (expressed as an apparent temperature in degrees Fahrenheit),
+
+T = ambient dry-bulb temperature (in degrees Fahrenheit),
+
+R = relative humidity (percentage value between 0 and 100),
+
+$c_1$ = -42.379,
+
+$c_2$ = 2.04901523,
+
+$c_3$ = 10.14333127,
+
+$c_4$ = -0.22475541,
+
+$c_5$ = -0.00683783,
+
+$c_6$ = -0.05481717,
+
+$c_7$ = 0.00122874,
+
+$c_8$ = 0.00085282,
+
+$c_9$ = -0.00000199.
+
+If the RH is less than 13\% and the temperature is between 80 and \IP{112}{\fahrenheit}, then
+the following adjustment is subtracted from HI:
+
+\begin{equation}  \label{eq:rm-2}
+HI = (13 - R) / 4 * ((17 - |T - 95|) / 17)^{0.5}
+\end{equation}
+
+Otherwise, if the RH is greater than 85\% and the temperature is between 80 and
+\IP{87}{\fahrenheit}, then the following adjustment is added to HI:
+
+\begin{equation}  \label{eq:rm-3}
+HI = (R - 85) / 10 * (87 - T) / 5
+\end{equation}
+
+The Rothfusz regression is not appropriate when conditions of temperature and
+humidity warrant a heat index value below about \IP{80}{\fahrenheit}. In those cases, a simpler
+formula is applied to calculate values consistent with Steadman's results:
+
+\begin{equation}  \label{eq:rm-4}
+HI = 0.5 * (T + 61.0 + (T - 68.0) * 1.2 + (R * 0.094))
+\end{equation}
+
+In practice, the simple formula is computed first based on the temperature and
+humidity. If this heat index value is \IP{80}{\fahrenheit} or higher, the full regression
+equation along with any adjustment as described above is applied. The Rothfusz
+regression is not valid for extreme temperature and relative humidity conditions
+beyond the range of data considered by Steadman.
 
 The Heat Index Hours (accumulated hours for a space) and Heat Index
 OccupantHours (accumulated hours for the sum of all occupants in a space) of
@@ -464,11 +516,3 @@ \subsection{References}
 
 {[}16{]} ACGIH, Threshold Limit Values (TLVs) and Biological Exposure Indices
 (BEIs), 2012. doi:10.1073/pnas.0703993104.
-
-{[}17{]} Lu, Yi-Chuan, and David M. Romps. ``Extending the heat index''. Journal
-of Applied Meteorology and Climatology 61, no. 10 (2022): 1367-1383.
-doi:10.1073/pnas.0703993104.
-
-{[}18{]} Lu, Yi-Chuan, and David M. Romps. ``Lu and Romps, Extending the heat index, JAMC, 2022'',
-Physics of Climate, February 23, 2023.
-https://romps.berkeley.edu/papers/pubs-2020-heatindex.html.

src/EnergyPlus/AirflowNetwork/src/Solver.cpp

@@ -11228,13 +11228,15 @@ namespace AirflowNetwork {
                                     found = true;
                                 }
                             }
-                            if (!found) {
-                                ShowSevereError(m_state, format("{}Fan:ZoneExhaust is not defined in {}", RoutineName, CurrentModuleObject));
-                                ShowContinueError(m_state,
-                                                  "Zone Air Exhaust Node in ZoneHVAC:EquipmentConnections =" +
-                                                      m_state.dataLoopNodes->NodeID(m_state.dataZoneEquip->ZoneEquipConfig(j).ExhaustNode(k)));
-                                ErrorsFound = true;
-                            }
+                        }
+                        if (!found) {
+                            ShowSevereError(m_state, format("{}Fan:ZoneExhaust is not defined in {}", RoutineName, CurrentModuleObject));
+                            ShowContinueError(
+                                m_state,
+                                format("The inlet node of the {} Fan:ZoneExhaust is not defined in the {}'s ZoneHVAC:EquipmentConnections",
+                                       m_state.dataZoneEquip->ZoneEquipList(j).EquipName,
+                                       m_state.dataZoneEquip->ZoneEquipConfig(j).ZoneName));
+                            ErrorsFound = true;
                         }
                     }
                 }

src/EnergyPlus/CMakeLists.txt

@@ -262,8 +262,6 @@ set(SRC
     EvaporativeFluidCoolers.hh
     ExhaustAirSystemManager.cc
     ExhaustAirSystemManager.hh
-    ExtendedHI.cc
-    ExtendedHI.hh
     ExteriorEnergyUse.cc
     ExteriorEnergyUse.hh
     ExternalInterface.cc

src/EnergyPlus/Coils/CoilCoolingDXCurveFitOperatingMode.cc

@@ -291,7 +291,7 @@ void CoilCoolingDXCurveFitOperatingMode::CalcOperatingMode(EnergyPlus::EnergyPlu
     // Currently speedNum is 1-based, while this->speeds are zero-based
     auto &thisspeed(this->speeds[max(speedNum - 1, 0)]);
 
-    if (((speedNum == 1) && (PLR == 0.0)) || (inletNode.MassFlowRate == 0.0)) {
+    if ((speedNum == 0) || ((speedNum == 1) && (PLR == 0.0)) || (inletNode.MassFlowRate == 0.0)) {
         outletNode.Temp = inletNode.Temp;
         outletNode.HumRat = inletNode.HumRat;
         outletNode.Enthalpy = inletNode.Enthalpy;

src/EnergyPlus/Coils/CoilCoolingDXCurveFitSpeed.cc

@@ -838,7 +838,7 @@ Real64 CoilCoolingDXCurveFitSpeed::calcEffectiveSHR(const DataLoopNode::NodeData
     To1 = aa + Tcl;
     Error = 1.0;
     while (Error > 0.001) {
-        To2 = aa - Tcl * (std::exp(-To1 / Tcl) - 1.0);
+        To2 = aa - Tcl * std::expm1(-To1 / Tcl);
         Error = std::abs((To2 - To1) / To1);
         To1 = To2;
     }

src/EnergyPlus/ConvectionCoefficients.cc

@@ -6301,7 +6301,7 @@ Real64 CalcClearRoof(EnergyPlusData &state,
 
     Real64 Rf = RoughnessMultiplier[(int)RoughnessIndex];
     if (Rex > 0.1) { // avoid zero and crazy small denominators
-        Real64 tmp = std::log(1.0 + GrLn / pow_2(Rex));
+        Real64 tmp = std::log1p(GrLn / pow_2(Rex));
         eta = tmp / (1.0 + tmp);
     } else {
         eta = 1.0; // forced convection gone because no wind

src/EnergyPlus/DElightManagerF.cc

@@ -333,7 +333,7 @@ namespace DElightManagerF {
 
                     // Write each opaque bounding Surface to the DElight input file
                     for (int spaceNum : zn.spaceIndexes) {
-                        auto &thisSpace = state.dataHeatBal->space(spaceNum);
+                        auto const &thisSpace = state.dataHeatBal->space(spaceNum);
                         int const iSurfaceFirst = thisSpace.HTSurfaceFirst;
                         int const iSurfaceLast = thisSpace.HTSurfaceLast;
                         for (int isurf = iSurfaceFirst; isurf <= iSurfaceLast; ++isurf) {

src/EnergyPlus/DXCoils.cc

@@ -16787,7 +16787,7 @@ void CalcVRFCoolingCoil_FluidTCtrl(EnergyPlusData &state,
         }
 
         // If cycling fan, send coil part-load fraction to on/off fan via HVACDataGlobals
-        if (fanOp == HVAC::FanOp::Cycling) state.dataHVACGlobal->OnOffFanPartLoadFraction = PLF;
+        if (fanOp == HVAC::FanOp::Cycling) state.dataHVACGlobal->OnOffFanPartLoadFraction = thisDXCoil.CoolingCoilRuntimeFraction;
 
         // Check for saturation error and modify temperature at constant enthalpy
         if (OutletAirTemp < PsyTsatFnHPb(state, OutletAirEnthalpy, OutdoorPressure)) {

src/EnergyPlus/DataAirLoop.hh

@@ -85,7 +85,7 @@ namespace DataAirLoop {
         Array1D_int TermUnitCoolSizingIndex;       // Air terminal sizing numbers for zones cooled by this air loop
         Array1D_int TermUnitHeatSizingIndex;       // Air terminal sizing numbers for zones heated by this air loop
         Array1D<HVAC::AirDuctType> SupplyDuctType; // 1=main, 2=cooling, 3=heating, 4=other
-        EPVector<int> SupplyDuctBranchNum;         // Supply duct branch number
+        EPVector<int> SupplyDuctBranchNum;         // Supply duct branch number (airloop branchnum, not the actual branch index)
         EPVector<int> SupplyAirPathNum;            // Supply air path indexes
         EPVector<int> ReturnAirPathNum;            // Return air path indexes
     };

src/EnergyPlus/DataAirSystems.hh

@@ -190,7 +190,7 @@ namespace DataAirSystems {
         AirLoopMixerData Mixer;            // Data for mixer (if any)
         Array1D_bool ControlConverged;     // Convergence Parameter for controllers
         int NumOutletBranches = 0;
-        std::array<int, 3> OutletBranchNum = {0}; // branch numbers of system outlets
+        std::array<int, 3> OutletBranchNum = {0}; // airloop branch numbers of system outlets (not the actual branch index)
         int NumInletBranches = 0;
         std::array<int, 3> InletBranchNum = {0}; // branch number of system inlets
         bool CentralHeatCoilExists = true;       // true if there are central heating coils