From 29ef97fef663db873b38a51fd814d6c977ad9a3c Mon Sep 17 00:00:00 2001 From: Alex Swindler Date: Thu, 26 Feb 2026 15:38:15 -0700 Subject: [PATCH 1/3] GLHEC integration --- CMakeLists.txt | 2 +- dictionary.dic | 1 + .../heat-exchangers.tex | 186 ++++ idd/Energy+.idd.in | 24 +- src/EnergyPlus/CMakeLists.txt | 2 + .../GroundHeatExchangers/GLHEC/Model.cc | 837 ++++++++++++++++++ .../GroundHeatExchangers/GLHEC/Model.hh | 145 +++ src/EnergyPlus/GroundHeatExchangers/README.md | 5 + .../GroundHeatExchangers/Vertical.cc | 148 +++- .../GroundHeatExchangers/Vertical.hh | 16 + tst/EnergyPlus/unit/CMakeLists.txt | 2 + .../unit/GroundHeatExchangers/GLHEC.unit.cc | 110 +++ .../GroundHeatExchangers/GLHECModel.unit.cc | 450 ++++++++++ 13 files changed, 1911 insertions(+), 17 deletions(-) create mode 100644 src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc create mode 100644 src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.hh create mode 100644 tst/EnergyPlus/unit/GroundHeatExchangers/GLHEC.unit.cc create mode 100644 tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 95390daec95..a08ed977d09 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.19...4.0.3) # policy_max is set, so that will impliclty invoke cmake_policy(VERSION [...]) +cmake_minimum_required(VERSION 3.19...4.0.3) # policy_max is set, so that will implicitly invoke cmake_policy(VERSION [...]) # Use ccache if available, has to be before "project()" find_program(CCACHE_PROGRAM NAMES ccache sccache) diff --git a/dictionary.dic b/dictionary.dic index 81124bc4dda..739395e8544 100644 --- a/dictionary.dic +++ b/dictionary.dic @@ -10,6 +10,7 @@ DetailedGeometry Evaporative EvaporativeCooler FlatPlate +GLHEC HeatExchanger HeatExchangerAssisted HVACDOAS diff --git a/doc/engineering-reference/src/simulation-models-encyclopedic-reference-002/heat-exchangers.tex b/doc/engineering-reference/src/simulation-models-encyclopedic-reference-002/heat-exchangers.tex index b8d3cf2fc6f..f964c16e590 100644 --- a/doc/engineering-reference/src/simulation-models-encyclopedic-reference-002/heat-exchangers.tex +++ b/doc/engineering-reference/src/simulation-models-encyclopedic-reference-002/heat-exchangers.tex @@ -1136,6 +1136,192 @@ \subsubsection{Summary of Variable Short Time Step Response Factor Model}\label{ \textbf{End do} +\subsubsection{GLHEC Transient Segment Model (GLHECcalc)}\label{glhec-transient-segment-model-glheccalc} + +When \emph{GroundHeatExchanger:System} is configured with \emph{g-Function Calculation Method = \mbox{GLHECcalc}}, EnergyPlus uses a transient segment model for short-time borehole dynamics and a dynamic/sub-hour load aggregation scheme for long-term superposition. + +The GLHEC formulation uses: + +\begin{enumerate} +\item The same response-factor (\emph{g}-function) data infrastructure used by the vertical GHE framework. +\item A coupled ODE model of fluid and grout temperatures along borehole segments. +\item Dynamic temporal aggregation to limit history growth while preserving long-horizon response. +\end{enumerate} + +\paragraph{Borehole wall temperature superposition}\label{borehole-wall-temperature-superposition} + +At each simulation step \(n\), borehole wall temperature is computed from far-field ground temperature and superposed load history: + +\begin{equation} +T_{bw,n} = T_{\infty,n} + \frac{1}{2\pi k_s}\sum_{i=1}^{N_h} \Delta q_i g\left(\ln\left(\frac{\tau_i}{t_s}\right)\right) +\end{equation} + +where: + +\begin{equation} +t_s = \frac{H^2}{9\alpha_s} +\end{equation} + +\(T_{bw,n}\) is borehole wall temperature for the current time step (\(^\circ\)C) + +\(T_{\infty,n}\) is the far-field ground temperature (\(^\circ\)C) + +\(k_s\) is soil thermal conductivity (W/m-K) + +\(H\) is borehole length (m) + +\(\alpha_s\) is soil thermal diffusivity (m\(^2\)/s) + +\(\Delta q_i\) is the change in aggregated load rate for history bin \(i\) (W/m) + +\(\tau_i\) is elapsed time between the current step and the implementation-defined cumulative time location for history bin \(i\) (s). + +\paragraph{Dynamic and sub-hour aggregation}\label{dynamic-and-sub-hour-aggregation} + +GLHEC uses two aggregation levels: + +\begin{enumerate} +\item A \textbf{sub-hour} tracker that retains irregular bins over the first hour. +\item A \textbf{dynamic} bin stack beginning at 1 hour, with level width growth controlled by user inputs. +\end{enumerate} + +Dynamic bins are stored from oldest to youngest. For each dynamic bin \(i\), with stored energy \(E_i\) (J/m) and width \(\Delta t_i\) (s), the shifted amount over current step \(\Delta t_n\) is: + +\begin{equation} +\delta_i = E_i \frac{\Delta t_n}{\Delta t_i} +\end{equation} + +and bin update is: + +\begin{equation} +E_i^{+} = E_i - \delta_i + \delta_{i+1}, \qquad i=0,\ldots,N_d-2,\qquad \delta_0 = 0 +\end{equation} + +\begin{equation} +E_{N_d-1}^{+} = E_{N_d-1} - \delta_{N_d-1} + E_{\text{subhr}\rightarrow\text{dyn}} +\end{equation} + +The sub-hour tracker exports energy older than one hour into the youngest dynamic bin. Superposition uses \(q_i = E_i/\Delta t_i\), then differences \(\Delta q_i = q_i - q_{i-1}\). + +\paragraph{Per-borehole flow and energy bookkeeping}\label{per-borehole-flow-and-energy-bookkeeping} + +Total loop flow is split uniformly across boreholes: + +\begin{equation} +\dot m_{bh} = \frac{\dot m_{loop}}{N_{bh}} +\end{equation} + +Loop-side heat transfer is: + +\begin{equation} +\dot Q_{loop} = \dot m_{loop} c_p(T_{in}) \left(T_{in} - T_{out}\right) +\end{equation} + +Segment borehole heat transfer is accumulated per borehole and then scaled by number of boreholes: + +\begin{equation} +\dot Q_{bh,total} = N_{bh}\sum_{s=1}^{N_{seg}} \dot Q_{bh,s} +\end{equation} + +Energy stored in history for the current time step is normalized by borehole length: + +\begin{equation} +E_n = \left(\frac{\dot Q_{bh,total}/N_{bh}}{H}\right)\Delta t_n +\end{equation} + +\paragraph{Segment ODE system}\label{segment-ode-system} + +Each borehole segment solves five coupled states: + +\begin{equation} +\left[T_{f,1},T_{f,2},T_{g,c},T_{g,1},T_{g,2}\right] +\end{equation} + +where \(T_{f,1}\) and \(T_{f,2}\) are fluid temperatures in the down/up legs, \(T_{g,c}\) is the central grout node, and \(T_{g,1}, T_{g,2}\) are grout nodes adjacent to each leg. + +The fluid-leg equations are: + +\begin{equation} +\frac{dT_{f,1}}{dt} = \frac{\left(T_{in,1}-T_{f,1}\right)/R_f + \left(T_{g,c}-T_{f,1}\right)L_s/(R_{12}/2) + \left(T_{g,1}-T_{f,1}\right)L_s/R_b}{C_f} +\end{equation} + +\begin{equation} +\frac{dT_{f,2}}{dt} = \frac{\left(T_{in,2}-T_{f,2}\right)/R_f + \left(T_{g,c}-T_{f,2}\right)L_s/(R_{12}/2) + \left(T_{g,2}-T_{f,2}\right)L_s/R_b}{C_f} +\end{equation} + +The grout-node equations are: + +\begin{equation} +\frac{dT_{g,c}}{dt} = \frac{\left(T_{f,1}-T_{g,c}\right)L_s/(R_{12}/2) + \left(T_{f,2}-T_{g,c}\right)L_s/(R_{12}/2)}{C_{g,1}} +\end{equation} + +\begin{equation} +\frac{dT_{g,1}}{dt} = \frac{\left(T_{f,1}-T_{g,1}\right)L_s/R_b + \left(T_{bw}-T_{g,1}\right)L_s/R_b}{C_{g,2}} +\end{equation} + +\begin{equation} +\frac{dT_{g,2}}{dt} = \frac{\left(T_{f,2}-T_{g,2}\right)L_s/R_b + \left(T_{bw}-T_{g,2}\right)L_s/R_b}{C_{g,2}} +\end{equation} + +where: + +\(L_s\) is segment length (m) + +\(R_f = 1/(\dot m_{bh} c_p)\) is fluid flow resistance term + +\(R_b\) is average borehole resistance (m-K/W) + +\(R_{12}\) is direct-coupling resistance between U-tube legs (m-K/W) + +\(C_f\), \(C_{g,1}\), and \(C_{g,2}\) are fluid and grout/pipe effective thermal capacitances (J/K). + +\paragraph{Pipe and borehole resistances}\label{pipe-and-borehole-resistances-glhec} + +Pipe resistance combines convection and pipe-wall conduction: + +\begin{equation} +R_p = \frac{1}{Nu \pi k_f} + \frac{\ln\left(d_o/d_i\right)}{2\pi k_p} +\end{equation} + +where \(Nu\) is computed from Reynolds/Prandtl number correlations with laminar-transition-turbulent smoothing. + +The borehole resistance and direct-coupling resistance use the multipole-style relations in Claesson and Hellstr\"om (2011), parameterized by: + +\begin{equation} +\theta_1 = \frac{s}{D_b}, \qquad \theta_2 = \frac{D_b}{d_o}, \qquad \theta_3 = \frac{1}{\theta_1\theta_2}, \qquad +\sigma = \frac{k_g-k_s}{k_g+k_s} +\end{equation} + +with \(s\) shank spacing, \(D_b\) borehole diameter, \(d_o\) outer pipe diameter, \(k_g\) grout conductivity, and \(k_s\) soil conductivity. + +\paragraph{Segment coupling iteration and integration}\label{segment-coupling-iteration-and-integration} + +Within each plant/system time step, segment equations are solved iteratively to couple downflow and upflow legs through the U-bend boundary condition. A fixed-step fourth-order Runge-Kutta (RK4) integrator advances each segment state. + +\subsubsection{GLHEC User Controls and Solver Defaults}\label{glhec-user-controls-and-solver-defaults} + +Table~\ref{tab:glhec-user-controls} summarizes user-facing controls in \emph{GroundHeatExchanger:System} and core solver defaults currently used by the GLHEC implementation. + +\begin{table}[hbtp] +\centering +\caption{GLHEC User Controls and Defaults\protect\label{tab:glhec-user-controls}} +\begin{tabular}{p{0.38\textwidth}p{0.18\textwidth}p{0.38\textwidth}} +\toprule +\textbf{Input/Setting} & \textbf{Default} & \textbf{Effect} \\ +\midrule +g-Function Calculation Method = \mbox{GLHECcalc} & N/A & Activates GLHEC transient segment model path. \\ +GLHEC Number of Borehole Segments & 10 & Number of axial thermal segments per borehole. \\ +GLHEC Number of Iterations & 2 & Number of within-step segment-coupling iterations. \\ +GLHEC Grout Fraction & 0.5 & Splits grout/pipe thermal mass between central and side grout nodes. \\ +GLHEC Aggregation Expansion Rate & 1.62 & Dynamic-bin width growth factor between aggregation levels. \\ +GLHEC Aggregation Bins Per Level & 9 & Number of bins before dynamic-bin width is expanded. \\ +RK4 substep control (internal) & 100 substeps per plant/system step & Fixed-step RK4 advancement for each segment state. \\ +\bottomrule +\end{tabular} +\end{table} + +When total loop mass flow is effectively zero, GLHEC returns \(T_{out}=T_{in}\) with zero heat transfer for the current step while retaining borehole wall state history. + \subsubsection{References}\label{references-2-006} \hangindent=2em diff --git a/idd/Energy+.idd.in b/idd/Energy+.idd.in index ed80e986352..3ab3a913bff 100644 --- a/idd/Energy+.idd.in +++ b/idd/Energy+.idd.in @@ -79979,6 +79979,7 @@ GroundHeatExchanger:System, \key UHFcalc \key UBHWTcalc \key FullDesign + \key GLHECcalc \default UHFcalc A8, \field GHE:Vertical:Sizing Object Type \type choice @@ -80288,9 +80289,30 @@ GroundHeatExchanger:System, A109, \field GHE:Vertical:Single Object Name 99 \type object-list \object-list GroundHeatExchangerVerticalSingleNames - A110; \field GHE:Vertical:Single Object Name 100 + A110, \field GHE:Vertical:Single Object Name 100 \type object-list \object-list GroundHeatExchangerVerticalSingleNames + N111, \field GLHEC Number of Borehole Segments + \type integer + \minimum 1 + \default 10 + N112, \field GLHEC Number of Iterations + \type integer + \minimum 1 + \default 2 + N113, \field GLHEC Grout Fraction + \type real + \minimum 0.0 + \maximum 1.0 + \default 0.5 + N114, \field GLHEC Aggregation Expansion Rate + \type real + \minimum 1.0 + \default 1.62 + N115; \field GLHEC Aggregation Bins Per Level + \type integer + \minimum 1 + \default 9 GroundHeatExchanger:Vertical:Sizing:Rectangle, \memo Specifies parameters to be used for Ground Heat Exchanger borehole diff --git a/src/EnergyPlus/CMakeLists.txt b/src/EnergyPlus/CMakeLists.txt index 7d412c8036c..de0897a19a6 100644 --- a/src/EnergyPlus/CMakeLists.txt +++ b/src/EnergyPlus/CMakeLists.txt @@ -309,6 +309,8 @@ set(SRC GroundHeatExchangers/BoreholeArray.hh GroundHeatExchangers/BoreholeSingle.cc GroundHeatExchangers/BoreholeSingle.hh + GroundHeatExchangers/GLHEC/Model.cc + GroundHeatExchangers/GLHEC/Model.hh GroundHeatExchangers/Properties.cc GroundHeatExchangers/Properties.hh GroundHeatExchangers/ResponseFactors.cc diff --git a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc new file mode 100644 index 00000000000..ef3868b2e49 --- /dev/null +++ b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc @@ -0,0 +1,837 @@ +// EnergyPlus, Copyright (c) 1996-present, The Board of Trustees of the University of Illinois, +// The Regents of the University of California, through Lawrence Berkeley National Laboratory +// (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge +// National Laboratory, managed by UT-Battelle, Alliance for Energy Innovation, LLC, and other +// contributors. All rights reserved. +// +// NOTICE: This Software was developed under funding from the U.S. Department of Energy and the +// U.S. Government consequently retains certain rights. As such, the U.S. Government has been +// granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable, +// worldwide license in the Software to reproduce, distribute copies to the public, prepare +// derivative works, and perform publicly and display publicly, and to permit others to do so. +// +// Redistribution and use in source and binary forms, with or without modification, are permitted +// provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this list of +// conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, this list of +// conditions and the following disclaimer in the documentation and/or other materials +// provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, +// the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without specific prior +// written permission. +// +// (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form +// without changes from the version obtained under this License, or (ii) Licensee makes a +// reference solely to the software portion of its product, Licensee must refer to the +// software as "EnergyPlus version X" software, where "X" is the version number Licensee +// obtained under this License and may not use a different name for the software. Except as +// specifically required in this Section (4), Licensee shall not use in a company name, a +// product name, in advertising, publicity, or other promotional activities any name, trade +// name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly +// similar designation, without the U.S. Department of Energy's prior written consent. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY +// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +#include + +#include +#include +#include +#include +#include + +#include + +namespace EnergyPlus::GroundHeatExchangers::GLHEC { + +namespace { + + constexpr Real64 small = 1.0e-12; + + template T clampMin(T value, T minimum) + { + return (value < minimum) ? minimum : value; + } + + Real64 linInterp(Real64 x, Real64 xLow, Real64 xHigh, Real64 yLow, Real64 yHigh) + { + if (std::abs(xHigh - xLow) < small) { + return yLow; + } + return (x - xLow) / (xHigh - xLow) * (yHigh - yLow) + yLow; + } + + class LinearInterpolator + { + public: + void setXY(std::vector xVals, std::vector yVals) + { + x = std::move(xVals); + y = std::move(yVals); + } + + [[nodiscard]] Real64 interpolate(Real64 q) const + { + if (x.empty() || y.empty()) { + return 0.0; + } + if (x.size() == 1) { + return y[0]; + } + if (q <= x.front()) { + return y.front(); + } + if (q >= x.back()) { + return y.back(); + } + auto const upper = std::upper_bound(x.begin(), x.end(), q); + std::size_t const idxH = static_cast(std::distance(x.begin(), upper)); + std::size_t const idxL = idxH - 1; + return linInterp(q, x[idxL], x[idxH], y[idxL], y[idxH]); + } + + [[nodiscard]] std::vector interpolate(std::vector const &queries) const + { + std::vector values; + values.reserve(queries.size()); + if (x.empty() || y.empty()) { + values.assign(queries.size(), 0.0); + return values; + } + if (x.size() == 1) { + values.assign(queries.size(), y[0]); + return values; + } + + for (Real64 const q : queries) { + if (q < x.front()) { + values.push_back(linInterp(q, x[0], x[1], y[0], y[1])); + } else if (q >= x.back()) { + std::size_t const n = x.size(); + values.push_back(linInterp(q, x[n - 2], x[n - 1], y[n - 2], y[n - 1])); + } else { + auto const upper = std::upper_bound(x.begin(), x.end(), q); + std::size_t const idxH = static_cast(std::distance(x.begin(), upper)); + std::size_t const idxL = idxH - 1; + values.push_back(linInterp(q, x[idxL], x[idxH], y[idxL], y[idxH])); + } + } + return values; + } + + private: + std::vector x; + std::vector y; + }; + + template void removeIndices(std::vector &values, std::vector const &indicesToRemove) + { + if (indicesToRemove.empty() || values.empty()) { + return; + } + + std::vector removeMask(values.size(), false); + for (std::size_t const idx : indicesToRemove) { + if (idx < values.size()) { + removeMask[idx] = true; + } + } + + std::vector filtered; + filtered.reserve(values.size()); + for (std::size_t i = 0; i < values.size(); ++i) { + if (!removeMask[i]) { + filtered.push_back(values[i]); + } + } + values = std::move(filtered); + } + + class SubHourAggregation + { + public: + void setup(std::vector const &lntts, std::vector const &gValues, std::vector const &gBValues, Real64 timeScale) + { + interpG.setXY(lntts, gValues); + interpGb.setXY(lntts, gBValues.empty() ? gValues : gBValues); + ts = clampMin(timeScale, small); + reset(); + } + + void reset() + { + energy.clear(); + dts.clear(); + energy.push_back(0.0); + dts.push_back(secondsPerHour); + prevUpdateTime = 0; + } + + [[nodiscard]] Real64 aggregate(int timeSeconds, Real64 energyJPerMeter) + { + if (timeSeconds <= prevUpdateTime) { + return 0.0; + } + + energy.push_back(energyJPerMeter); + dts.push_back(timeSeconds - prevUpdateTime); + + std::vector dtsFlipped = dts; + std::reverse(dtsFlipped.begin(), dtsFlipped.end()); + + std::vector dtUpper(dts.size(), 0.0); + Real64 cumulative = 0.0; + for (std::size_t i = 0; i < dtsFlipped.size(); ++i) { + cumulative += static_cast(dtsFlipped[i]); + dtUpper[i] = cumulative; + } + std::reverse(dtUpper.begin(), dtUpper.end()); + + std::vector dtLower(dts.size(), 0.0); + for (std::size_t i = 0; i < dts.size(); ++i) { + dtLower[i] = dtUpper[i] - static_cast(dts[i]); + } + + std::vector idxFull; + std::vector idxPartial; + for (std::size_t i = 0; i < dtLower.size(); ++i) { + if (dtLower[i] >= secondsPerHour) { + idxFull.push_back(i); + } else if (dtUpper[i] > secondsPerHour && dtLower[i] < secondsPerHour) { + idxPartial.push_back(i); + } + } + + Real64 loadToShift = 0.0; + for (std::size_t const i : idxFull) { + loadToShift += energy[i]; + } + + if (!idxPartial.empty()) { + std::size_t const idx = idxPartial.front(); + Real64 const upperEdge = dtUpper[idx]; + Real64 const lowerEdge = dtLower[idx]; + Real64 const denom = clampMin(upperEdge - lowerEdge, small); + Real64 const fraction = (upperEdge - secondsPerHour) / denom; + loadToShift += fraction * energy[idx]; + energy[idx] = (1.0 - fraction) * energy[idx]; + dts[idx] = std::max(1, static_cast(std::llround((1.0 - fraction) * dts[idx]))); + } + + removeIndices(energy, idxFull); + removeIndices(dts, idxFull); + prevUpdateTime = timeSeconds; + return loadToShift; + } + + [[nodiscard]] std::vector getLoadRates() const + { + std::vector rates; + rates.reserve(energy.size()); + for (std::size_t i = 0; i < energy.size(); ++i) { + rates.push_back(energy[i] / clampMin(static_cast(dts[i]), small)); + } + return rates; + } + + [[nodiscard]] std::vector const &getTimeBins() const + { + return dts; + } + + [[nodiscard]] int getPrevUpdateTime() const + { + return prevUpdateTime; + } + + private: + static constexpr int secondsPerHour = 3600; + LinearInterpolator interpG; + LinearInterpolator interpGb; + Real64 ts = 1.0; + int prevUpdateTime = 0; + std::vector energy; + std::vector dts; + }; + + class DynamicAggregation + { + public: + void setup(ModelConfig const &config, Real64 timeScale) + { + interpG.setXY(config.lntts, config.gValues); + std::vector gB = config.gBValues; + if (gB.empty()) { + gB = config.gValues; + } + interpGb.setXY(config.lntts, gB); + subHour.setup(config.lntts, config.gValues, gB, timeScale); + ts = clampMin(timeScale, small); + expansionRate = std::max(1.0, config.aggregation.expansionRate); + binsPerLevel = std::max(1, config.aggregation.binsPerLevel); + simulationHorizonSeconds = std::max(static_cast(secondsPerHour), config.aggregation.simulationHorizonSeconds); + dts.clear(); + energy.clear(); + + Real64 dt = static_cast(secondsPerHour); + Real64 t = static_cast(secondsPerHour); + while (true) { + for (int i = 0; i < binsPerLevel; ++i) { + t += dt; + energy.insert(energy.begin(), 0.0); + dts.insert(dts.begin(), std::max(1, static_cast(std::llround(dt)))); + if (t >= simulationHorizonSeconds) { + reset(); + return; + } + } + dt *= expansionRate; + } + } + + void reset() + { + std::fill(energy.begin(), energy.end(), 0.0); + subHour.reset(); + prevUpdateTime = 0; + } + + void aggregate(int timeSeconds, Real64 energyJPerMeter) + { + if (timeSeconds < prevUpdateTime) { + reset(); + } + if (timeSeconds <= prevUpdateTime) { + return; + } + + Real64 const shiftedFromSubHour = subHour.aggregate(timeSeconds, energyJPerMeter); + int const elapsedTime = timeSeconds - prevUpdateTime; + + std::vector fractionShift; + fractionShift.reserve(dts.size()); + for (int const dt : dts) { + fractionShift.push_back(static_cast(elapsedTime) / static_cast(dt)); + } + if (!fractionShift.empty()) { + fractionShift[0] = 0.0; + } + + std::vector delta; + delta.reserve(energy.size()); + for (std::size_t i = 0; i < energy.size(); ++i) { + Real64 const thisDelta = energy[i] * fractionShift[i]; + delta.push_back(thisDelta); + energy[i] -= thisDelta; + } + + std::vector rolledDelta(delta.size(), 0.0); + for (std::size_t i = 0; i < delta.size(); ++i) { + rolledDelta[i] = delta[(i + 1) % delta.size()]; + } + + for (std::size_t i = 0; i < energy.size(); ++i) { + energy[i] += rolledDelta[i]; + } + if (!energy.empty()) { + energy.back() += shiftedFromSubHour; + } + + prevUpdateTime = timeSeconds; + } + + [[nodiscard]] std::pair temporalSuperposition(int dtSeconds) const + { + if (dts.empty() || dtSeconds <= 0) { + return {0.0, 0.0}; + } + + std::vector longTermRates; + longTermRates.reserve(energy.size()); + for (std::size_t i = 0; i < energy.size(); ++i) { + longTermRates.push_back(energy[i] / clampMin(static_cast(dts[i]), small)); + } + + std::vector shortTermRates = subHour.getLoadRates(); + + std::vector q = longTermRates; + q.insert(q.end(), shortTermRates.begin(), shortTermRates.end()); + if (q.empty()) { + return {0.0, 0.0}; + } + + std::vector dq(q.size(), 0.0); + dq[0] = q[0]; + for (std::size_t i = 1; i < q.size(); ++i) { + dq[i] = q[i] - q[i - 1]; + } + + std::vector allDts = dts; + std::vector const &subHourDts = subHour.getTimeBins(); + allDts.insert(allDts.end(), subHourDts.begin(), subHourDts.end()); + allDts.push_back(dtSeconds); + + std::vector reversedDts = allDts; + std::reverse(reversedDts.begin(), reversedDts.end()); + + std::vector cumulativeTimes(reversedDts.size(), 0.0); + Real64 cumulative = 0.0; + for (std::size_t i = 0; i < reversedDts.size(); ++i) { + cumulative += static_cast(reversedDts[i]); + cumulativeTimes[i] = cumulative; + } + std::reverse(cumulativeTimes.begin(), cumulativeTimes.end()); + cumulativeTimes.pop_back(); + + std::vector lntts; + lntts.reserve(cumulativeTimes.size()); + for (Real64 const t : cumulativeTimes) { + lntts.push_back(std::log(clampMin(t / ts, small))); + } + + std::vector g = interpG.interpolate(lntts); + std::vector gB = interpGb.interpolate(lntts); + + Real64 dotG = 0.0; + Real64 dotGb = 0.0; + for (std::size_t i = 0; i < dq.size(); ++i) { + dotG += dq[i] * g[i]; + dotGb += dq[i] * gB[i]; + } + return {dotG, dotGb}; + } + + private: + static constexpr int secondsPerHour = 3600; + LinearInterpolator interpG; + LinearInterpolator interpGb; + SubHourAggregation subHour; + Real64 ts = 1.0; + Real64 expansionRate = 1.62; + int binsPerLevel = 9; + Real64 simulationHorizonSeconds = 50.0 * 365.0 * 24.0 * static_cast(secondsPerHour); + int prevUpdateTime = 0; + std::vector energy; + std::vector dts; + }; + + struct SegmentInputs + { + Real64 flowRate = 0.0; + Real64 inletTemp1 = 0.0; + Real64 inletTemp2 = 0.0; + Real64 boundaryTemp = 0.0; + Real64 bhResistance = 0.0; + Real64 dcResistance = 0.0; + }; + + class Segment + { + public: + Segment(ModelConfig const &config, Real64 segmentLength) + : length(segmentLength), diameter(config.boreholeDiameter), groutFraction(config.groutFraction), groutDensity(config.groutDensity), + groutSpecificHeat(config.groutSpecificHeat), pipeDensity(config.pipeDensity), pipeSpecificHeat(config.pipeSpecificHeat), + pipeInnerDiameter(config.pipeInnerDiameter), pipeOuterDiameter(config.pipeOuterDiameter), odeSettings(config.ode) + { + reset(20.0); + } + + void reset(Real64 temp) + { + y = {temp, temp, temp, temp, temp}; + heatRate = 0.0; + } + + void step(Real64 timeStepSeconds, SegmentInputs const &inputs, FluidPropertyFunctions const &fluidProps) + { + if (timeStepSeconds <= 0.0) { + return; + } + Context ctx; + ctx.seg = this; + ctx.inputs = inputs; + ctx.fluidCp = clampMin(fluidProps.cp(inputs.inletTemp1), small); + ctx.fluidRho = clampMin(fluidProps.rho(inputs.inletTemp1), small); + + rk4Step(timeStepSeconds, ctx); + Real64 const rB = clampMin(inputs.bhResistance, small); + heatRate = ((y[3] - inputs.boundaryTemp) / rB + (y[4] - inputs.boundaryTemp) / rB) * length; + } + + [[nodiscard]] Real64 outlet1() const + { + return y[0]; + } + + [[nodiscard]] Real64 outlet2() const + { + return y[1]; + } + + [[nodiscard]] Real64 boreholeHeatRate() const + { + return heatRate; + } + + private: + struct Context + { + Segment *seg = nullptr; + SegmentInputs inputs; + Real64 fluidCp = 0.0; + Real64 fluidRho = 0.0; + }; + + void rhs(Context const &ctx, double const yvals[], double dydt[]) const + { + Real64 const rF = 1.0 / clampMin(ctx.inputs.flowRate * ctx.fluidCp, small); + Real64 const rB = clampMin(ctx.inputs.bhResistance, small); + Real64 const r12 = clampMin(ctx.inputs.dcResistance, small); + Real64 const cF = clampMin(ctx.fluidRho * ctx.fluidCp * pipeFluidVolume(), small); + Real64 const cG1 = + clampMin(groutFraction * groutSpecificHeat * groutDensity * groutVolume() + pipeSpecificHeat * pipeDensity * pipeWallVolume(), small); + Real64 const cG2 = clampMin( + ((1.0 - groutFraction) * groutSpecificHeat * groutDensity * groutVolume() + pipeSpecificHeat * pipeDensity * pipeWallVolume()) / 2.0, + small); + + dydt[0] = + ((ctx.inputs.inletTemp1 - yvals[0]) / rF + (yvals[2] - yvals[0]) * length / (r12 / 2.0) + (yvals[3] - yvals[0]) * length / rB) / cF; + dydt[1] = + ((ctx.inputs.inletTemp2 - yvals[1]) / rF + (yvals[2] - yvals[1]) * length / (r12 / 2.0) + (yvals[4] - yvals[1]) * length / rB) / cF; + dydt[2] = ((yvals[0] - yvals[2]) * length / (r12 / 2.0) + (yvals[1] - yvals[2]) * length / (r12 / 2.0)) / cG1; + dydt[3] = ((yvals[0] - yvals[3]) * length / rB + (ctx.inputs.boundaryTemp - yvals[3]) * length / rB) / cG2; + dydt[4] = ((yvals[1] - yvals[4]) * length / rB + (ctx.inputs.boundaryTemp - yvals[4]) * length / rB) / cG2; + } + + void rk4Step(Real64 timeStepSeconds, Context const &ctx) + { + int const n = std::max(20, std::max(1, odeSettings.stepsPerTimeStep) * 5); + Real64 const h = timeStepSeconds / n; + std::array k1{}, k2{}, k3{}, k4{}, y2{}, y3{}, y4{}; + for (int i = 0; i < n; ++i) { + rhs(ctx, y.data(), k1.data()); + for (std::size_t j = 0; j < y.size(); ++j) { + y2[j] = y[j] + h * 0.5 * k1[j]; + } + rhs(ctx, y2.data(), k2.data()); + for (std::size_t j = 0; j < y.size(); ++j) { + y3[j] = y[j] + h * 0.5 * k2[j]; + } + rhs(ctx, y3.data(), k3.data()); + for (std::size_t j = 0; j < y.size(); ++j) { + y4[j] = y[j] + h * k3[j]; + } + rhs(ctx, y4.data(), k4.data()); + for (std::size_t j = 0; j < y.size(); ++j) { + y[j] += h * (k1[j] + 2.0 * k2[j] + 2.0 * k3[j] + k4[j]) / 6.0; + } + } + } + + [[nodiscard]] Real64 pipeFluidVolume() const + { + return Constant::Pi * std::pow(pipeInnerDiameter, 2) / 4.0 * length; + } + + [[nodiscard]] Real64 pipeWallVolume() const + { + return (Constant::Pi * std::pow(pipeOuterDiameter, 2) / 4.0 - Constant::Pi * std::pow(pipeInnerDiameter, 2) / 4.0) * length; + } + + [[nodiscard]] Real64 groutVolume() const + { + Real64 const seg = Constant::Pi * std::pow(diameter, 2) / 4.0 * length; + Real64 const pipe = 2.0 * Constant::Pi * std::pow(pipeOuterDiameter, 2) / 4.0 * length; + return std::max(0.0, seg - pipe); + } + + Real64 length = 0.0; + Real64 diameter = 0.0; + Real64 groutFraction = 0.5; + Real64 groutDensity = 0.0; + Real64 groutSpecificHeat = 0.0; + Real64 pipeDensity = 0.0; + Real64 pipeSpecificHeat = 0.0; + Real64 pipeInnerDiameter = 0.0; + Real64 pipeOuterDiameter = 0.0; + OdeSolverSettings odeSettings{}; + std::array y{}; + Real64 heatRate = 0.0; + }; + +} // namespace + +struct Model::Impl +{ + explicit Impl(ModelConfig cfg, FluidPropertyFunctions props) : config(std::move(cfg)), fluidProps(std::move(props)) + { + if (config.lntts.empty() || config.gValues.empty() || config.lntts.size() != config.gValues.size()) { + throw std::runtime_error("GLHEC requires g-function data"); + } + if (!config.gBValues.empty() && config.gBValues.size() != config.lntts.size()) { + throw std::runtime_error("GLHEC gB-function data size mismatch"); + } + if (!fluidProps.cp || !fluidProps.rho || !fluidProps.viscosity || !fluidProps.conductivity) { + throw std::runtime_error("GLHEC requires fluid property callbacks"); + } + config.numSegments = std::max(1, config.numSegments); + config.numIterations = std::max(1, config.numIterations); + config.numBoreholes = std::max(1u, config.numBoreholes); + if (config.boreholeLength <= 0.0 || config.boreholeDiameter <= 0.0 || config.pipeInnerDiameter <= 0.0 || + config.pipeOuterDiameter <= config.pipeInnerDiameter) { + throw std::runtime_error("GLHEC geometry inputs are invalid"); + } + int const nSeg = config.numSegments; + Real64 const segLength = config.boreholeLength / static_cast(nSeg); + for (int i = 0; i < nSeg; ++i) { + segments.emplace_back(config, segLength); + } + theta1 = config.shankSpacing / clampMin(config.boreholeDiameter, small); + theta2 = config.boreholeDiameter / clampMin(config.pipeOuterDiameter, small); + theta3 = 1.0 / clampMin(theta1 * theta2, small); + sigma = (config.groutConductivity - config.soilConductivity) / clampMin(config.groutConductivity + config.soilConductivity, small); + ts = std::pow(config.boreholeLength, 2) / clampMin(9.0 * config.soilDiffusivity, small); + c0 = 1.0 / (2.0 * Constant::Pi * clampMin(config.soilConductivity, small)); + aggregation.setup(config, ts); + } + + void reset(Real64 temp) + { + lastSimTime = -1; + haveStepStartSnapshot = false; + outletTemp = temp; + bhWallTemp = temp; + energy = 0.0; + uBendTemp = temp; + aggregation.reset(); + for (auto &seg : segments) { + seg.reset(temp); + } + } + + [[nodiscard]] ModelStepOutputs simulate(ModelStepInputs const &inputs) + { + if (inputs.timeStepSeconds <= 0) { + ModelStepOutputs outputs; + outputs.outletTemp = inputs.inletTemp; + outputs.heatRate = 0.0; + outputs.boreholeHeatRate = 0.0; + outputs.boreholeWallTemp = bhWallTemp; + outputs.avgFluidTemp = inputs.inletTemp; + return outputs; + } + + if (inputs.timeSeconds < lastSimTime) { + reset(inputs.farFieldGroundTemp); + } + + if (lastSimTime < inputs.timeSeconds) { + aggregationAtStepStart = aggregation; + segmentsAtStepStart = segments; + energyAtStepStart = energy; + uBendTempAtStepStart = uBendTemp; + outletTempAtStepStart = outletTemp; + bhWallTempAtStepStart = bhWallTemp; + haveStepStartSnapshot = true; + lastSimTime = inputs.timeSeconds; + } else if (lastSimTime == inputs.timeSeconds && haveStepStartSnapshot) { + aggregation = aggregationAtStepStart; + segments = segmentsAtStepStart; + energy = energyAtStepStart; + uBendTemp = uBendTempAtStepStart; + outletTemp = outletTempAtStepStart; + bhWallTemp = bhWallTempAtStepStart; + } + + aggregation.aggregate(inputs.timeSeconds, energy); + auto const hist = aggregation.temporalSuperposition(inputs.timeStepSeconds); + bhWallTemp = inputs.farFieldGroundTemp + hist.first * c0; + + Real64 const numBoreholes = static_cast(std::max(1u, config.numBoreholes)); + Real64 const massFlowRatePerBorehole = inputs.massFlowRate / numBoreholes; + + if (massFlowRatePerBorehole <= small) { + energy = 0.0; + ModelStepOutputs outputs; + outputs.outletTemp = inputs.inletTemp; + outputs.heatRate = 0.0; + outputs.boreholeHeatRate = 0.0; + outputs.boreholeWallTemp = bhWallTemp; + outputs.avgFluidTemp = inputs.inletTemp; + return outputs; + } + + Real64 const rB = calcBHAverageResistance(inputs.inletTemp, massFlowRatePerBorehole); + Real64 const r12 = calcDirectCouplingResistance(inputs.inletTemp, massFlowRatePerBorehole, rB); + + SegmentInputs segIn; + segIn.flowRate = massFlowRatePerBorehole; + segIn.boundaryTemp = bhWallTemp; + segIn.bhResistance = rB; + segIn.dcResistance = r12; + + for (int iter = 0; iter < std::max(1, config.numIterations); ++iter) { + for (std::size_t i = 0; i < segments.size(); ++i) { + segIn.inletTemp1 = (i == 0) ? inputs.inletTemp : segments[i - 1].outlet1(); + segIn.inletTemp2 = (i == segments.size() - 1) ? uBendTemp : segments[i + 1].outlet2(); + segments[i].step(inputs.timeStepSeconds, segIn, fluidProps); + } + uBendTemp = segments.back().outlet1(); + } + + outletTemp = segments.front().outlet2(); + + Real64 boreholeHeatRatePerBorehole = 0.0; + for (auto const &seg : segments) { + boreholeHeatRatePerBorehole += seg.boreholeHeatRate(); + } + + Real64 const cp = fluidProps.cp(inputs.inletTemp); + Real64 const heatRate = inputs.massFlowRate * cp * (inputs.inletTemp - outletTemp); + energy = boreholeHeatRatePerBorehole / clampMin(config.boreholeLength, small) * static_cast(inputs.timeStepSeconds); + + ModelStepOutputs outputs; + outputs.outletTemp = outletTemp; + outputs.heatRate = heatRate; + outputs.boreholeHeatRate = boreholeHeatRatePerBorehole * numBoreholes; + outputs.boreholeWallTemp = bhWallTemp; + outputs.avgFluidTemp = 0.5 * (inputs.inletTemp + outletTemp); + return outputs; + } + + [[nodiscard]] Real64 pipeResistance(Real64 massFlowRate, Real64 temp) const + { + if (massFlowRate <= 0.0) { + return 0.0; + } + + auto frictionFactor = [](Real64 reynoldsNum) { + constexpr Real64 lowerLimit = 1500.0; + constexpr Real64 upperLimit = 5000.0; + if (reynoldsNum < lowerLimit) { + return 64.0 / clampMin(reynoldsNum, 1.0); + } + if (reynoldsNum < upperLimit) { + Real64 const fLow = 64.0 / clampMin(reynoldsNum, 1.0); + Real64 const fHigh = std::pow(0.79 * std::log(reynoldsNum) - 1.64, -2.0); + Real64 const sf = 1.0 / (1.0 + std::exp(-(reynoldsNum - 3000.0) / 450.0)); + return (1.0 - sf) * fLow + sf * fHigh; + } + return std::pow(0.79 * std::log(reynoldsNum) - 1.64, -2.0); + }; + + auto turbulentNusselt = [&](Real64 reynoldsNum, Real64 prandtlNum) { + Real64 const f = frictionFactor(reynoldsNum); + return (f / 8.0) * (reynoldsNum - 1000.0) * prandtlNum / (1.0 + 12.7 * std::sqrt(f / 8.0) * (std::pow(prandtlNum, 2.0 / 3.0) - 1.0)); + }; + + Real64 const mu = clampMin(fluidProps.viscosity(temp), small); + Real64 const cp = clampMin(fluidProps.cp(temp), small); + Real64 const k = clampMin(fluidProps.conductivity(temp), small); + Real64 const di = clampMin(config.pipeInnerDiameter, small); + Real64 const doPipe = std::max(config.pipeOuterDiameter, di + small); + Real64 const re = 4.0 * massFlowRate / (mu * Constant::Pi * di); + Real64 const pr = cp * mu / k; + + Real64 nu = 0.0; + if (re < 2000.0) { + nu = 4.01; + } else if (re < 4000.0) { + Real64 const nuLow = 4.01; + Real64 const nuHigh = turbulentNusselt(re, pr); + Real64 const sf = 1.0 / (1.0 + std::exp(-(re - 3000.0) / 150.0)); + nu = (1.0 - sf) * nuLow + sf * nuHigh; + } else { + nu = turbulentNusselt(re, pr); + } + nu = clampMin(nu, small); + + Real64 const conv = 1.0 / (nu * Constant::Pi * k); + Real64 const cond = std::log(doPipe / di) / (2.0 * Constant::Pi * clampMin(config.pipeConductivity, small)); + return conv + cond; + } + + [[nodiscard]] Real64 calcBHAverageResistance(Real64 temp, Real64 flowRate) const + { + Real64 const beta = std::max(0.0, std::min(0.999, 2.0 * Constant::Pi * config.groutConductivity * pipeResistance(flowRate, temp))); + Real64 const t14 = std::pow(theta1, 4); + Real64 const d = clampMin(1.0 - t14, small); + Real64 const term1 = std::log(theta2 / clampMin(2.0 * theta1 * std::pow(d, sigma), small)); + Real64 const num2 = std::pow(theta3, 2) * std::pow(1.0 - (4.0 * sigma * t14) / d, 2); + Real64 const den2 = + clampMin((1.0 + beta) / clampMin(1.0 - beta, small) + std::pow(theta3, 2) * (1.0 + (16.0 * sigma * t14) / std::pow(d, 2)), small); + return (1.0 / (4.0 * Constant::Pi * config.groutConductivity)) * (beta + term1 - num2 / den2); + } + + [[nodiscard]] Real64 calcDirectCouplingResistance(Real64 temp, Real64 flowRate, Real64 rB) const + { + Real64 const beta = std::max(0.0, std::min(0.999, 2.0 * Constant::Pi * config.groutConductivity * pipeResistance(flowRate, temp))); + Real64 const t12 = std::pow(theta1, 2); + Real64 const t14 = std::pow(theta1, 4); + Real64 const num2 = std::pow(theta3, 2) * std::pow(1.0 - t14 + 4.0 * sigma * t12, 2); + Real64 const den2 = clampMin(((1.0 + beta) / clampMin(1.0 - beta, small)) * std::pow(1.0 - t14, 2) - + std::pow(theta3, 2) * std::pow(1.0 - t14, 2) + 8.0 * sigma * t12 * std::pow(theta3, 2) * (1.0 + t14), + small); + Real64 const rA = (1.0 / (Constant::Pi * config.groutConductivity)) * + (beta + std::log(std::pow(1.0 + t12, sigma) / clampMin(theta3 * std::pow(1.0 - t12, sigma), small)) - num2 / den2); + Real64 const denom = 4.0 * rB - rA; + if (std::abs(denom) < small) { + return 70.0; + } + Real64 const r12 = (4.0 * rA * rB) / denom; + return (r12 > small && std::isfinite(r12)) ? r12 : 70.0; + } + + ModelConfig config; + FluidPropertyFunctions fluidProps; + DynamicAggregation aggregation; + std::vector segments; + Real64 theta1 = 0.0; + Real64 theta2 = 0.0; + Real64 theta3 = 0.0; + Real64 sigma = 0.0; + Real64 ts = 0.0; + Real64 c0 = 0.0; + Real64 energy = 0.0; + Real64 uBendTemp = 20.0; + Real64 outletTemp = 20.0; + Real64 bhWallTemp = 20.0; + int lastSimTime = -1; + bool haveStepStartSnapshot = false; + DynamicAggregation aggregationAtStepStart; + std::vector segmentsAtStepStart; + Real64 energyAtStepStart = 0.0; + Real64 uBendTempAtStepStart = 20.0; + Real64 outletTempAtStepStart = 20.0; + Real64 bhWallTempAtStepStart = 20.0; +}; + +Model::Model(ModelConfig config, FluidPropertyFunctions fluidProps) : impl(std::make_unique(std::move(config), std::move(fluidProps))) +{ +} + +Model::~Model() = default; + +void Model::reset(Real64 initialTemp) +{ + impl->reset(initialTemp); +} + +ModelStepOutputs Model::simulate(ModelStepInputs const &inputs) +{ + return impl->simulate(inputs); +} + +} // namespace EnergyPlus::GroundHeatExchangers::GLHEC diff --git a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.hh b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.hh new file mode 100644 index 00000000000..a2707ce2602 --- /dev/null +++ b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.hh @@ -0,0 +1,145 @@ +// EnergyPlus, Copyright (c) 1996-present, The Board of Trustees of the University of Illinois, +// The Regents of the University of California, through Lawrence Berkeley National Laboratory +// (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge +// National Laboratory, managed by UT-Battelle, Alliance for Energy Innovation, LLC, and other +// contributors. All rights reserved. +// +// NOTICE: This Software was developed under funding from the U.S. Department of Energy and the +// U.S. Government consequently retains certain rights. As such, the U.S. Government has been +// granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable, +// worldwide license in the Software to reproduce, distribute copies to the public, prepare +// derivative works, and perform publicly and display publicly, and to permit others to do so. +// +// Redistribution and use in source and binary forms, with or without modification, are permitted +// provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this list of +// conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, this list of +// conditions and the following disclaimer in the documentation and/or other materials +// provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, +// the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without specific prior +// written permission. +// +// (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form +// without changes from the version obtained under this License, or (ii) Licensee makes a +// reference solely to the software portion of its product, Licensee must refer to the +// software as "EnergyPlus version X" software, where "X" is the version number Licensee +// obtained under this License and may not use a different name for the software. Except as +// specifically required in this Section (4), Licensee shall not use in a company name, a +// product name, in advertising, publicity, or other promotional activities any name, trade +// name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly +// similar designation, without the U.S. Department of Energy's prior written consent. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY +// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +#pragma once + +#include +#include +#include + +#include + +namespace EnergyPlus::GroundHeatExchangers::GLHEC { + +struct FluidPropertyFunctions +{ + std::function cp; + std::function rho; + std::function viscosity; + std::function conductivity; +}; + +struct DynamicAggregationSettings +{ + Real64 expansionRate = 1.62; + int binsPerLevel = 9; + Real64 simulationHorizonSeconds = 50.0 * 365.0 * 24.0 * 3600.0; +}; + +struct OdeSolverSettings +{ + Real64 absoluteTolerance = 1.0e-10; + Real64 relativeTolerance = 1.0e-8; + int stepsPerTimeStep = 20; + Real64 minInitialStep = 1.0; +}; + +struct ModelConfig +{ + Real64 boreholeLength = 0.0; + Real64 boreholeDiameter = 0.0; + Real64 shankSpacing = 0.0; + + Real64 groutConductivity = 0.0; + Real64 groutDensity = 0.0; + Real64 groutSpecificHeat = 0.0; + + Real64 pipeConductivity = 0.0; + Real64 pipeDensity = 0.0; + Real64 pipeSpecificHeat = 0.0; + Real64 pipeInnerDiameter = 0.0; + Real64 pipeOuterDiameter = 0.0; + + Real64 soilConductivity = 0.0; + Real64 soilDiffusivity = 0.0; + + unsigned int numBoreholes = 1; + int numSegments = 10; + int numIterations = 2; + Real64 groutFraction = 0.5; + int pipeTransitCells = 16; + bool applyPipeTransitDelay = true; + + DynamicAggregationSettings aggregation{}; + OdeSolverSettings ode{}; + std::vector lntts; + std::vector gValues; + std::vector gBValues; +}; + +struct ModelStepInputs +{ + int timeSeconds = 0; + int timeStepSeconds = 0; + Real64 massFlowRate = 0.0; + Real64 inletTemp = 0.0; + Real64 farFieldGroundTemp = 0.0; +}; + +struct ModelStepOutputs +{ + Real64 outletTemp = 0.0; + Real64 heatRate = 0.0; + Real64 boreholeHeatRate = 0.0; + Real64 boreholeWallTemp = 0.0; + Real64 avgFluidTemp = 0.0; +}; + +class Model +{ +public: + Model(ModelConfig config, FluidPropertyFunctions fluidProps); + ~Model(); + + void reset(Real64 initialTemp); + [[nodiscard]] ModelStepOutputs simulate(ModelStepInputs const &inputs); + +private: + struct Impl; + std::unique_ptr impl; +}; + +} // namespace EnergyPlus::GroundHeatExchangers::GLHEC diff --git a/src/EnergyPlus/GroundHeatExchangers/README.md b/src/EnergyPlus/GroundHeatExchangers/README.md index 56b6857dba6..0a7446a0c7a 100644 --- a/src/EnergyPlus/GroundHeatExchangers/README.md +++ b/src/EnergyPlus/GroundHeatExchangers/README.md @@ -23,6 +23,11 @@ The data defining this function is read from input. The heat pulse histories need to be recorded over an extended period (months). To aid computational efficiency past pulses are continuously aggregated into equivalent heat pulses of longer duration, as each pulse becomes less recent. +## GLHEC Mode +`GroundHeatExchanger:System` also supports `g-Function Calculation Method = GLHECcalc` to run the GLHEC transient borehole model located under `GroundHeatExchangers/GLHEC/`. + +This mode uses in-tree RK4 ODE integration for segment state advancement and consumes the same long-term g-function curves used by the vertical GLHE framework. + ## References - Eskilson, P. 'Thermal Analysis of Heat Extraction Boreholes' Ph.D. Thesis: Dept. of Mathematical Physics, University of Lund, Sweden, June 1987. - Yavuzturk, C., J.D. Spitler. 1999. 'A Short Time Step Response Factor Model for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2): 475-485. diff --git a/src/EnergyPlus/GroundHeatExchangers/Vertical.cc b/src/EnergyPlus/GroundHeatExchangers/Vertical.cc index a7ea9305006..5b9a6df7a52 100644 --- a/src/EnergyPlus/GroundHeatExchangers/Vertical.cc +++ b/src/EnergyPlus/GroundHeatExchangers/Vertical.cc @@ -61,6 +61,8 @@ #include #include +#include + namespace EnergyPlus::GroundHeatExchangers { GLHEVert::GLHEVert(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j) { @@ -111,6 +113,40 @@ GLHEVert::GLHEVert(EnergyPlusData &state, std::string const &objName, nlohmann:: this->soil.k = j["ground_thermal_conductivity"].get(); this->soil.rhoCp = j["ground_thermal_heat_capacity"].get(); + this->gFuncCalcMethod = GFuncCalcMethod::UniformHeatFlux; + if (j.find("g_function_calculation_method") != j.end()) { + std::string const gFunctionMethodStr = Util::makeUPPER(j["g_function_calculation_method"].get()); + if (gFunctionMethodStr == "UHFCALC") { + this->gFuncCalcMethod = GFuncCalcMethod::UniformHeatFlux; + } else if (gFunctionMethodStr == "UBHWTCALC") { + this->gFuncCalcMethod = GFuncCalcMethod::UniformBoreholeWallTemp; + } else if (gFunctionMethodStr == "FULLDESIGN") { + this->gFuncCalcMethod = GFuncCalcMethod::FullDesign; + } else if (gFunctionMethodStr == "GLHECCALC") { + this->gFuncCalcMethod = GFuncCalcMethod::UniformHeatFlux; + this->useGLHEC = true; + } else { + errorsFound = true; + ShowSevereError(state, fmt::format("g-Function Calculation Method: \"{}\" is invalid", gFunctionMethodStr)); + } + } + + if (auto const it = j.find("glhec_number_of_borehole_segments"); it != j.end()) { + this->glhecNumSegments = std::max(1, it.value().get()); + } + if (auto const it = j.find("glhec_number_of_iterations"); it != j.end()) { + this->glhecNumIterations = std::max(1, it.value().get()); + } + if (auto const it = j.find("glhec_grout_fraction"); it != j.end()) { + this->glhecGroutFraction = it.value().get(); + } + if (auto const it = j.find("glhec_aggregation_expansion_rate"); it != j.end()) { + this->glhecAggExpansionRate = it.value().get(); + } + if (auto const it = j.find("glhec_aggregation_bins_per_level"); it != j.end()) { + this->glhecAggBinsPerLevel = std::max(1, it.value().get()); + } + if (j.find("ghe_vertical_responsefactors_object_name") != j.end()) { // Response factors come from IDF object this->myRespFactors = GetResponseFactor(state, Util::makeUPPER(j["ghe_vertical_responsefactors_object_name"].get())); @@ -120,21 +156,6 @@ GLHEVert::GLHEVert(EnergyPlusData &state, std::string const &objName, nlohmann:: // no g-functions in the input file, so they need to be calculated if (!this->gFunctionsExist) { - // g-function calculation method - if (j.find("g_function_calculation_method") != j.end()) { - std::string gFunctionMethodStr = Util::makeUPPER(j["g_function_calculation_method"].get()); - if (gFunctionMethodStr == "UHFCALC") { - this->gFuncCalcMethod = GFuncCalcMethod::UniformHeatFlux; - } else if (gFunctionMethodStr == "UBHWTCALC") { - this->gFuncCalcMethod = GFuncCalcMethod::UniformBoreholeWallTemp; - } else if (gFunctionMethodStr == "FULLDESIGN") { - this->gFuncCalcMethod = GFuncCalcMethod::FullDesign; - } else { - errorsFound = true; - ShowSevereError(state, fmt::format("g-Function Calculation Method: \"{}\" is invalid", gFunctionMethodStr)); - } - } - // get borehole data from array or individual borehole instance objects if (this->gFuncCalcMethod == GFuncCalcMethod::FullDesign) { #ifndef PYTHON_CLI @@ -407,6 +428,42 @@ void GLHEVert::simulate(EnergyPlusData &state, return; } + if (this->useGLHEC) { + if (!this->gFunctionsExist) { + this->calcGFunctions(state); + this->gFunctionsExist = true; + } + if (!this->glhecModel) { + this->setupGLHECModel(state); + } + + this->inletTemp = state.dataLoopNodes->Node(this->inletNodeNum).Temp; + + Real64 const currTimeSeconds = ((state.dataGlobal->DayOfSim - 1) * 24 + (state.dataGlobal->HourOfDay - 1) + + (state.dataGlobal->TimeStep - 1) * state.dataGlobal->TimeStepZone + state.dataHVACGlobal->SysTimeElapsed) * + Constant::rSecsInHour; + + GLHEC::ModelStepInputs stepInputs; + stepInputs.timeSeconds = std::max(0, static_cast(std::llround(currTimeSeconds))); + stepInputs.timeStepSeconds = std::max(1, static_cast(std::llround(state.dataHVACGlobal->TimeStepSysSec))); + stepInputs.massFlowRate = this->massFlowRate; + stepInputs.inletTemp = this->inletTemp; + stepInputs.farFieldGroundTemp = this->tempGround; + + auto const outputs = this->glhecModel->simulate(stepInputs); + this->outletTemp = outputs.outletTemp; + this->QGLHE = outputs.heatRate; + this->bhTemp = outputs.boreholeWallTemp; + this->aveFluidTemp = outputs.avgFluidTemp; + if (this->totalTubeLength > 0.0) { + this->lastQnSubHr = outputs.boreholeHeatRate / this->totalTubeLength; + } else { + this->lastQnSubHr = 0.0; + } + this->updateGHX(state); + return; + } + if (this->gFuncCalcMethod == GFuncCalcMethod::FullDesign) { // we need to do some special things for the full design mode this->outletTemp = this->tempGround; @@ -1499,6 +1556,64 @@ void GLHEVert::initGLHESimVars(EnergyPlusData &state) } } +void GLHEVert::setupGLHECModel(EnergyPlusData &state) +{ + GLHEC::ModelConfig cfg; + cfg.boreholeLength = this->bhLength; + cfg.boreholeDiameter = this->bhDiameter; + cfg.shankSpacing = this->bhUTubeDist; + cfg.groutConductivity = this->grout.k; + cfg.soilConductivity = this->soil.k; + cfg.soilDiffusivity = this->soil.diffusivity; + cfg.pipeConductivity = this->pipe.k; + cfg.pipeInnerDiameter = this->pipe.innerDia; + cfg.pipeOuterDiameter = this->pipe.outDia; + cfg.numBoreholes = this->myRespFactors->numBoreholes; + cfg.numSegments = this->glhecNumSegments; + cfg.numIterations = this->glhecNumIterations; + cfg.groutFraction = this->glhecGroutFraction; + cfg.pipeTransitCells = 16; + cfg.applyPipeTransitDelay = true; + cfg.ode.absoluteTolerance = 1.0e-10; + cfg.ode.relativeTolerance = 1.0e-8; + cfg.ode.stepsPerTimeStep = 20; + cfg.ode.minInitialStep = 1.0; + + // Input objects provide volumetric heat capacity. These default density assumptions preserve rho*cp. + constexpr Real64 assumedGroutDensity = 2200.0; + constexpr Real64 assumedPipeDensity = 950.0; + cfg.groutDensity = assumedGroutDensity; + cfg.groutSpecificHeat = this->grout.rhoCp / assumedGroutDensity; + cfg.pipeDensity = assumedPipeDensity; + cfg.pipeSpecificHeat = this->pipe.rhoCp / assumedPipeDensity; + + cfg.aggregation.expansionRate = this->glhecAggExpansionRate; + cfg.aggregation.binsPerLevel = this->glhecAggBinsPerLevel; + cfg.aggregation.simulationHorizonSeconds = this->myRespFactors->maxSimYears * 365.0 * 24.0 * 3600.0; + + cfg.lntts = this->myRespFactors->LNTTS; + cfg.gValues = this->myRespFactors->GFNC; + cfg.gBValues = this->myRespFactors->GFNC; + + EnergyPlusData *statePtr = &state; + GLHEC::FluidPropertyFunctions fluidFuncs; + fluidFuncs.cp = [this, statePtr](Real64 temp) { + return statePtr->dataPlnt->PlantLoop(this->plantLoc.loopNum).glycol->getSpecificHeat(*statePtr, temp, "GLHEC::cp"); + }; + fluidFuncs.rho = [this, statePtr](Real64 temp) { + return statePtr->dataPlnt->PlantLoop(this->plantLoc.loopNum).glycol->getDensity(*statePtr, temp, "GLHEC::rho"); + }; + fluidFuncs.viscosity = [this, statePtr](Real64 temp) { + return statePtr->dataPlnt->PlantLoop(this->plantLoc.loopNum).glycol->getViscosity(*statePtr, temp, "GLHEC::mu"); + }; + fluidFuncs.conductivity = [this, statePtr](Real64 temp) { + return statePtr->dataPlnt->PlantLoop(this->plantLoc.loopNum).glycol->getConductivity(*statePtr, temp, "GLHEC::k"); + }; + + this->glhecModel = std::make_unique(std::move(cfg), std::move(fluidFuncs)); + this->glhecModel->reset(this->tempGround); +} + void GLHEVert::initEnvironment(EnergyPlusData &state, [[maybe_unused]] Real64 const CurTime) { constexpr std::string_view RoutineName = "initEnvironment"; @@ -1521,6 +1636,9 @@ void GLHEVert::initEnvironment(EnergyPlusData &state, [[maybe_unused]] Real64 co this->currentSimTime = 0.0; this->QGLHE = 0.0; this->prevHour = 1; + if (this->glhecModel) { + this->glhecModel->reset(this->tempGround); + } } void GLHEVert::oneTimeInit_new(EnergyPlusData &state) diff --git a/src/EnergyPlus/GroundHeatExchangers/Vertical.hh b/src/EnergyPlus/GroundHeatExchangers/Vertical.hh index d38ceb2b54c..1c2a5fe1931 100644 --- a/src/EnergyPlus/GroundHeatExchangers/Vertical.hh +++ b/src/EnergyPlus/GroundHeatExchangers/Vertical.hh @@ -47,7 +47,10 @@ #pragma once +#include + #include +#include namespace EnergyPlus { @@ -81,6 +84,13 @@ namespace GroundHeatExchangers { Real64 bhLength = 0.0; // Length of borehole {m} Real64 bhUTubeDist = 0.0; // Distance between u-tube legs {m} GFuncCalcMethod gFuncCalcMethod = GFuncCalcMethod::Invalid; + bool useGLHEC = false; + int glhecNumSegments = 10; + int glhecNumIterations = 2; + Real64 glhecGroutFraction = 0.5; + Real64 glhecAggExpansionRate = 1.62; + int glhecAggBinsPerLevel = 9; + std::unique_ptr glhecModel; // Parameters for the multipole method Real64 theta_1 = 0.0; @@ -96,6 +106,10 @@ namespace GroundHeatExchangers { GLHEVert() = default; GLHEVert(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j); + GLHEVert(GLHEVert const &) = delete; + GLHEVert &operator=(GLHEVert const &) = delete; + GLHEVert(GLHEVert &&) noexcept = default; + GLHEVert &operator=(GLHEVert &&) noexcept = default; ~GLHEVert() override = default; void simulate([[maybe_unused]] EnergyPlusData &state, @@ -161,6 +175,8 @@ namespace GroundHeatExchangers { void oneTimeInit_new(EnergyPlusData &state) override; void setupTimeVectors() const; + + void setupGLHECModel(EnergyPlusData &state); }; } // namespace GroundHeatExchangers diff --git a/tst/EnergyPlus/unit/CMakeLists.txt b/tst/EnergyPlus/unit/CMakeLists.txt index 75c1fd20184..464fea0b971 100644 --- a/tst/EnergyPlus/unit/CMakeLists.txt +++ b/tst/EnergyPlus/unit/CMakeLists.txt @@ -132,6 +132,8 @@ set(test_src GroundHeatExchangers/Base.unit.cc GroundHeatExchangers/BoreholeArray.unit.cc GroundHeatExchangers/BoreholeSingle.unit.cc + GroundHeatExchangers/GLHEC.unit.cc + GroundHeatExchangers/GLHECModel.unit.cc GroundHeatExchangers/Properties.unit.cc GroundHeatExchangers/ResponseFactors.unit.cc GroundHeatExchangers/Slinky.unit.cc diff --git a/tst/EnergyPlus/unit/GroundHeatExchangers/GLHEC.unit.cc b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHEC.unit.cc new file mode 100644 index 00000000000..9cba2ca56bf --- /dev/null +++ b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHEC.unit.cc @@ -0,0 +1,110 @@ +// EnergyPlus, Copyright (c) 1996-present, The Board of Trustees of the University of Illinois, +// The Regents of the University of California, through Lawrence Berkeley National Laboratory +// (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge +// National Laboratory, managed by UT-Battelle, Alliance for Energy Innovation, LLC, and other +// contributors. All rights reserved. +// +// NOTICE: This Software was developed under funding from the U.S. Department of Energy and the +// U.S. Government consequently retains certain rights. As such, the U.S. Government has been +// granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable, +// worldwide license in the Software to reproduce, distribute copies to the public, prepare +// derivative works, and perform publicly and display publicly, and to permit others to do so. +// +// Redistribution and use in source and binary forms, with or without modification, are permitted +// provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this list of +// conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, this list of +// conditions and the following disclaimer in the documentation and/or other materials +// provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, +// the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without specific prior +// written permission. +// +// (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form +// without changes from the version obtained under this License, or (ii) Licensee makes a +// reference solely to the software portion of its product, Licensee must refer to the +// software as "EnergyPlus version X" software, where "X" is the version number Licensee +// obtained under this License and may not use a different name for the software. Except as +// specifically required in this Section (4), Licensee shall not use in a company name, a +// product name, in advertising, publicity, or other promotional activities any name, trade +// name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly +// similar designation, without the U.S. Department of Energy's prior written consent. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY +// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +#include "../Fixtures/EnergyPlusFixture.hh" + +#include +#include +#include +#include + +using namespace EnergyPlus; +using namespace EnergyPlus::GroundHeatExchangers; + +TEST_F(EnergyPlusFixture, GroundHeatExchanger_GLHEC_Mode_Input_Parses) +{ + std::string const idf_objects = delimited_string({ + "Site:GroundTemperature:Undisturbed:KusudaAchenbach,", + " KATemps, !- Name", + " 1.8, !- Soil Thermal Conductivity {W/m-K}", + " 920, !- Soil Density {kg/m3}", + " 2200, !- Soil Specific Heat {J/kg-K}", + " 15.5, !- Average Soil Surface Temperature {C}", + " 3.2, !- Average Amplitude of Surface Temperature {deltaC}", + " 8; !- Phase Shift of Minimum Surface Temperature {days}", + + "GroundHeatExchanger:Vertical:Properties,", + " GHE-1 Props, !- Name", + " 1, !- Depth of Top of Borehole {m}", + " 100, !- Borehole Length {m}", + " 0.109982, !- Borehole Diameter {m}", + " 0.744, !- Grout Thermal Conductivity {W/m-K}", + " 3.90E+06, !- Grout Thermal Heat Capacity {J/m3-K}", + " 0.389, !- Pipe Thermal Conductivity {W/m-K}", + " 1.77E+06, !- Pipe Thermal Heat Capacity {J/m3-K}", + " 0.0267, !- Pipe Outer Diameter {m}", + " 0.00243, !- Pipe Thickness {m}", + " 0.04556; !- U-Tube Distance {m}", + + "GroundHeatExchanger:Vertical:Single,", + " GHE-Single-1, !- Name", + " GHE-1 Props, !- GHE:Vertical:Properties Object Name", + " 0.0, !- X-Location", + " 0.0; !- Y-Location", + + "GroundHeatExchanger:System,", + " GLHEC Test, !- Name", + " GHLE Inlet, !- Inlet Node Name", + " GHLE Outlet, !- Outlet Node Name", + " 0.0007571, !- Design Flow Rate {m3/s}", + " Site:GroundTemperature:Undisturbed:KusudaAchenbach, !- Undisturbed Ground Temperature Model Type", + " KATemps, !- Undisturbed Ground Temperature Model Name", + " 2.423, !- Ground Thermal Conductivity {W/m-K}", + " 2.343E+06, !- Ground Thermal Heat Capacity {J/m3-K}", + " , !- GHE:Vertical:ResponseFactors Object Name", + " GLHECcalc, !- g-Function Calculation Method", + " , !- GHE:Vertical:Sizing Object Type", + " , !- GHE:Vertical:Sizing Object Name", + " , !- GHE:Vertical:Array Object Name", + " GHE-Single-1; !- GHE:Vertical:Single Object Name 1", + }); + + ASSERT_TRUE(process_idf(idf_objects)); + GetGroundHeatExchangerInput(*state); + ASSERT_EQ(1u, state->dataGroundHeatExchanger->verticalGLHE.size()); + auto &thisGLHE(state->dataGroundHeatExchanger->verticalGLHE[0]); + EXPECT_TRUE(thisGLHE.useGLHEC); +} diff --git a/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc new file mode 100644 index 00000000000..05f4168b306 --- /dev/null +++ b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc @@ -0,0 +1,450 @@ +// EnergyPlus, Copyright (c) 1996-present, The Board of Trustees of the University of Illinois, +// The Regents of the University of California, through Lawrence Berkeley National Laboratory +// (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge +// National Laboratory, managed by UT-Battelle, Alliance for Energy Innovation, LLC, and other +// contributors. All rights reserved. +// +// NOTICE: This Software was developed under funding from the U.S. Department of Energy and the +// U.S. Government consequently retains certain rights. As such, the U.S. Government has been +// granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable, +// worldwide license in the Software to reproduce, distribute copies to the public, prepare +// derivative works, and perform publicly and display publicly, and to permit others to do so. +// +// Redistribution and use in source and binary forms, with or without modification, are permitted +// provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this list of +// conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, this list of +// conditions and the following disclaimer in the documentation and/or other materials +// provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, +// the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without specific prior +// written permission. +// +// (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form +// without changes from the version obtained under this License, or (ii) Licensee makes a +// reference solely to the software portion of its product, Licensee must refer to the +// software as "EnergyPlus version X" software, where "X" is the version number Licensee +// obtained under this License and may not use a different name for the software. Except as +// specifically required in this Section (4), Licensee shall not use in a company name, a +// product name, in advertising, publicity, or other promotional activities any name, trade +// name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly +// similar designation, without the U.S. Department of Energy's prior written consent. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY +// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +#include + +#include +#include + +#include + +using namespace EnergyPlus; +using namespace EnergyPlus::GroundHeatExchangers::GLHEC; + +namespace { + +ModelConfig makeConfig() +{ + ModelConfig cfg; + cfg.boreholeLength = 100.0; + cfg.boreholeDiameter = 0.109982; + cfg.shankSpacing = 0.04556; + + cfg.groutConductivity = 0.744; + cfg.groutDensity = 2200.0; + cfg.groutSpecificHeat = 1772.7272727273; + + cfg.pipeConductivity = 0.389; + cfg.pipeDensity = 950.0; + cfg.pipeSpecificHeat = 1863.1578947368; + cfg.pipeInnerDiameter = 0.02184; + cfg.pipeOuterDiameter = 0.02670; + + cfg.soilConductivity = 2.423; + cfg.soilDiffusivity = 1.03e-6; + + cfg.numBoreholes = 10; + cfg.numSegments = 8; + cfg.numIterations = 2; + cfg.groutFraction = 0.5; + + cfg.lntts = {-12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0, 2.0}; + cfg.gValues = {0.0, 0.10, 0.30, 0.65, 1.05, 1.45, 1.90, 2.20}; + cfg.gBValues = cfg.gValues; + + return cfg; +} + +FluidPropertyFunctions makeFluidFuncs() +{ + FluidPropertyFunctions funcs; + funcs.cp = [](Real64) { return 4180.0; }; + funcs.rho = [](Real64) { return 997.0; }; + funcs.viscosity = [](Real64) { return 1.0e-3; }; + funcs.conductivity = [](Real64) { return 0.60; }; + return funcs; +} + +ModelConfig makePrototypeLikeConfig() +{ + ModelConfig cfg; + cfg.boreholeLength = 76.2; + cfg.boreholeDiameter = 0.114; + cfg.shankSpacing = 0.0469; + + cfg.groutConductivity = 0.85; + cfg.groutDensity = 2500.0; + cfg.groutSpecificHeat = 1560.0; + + cfg.pipeConductivity = 0.39; + cfg.pipeDensity = 950.0; + cfg.pipeSpecificHeat = 1900.0; + cfg.pipeInnerDiameter = 0.0218; + cfg.pipeOuterDiameter = 0.0267; + + cfg.soilConductivity = 2.7; + cfg.soilDiffusivity = 2.7 / (2500.0 * 880.0); + + cfg.numBoreholes = 1; + cfg.numSegments = 10; + cfg.numIterations = 2; + cfg.groutFraction = 0.5; + + cfg.aggregation.expansionRate = 1.5; + cfg.aggregation.binsPerLevel = 9; + cfg.aggregation.simulationHorizonSeconds = 8760.0 * 3600.0; + + cfg.lntts = {-17.075114469809652, -16.627802251765988, -16.331956868675046, -16.153708637268725, -15.959552622827768, -15.77583148567939, + -15.682430310590989, -15.591181639125844, -15.405957322673929, -15.313726098457895, -15.217215198077051, -15.125548009551228, + -15.034893641283096, -14.939583461478772, -14.847031904113527, -14.754739835252108, -14.661033872590263, -14.569065383447843, + -14.475277853347674, -14.380833849431198, -14.288205693795122, -14.196204155432493, -14.102729560252651, -14.008852045449823, + -13.916478725318807, -13.823945223220736, -13.731075501987444, -13.637824206995248, -13.544239264838964, -13.451115475629575, + -13.358420327512839, -13.2656070279071, -13.172365101166834, -13.079504748273239, -12.986770782860688, -12.893651698697814, + -12.800587389042345, -12.70761593569868, -12.614525512632667, -12.52123757820911, -12.428015210032074, -12.335099651310028, + -12.242093878427184, -12.148976233736313, -12.055968486045487, -11.962972229371559, -11.869981339358326, -11.777053546396756, + -11.684054318962536, -11.590956105944782, -11.497903536842914, -11.404879940036198, -11.311854335522456, -11.218794337158441, + -11.125741419279972, -11.032725286003604, -10.939716092757104, -10.846654041512448, -10.753608576410747, -10.660605202704337, + -10.567568901988443, -10.47451163808514, -10.381473119346294, -10.28845540129876, -10.19541578861004, -10.102388070966573, + -10.009347841117759, -9.9163012056817497, -9.8232786085743538, -9.73024635104019, -9.6372137108808698, -9.5441806323696312, + -9.4511467373365736, -9.3581170082252108, -9.2650841704847497, -9.172048759581525, -9.0790130948574852, -8.9859781305394382, + -8.8929422609011795, -8.7999079175091044, -8.7068735276925988, -8.6138406411697463, -8.5208083076518033, -8.427775096313507, + -8.3347408046839213, -8.2416330604701553}; + + cfg.gValues = {0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.028697692672201947, + 0.16339564189400591, + 0.39327037984064112, + 0.50389456883620731, + 0.54695324690661495, + 0.58875261192606554, + 0.61783257574788653, + 0.65801278168159738, + 0.78873155164181841, + 0.88228609685010839, + 0.93910764585584006, + 1.0012676082691361, + 1.1205677647618273, + 1.1987460294171912, + 1.2574993816078572, + 1.3378092865578339, + 1.4281407217225641, + 1.5151395142688655, + 1.5852001726168479, + 1.6613621044219196, + 1.7467449063896214, + 1.8480244254278784, + 1.9209020018623595, + 1.9902883975365804, + 2.084705199593365, + 2.1550243702289102, + 2.2014434008991048, + 2.289503427220219, + 2.3590624939083842, + 2.4375690842178557, + 2.5160768546333778, + 2.5924863543159158, + 2.6579691104092289, + 2.7150337596936618, + 2.7972212977700956, + 2.8467817516534972, + 2.9329572415227858, + 3.008453856740366, + 3.106681865004925, + 3.1490749450787567, + 3.191614994861546, + 3.2520549577221134, + 3.3094300605348148, + 3.3735590368801427, + 3.455360812349086, + 3.4946273633674028, + 3.5586651615672968, + 3.6160827173698085, + 3.6667717710707315, + 3.710886403325091, + 3.7575826652378361, + 3.8116730387671018, + 3.8477486812218298, + 3.9010544628374566, + 3.9570877386851762, + 4.0000787914203375, + 4.0458294516432316, + 4.1028019631225465, + 4.1501005489265745, + 4.1948909786978765, + 4.2261122235003485, + 4.2748815807902218, + 4.3197052900799413, + 4.3675379331172488, + 4.4075060875664631, + 4.4668499186316284, + 4.5120700271156817}; + cfg.gBValues = cfg.gValues; + return cfg; +} + +FluidPropertyFunctions makePrototypeLikeFluidFuncs() +{ + auto waterCp = [](Real64 t) { + constexpr Real64 acp0 = 4.21534; + constexpr Real64 acp1 = -0.00287819; + constexpr Real64 acp2 = 7.4729e-05; + constexpr Real64 acp3 = -7.79624e-07; + constexpr Real64 acp4 = 3.220424e-09; + return (acp0 + t * acp1 + std::pow(t, 2) * acp2 + std::pow(t, 3) * acp3 + std::pow(t, 4) * acp4) * 1000.0; + }; + + auto waterMu = [](Real64 t) { + if (t < 20.0) { + constexpr Real64 am0 = -3.30233; + constexpr Real64 am1 = 1301.0; + constexpr Real64 am2 = 998.333; + constexpr Real64 am3 = 8.1855; + constexpr Real64 am4 = 0.00585; + Real64 const exponent = am0 + am1 / (am2 + (t - 20.0) * (am3 + am4 * (t - 20.0))); + return std::pow(10.0, exponent) * 0.1; + } + if (t > 100.0) { + constexpr Real64 am10 = 0.68714; + constexpr Real64 am11 = -0.0059231; + constexpr Real64 am12 = 2.1249e-05; + constexpr Real64 am13 = -2.69575e-08; + return (am10 + t * am11 + std::pow(t, 2) * am12 + std::pow(t, 3) * am13) * 0.001; + } + constexpr Real64 am5 = 1.002; + constexpr Real64 am6 = -1.3272; + constexpr Real64 am7 = -0.001053; + constexpr Real64 am8 = 105.0; + Real64 const exponent = (t - 20.0) * (am6 + (t - 20.0) * am7) / (t + am8); + return am5 * std::pow(10.0, exponent) * 0.001; + }; + + auto waterRho = [](Real64 t) { + constexpr Real64 ar0 = 999.83952; + constexpr Real64 ar1 = 16.945176; + constexpr Real64 ar2 = -0.0079870401; + constexpr Real64 ar3 = -4.6170461e-05; + constexpr Real64 ar4 = 1.0556302e-07; + constexpr Real64 ar5 = -2.8054253e-10; + constexpr Real64 ar6 = 0.01687985; + return (ar0 + t * ar1 + std::pow(t, 2) * ar2 + std::pow(t, 3) * ar3 + std::pow(t, 4) * ar4 + std::pow(t, 5) * ar5) / (1.0 + ar6 * t); + }; + + auto waterK = [](Real64 t) { + constexpr Real64 ak0 = 0.560101; + constexpr Real64 ak1 = 0.00211703; + constexpr Real64 ak2 = -1.05172e-05; + constexpr Real64 ak3 = 1.497323e-08; + constexpr Real64 ak4 = -1.48553e-11; + return ak0 + t * ak1 + std::pow(t, 2) * ak2 + std::pow(t, 3) * ak3 + std::pow(t, 4) * ak4; + }; + + FluidPropertyFunctions funcs; + funcs.cp = std::move(waterCp); + funcs.rho = std::move(waterRho); + funcs.viscosity = std::move(waterMu); + funcs.conductivity = std::move(waterK); + return funcs; +} + +} // namespace + +TEST(GLHECModel, ZeroFlowMaintainsInletTemperature) +{ + Model model(makeConfig(), makeFluidFuncs()); + model.reset(15.0); + + ModelStepInputs inputs; + inputs.timeSeconds = 600; + inputs.timeStepSeconds = 600; + inputs.massFlowRate = 0.0; + inputs.inletTemp = 18.0; + inputs.farFieldGroundTemp = 15.0; + + auto const outputs = model.simulate(inputs); + + EXPECT_NEAR(outputs.outletTemp, inputs.inletTemp, 1.0e-10); + EXPECT_NEAR(outputs.heatRate, 0.0, 1.0e-10); + EXPECT_NEAR(outputs.boreholeHeatRate, 0.0, 1.0e-10); +} + +TEST(GLHECModel, RepeatedCallsAtSameTimeDoNotAccumulate) +{ + Model model(makeConfig(), makeFluidFuncs()); + model.reset(15.0); + + ModelStepInputs inputs; + inputs.timeSeconds = 600; + inputs.timeStepSeconds = 600; + inputs.massFlowRate = 1.0; + inputs.inletTemp = 20.0; + inputs.farFieldGroundTemp = 15.0; + + auto const first = model.simulate(inputs); + auto const second = model.simulate(inputs); + + EXPECT_NEAR(first.outletTemp, second.outletTemp, 1.0e-9); + EXPECT_NEAR(first.heatRate, second.heatRate, 1.0e-6); + EXPECT_NEAR(first.boreholeHeatRate, second.boreholeHeatRate, 1.0e-6); + EXPECT_NEAR(first.boreholeWallTemp, second.boreholeWallTemp, 1.0e-9); +} + +TEST(GLHECModel, AggregationTuningChangesLongTermResponse) +{ + ModelConfig coarse = makeConfig(); + coarse.aggregation.expansionRate = 2.0; + coarse.aggregation.binsPerLevel = 3; + coarse.aggregation.simulationHorizonSeconds = 20.0 * 24.0 * 3600.0; + + ModelConfig fine = makeConfig(); + fine.aggregation.expansionRate = 1.2; + fine.aggregation.binsPerLevel = 16; + fine.aggregation.simulationHorizonSeconds = 20.0 * 24.0 * 3600.0; + + Model coarseModel(coarse, makeFluidFuncs()); + Model fineModel(fine, makeFluidFuncs()); + coarseModel.reset(15.0); + fineModel.reset(15.0); + + ModelStepInputs inputs; + inputs.timeStepSeconds = 600; + inputs.massFlowRate = 1.0; + inputs.inletTemp = 20.0; + inputs.farFieldGroundTemp = 15.0; + + ModelStepOutputs coarseOut; + ModelStepOutputs fineOut; + for (int i = 1; i <= 300; ++i) { + inputs.timeSeconds = i * inputs.timeStepSeconds; + coarseOut = coarseModel.simulate(inputs); + fineOut = fineModel.simulate(inputs); + } + + EXPECT_GT(std::abs(coarseOut.boreholeWallTemp - fineOut.boreholeWallTemp), 1.0e-6); +} + +TEST(GLHECModel, PrototypeTrajectoryRegressionCheckpoints) +{ + struct Checkpoint + { + int step; + Real64 loadInletTemp; + Real64 loadOutletTemp; + Real64 heatRate; + Real64 boreholeHeatRate; + Real64 boreholeWallTemp; + }; + + // Reference checkpoints from prototype GLHEC run (temps2.csv) at 600 s timesteps. + std::vector const checkpoints = { + {0, 16.1, 19.6839, 2995.95, 444.649, 16.1}, + {1, 16.1021, 19.6861, 2983.52, 422.268, 16.5592}, + {2, 16.1191, 19.7031, 2975.57, 541.087, 16.6577}, + {5, 16.2218, 19.8059, 1530.6, 667.578, 17.1253}, + {10, 18.6001, 22.1861, 2910.73, 1064.14, 17.9731}, + {20, 21.5371, 25.1248, 2855.75, 1552.08, 19.3305}, + {50, 26.1979, 29.7875, 2909.74, 2279.93, 22.0499}, + {100, 29.6169, 33.2072, 2965.61, 2727.24, 24.2576}, + {200, 31.7675, 35.3581, 2990.76, 2926.44, 25.8696}, + {500, 33.3051, 36.8957, 2997.83, 2981.82, 27.2588}, + {1000, 34.1981, 37.7888, 2999.15, 2991.89, 28.1256}, + {2000, 35.0311, 38.6218, 2999.76, 2996.23, 28.9477}, + {5000, 36.0945, 39.6851, 3000.14, 2998.68, 30.0055}, + {10000, 36.8861, 40.4766, 3000.29, 2999.5, 30.7955}, + {20000, 37.6724, 41.2628, 3000.39, 2999.92, 31.5815}, + {50000, 38.7114, 42.3017, 3000.48, 3000.19, 32.6206}, + {52559, 38.7701, 42.3604, 3000.48, 3000.2, 32.6793}, + }; + + auto const fluid = makePrototypeLikeFluidFuncs(); + Model model(makePrototypeLikeConfig(), fluid); + model.reset(16.1); + + constexpr Real64 flowRate = 0.2; + constexpr Real64 loadRate = 3000.0; + constexpr int dtSeconds = 600; + constexpr Real64 farFieldGroundTemp = 16.1; + + Real64 demandInletTemp = 16.1; + std::size_t nextCheckpoint = 0; + + for (int step = 0; step <= checkpoints.back().step; ++step) { + Real64 const cp = fluid.cp(demandInletTemp); + Real64 const glheInletTemp = demandInletTemp + loadRate / (flowRate * cp); + + ModelStepInputs inputs; + inputs.timeSeconds = step * dtSeconds; + inputs.timeStepSeconds = dtSeconds; + inputs.massFlowRate = flowRate; + inputs.inletTemp = glheInletTemp; + inputs.farFieldGroundTemp = farFieldGroundTemp; + + auto const outputs = model.simulate(inputs); + + if (nextCheckpoint < checkpoints.size() && step == checkpoints[nextCheckpoint].step) { + auto const &cpRef = checkpoints[nextCheckpoint]; + + EXPECT_NEAR(demandInletTemp, cpRef.loadInletTemp, 1.5); + EXPECT_NEAR(glheInletTemp, cpRef.loadOutletTemp, 1.5); + EXPECT_NEAR(outputs.boreholeWallTemp, cpRef.boreholeWallTemp, 0.5); + + if (step >= 500) { + EXPECT_NEAR(outputs.heatRate, cpRef.heatRate, 150.0); + EXPECT_NEAR(outputs.boreholeHeatRate, cpRef.boreholeHeatRate, 150.0); + } + ++nextCheckpoint; + } + + demandInletTemp = outputs.outletTemp; + } +} From 88400caccd0716a401e942b28a671704e524e4d2 Mon Sep 17 00:00:00 2001 From: Alex Swindler Date: Mon, 30 Mar 2026 21:02:46 -0600 Subject: [PATCH 2/3] Fixes GLHECcalc instability for multi-borehole EnergyPlus runs by aligning the equivalent-borehole flow/history treatment with the original GLHEC field model. Also makes the segment RK4 integrator choose stable substeps from local thermal time constants to avoid low-flow numerical blow-up --- .../GroundHeatExchangers/GLHEC/Model.cc | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc index ef3868b2e49..a2e44fbac00 100644 --- a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc +++ b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc @@ -517,7 +517,10 @@ namespace { void rk4Step(Real64 timeStepSeconds, Context const &ctx) { - int const n = std::max(20, std::max(1, odeSettings.stepsPerTimeStep) * 5); + Real64 const maxStableSubstep = stableSubstepSeconds(ctx); + int const minSubsteps = std::max(20, std::max(1, odeSettings.stepsPerTimeStep) * 5); + int const stableSubsteps = static_cast(std::ceil(timeStepSeconds / clampMin(maxStableSubstep, small))); + int const n = std::max(minSubsteps, stableSubsteps); Real64 const h = timeStepSeconds / n; std::array k1{}, k2{}, k3{}, k4{}, y2{}, y3{}, y4{}; for (int i = 0; i < n; ++i) { @@ -540,6 +543,30 @@ namespace { } } + [[nodiscard]] Real64 stableSubstepSeconds(Context const &ctx) const + { + Real64 const rF = 1.0 / clampMin(ctx.inputs.flowRate * ctx.fluidCp, small); + Real64 const rB = clampMin(ctx.inputs.bhResistance, small); + Real64 const r12 = clampMin(ctx.inputs.dcResistance, small); + Real64 const cF = clampMin(ctx.fluidRho * ctx.fluidCp * pipeFluidVolume(), small); + Real64 const cG1 = + clampMin(groutFraction * groutSpecificHeat * groutDensity * groutVolume() + pipeSpecificHeat * pipeDensity * pipeWallVolume(), small); + Real64 const cG2 = clampMin( + ((1.0 - groutFraction) * groutSpecificHeat * groutDensity * groutVolume() + pipeSpecificHeat * pipeDensity * pipeWallVolume()) / 2.0, + small); + + Real64 const fluidConductance = 1.0 / rF + 2.0 * length / r12 + length / rB; + Real64 const coreConductance = 4.0 * length / r12; + Real64 const wallConductance = 2.0 * length / rB; + Real64 const tauFluid = cF / clampMin(fluidConductance, small); + Real64 const tauCore = cG1 / clampMin(coreConductance, small); + Real64 const tauWall = cG2 / clampMin(wallConductance, small); + Real64 const tauMin = std::min({tauFluid, tauCore, tauWall}); + + // Keep the explicit RK4 step well below the fastest segment time constant. + return std::max(0.1, 0.25 * tauMin); + } + [[nodiscard]] Real64 pipeFluidVolume() const { return Constant::Pi * std::pow(pipeInnerDiameter, 2) / 4.0 * length; @@ -660,7 +687,9 @@ struct Model::Impl bhWallTemp = inputs.farFieldGroundTemp + hist.first * c0; Real64 const numBoreholes = static_cast(std::max(1u, config.numBoreholes)); - Real64 const massFlowRatePerBorehole = inputs.massFlowRate / numBoreholes; + // Match the original GLHEC field behavior. The prototype applies the field loop flow to + // each borehole path and then normalizes the borehole-wall load history by borehole count. + Real64 const massFlowRatePerBorehole = inputs.massFlowRate; if (massFlowRatePerBorehole <= small) { energy = 0.0; @@ -700,7 +729,8 @@ struct Model::Impl Real64 const cp = fluidProps.cp(inputs.inletTemp); Real64 const heatRate = inputs.massFlowRate * cp * (inputs.inletTemp - outletTemp); - energy = boreholeHeatRatePerBorehole / clampMin(config.boreholeLength, small) * static_cast(inputs.timeStepSeconds); + Real64 const historyHeatRatePerBorehole = (config.numBoreholes > 1) ? (heatRate / numBoreholes) : boreholeHeatRatePerBorehole; + energy = historyHeatRatePerBorehole / clampMin(config.boreholeLength, small) * static_cast(inputs.timeStepSeconds); ModelStepOutputs outputs; outputs.outletTemp = outletTemp; From 3acbd5d8750b3dd7535b76452e399ac25d5ffbe3 Mon Sep 17 00:00:00 2001 From: Matt Mitchell Date: Tue, 31 Mar 2026 10:45:19 -0600 Subject: [PATCH 3/3] fix GLHEC sign GHE heat rate sign convention, fix flow distribution --- .../GroundHeatExchangers/GLHEC/Model.cc | 20 ++++---- .../GroundHeatExchangers/GLHECModel.unit.cc | 51 ++++++++++++++++++- 2 files changed, 58 insertions(+), 13 deletions(-) diff --git a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc index a2e44fbac00..5c7dfda8071 100644 --- a/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc +++ b/src/EnergyPlus/GroundHeatExchangers/GLHEC/Model.cc @@ -687,9 +687,7 @@ struct Model::Impl bhWallTemp = inputs.farFieldGroundTemp + hist.first * c0; Real64 const numBoreholes = static_cast(std::max(1u, config.numBoreholes)); - // Match the original GLHEC field behavior. The prototype applies the field loop flow to - // each borehole path and then normalizes the borehole-wall load history by borehole count. - Real64 const massFlowRatePerBorehole = inputs.massFlowRate; + Real64 const massFlowRatePerBorehole = inputs.massFlowRate / numBoreholes; if (massFlowRatePerBorehole <= small) { energy = 0.0; @@ -722,20 +720,20 @@ struct Model::Impl outletTemp = segments.front().outlet2(); - Real64 boreholeHeatRatePerBorehole = 0.0; + Real64 boreholeHeatRateToGroundPerBorehole = 0.0; for (auto const &seg : segments) { - boreholeHeatRatePerBorehole += seg.boreholeHeatRate(); + boreholeHeatRateToGroundPerBorehole += seg.boreholeHeatRate(); } Real64 const cp = fluidProps.cp(inputs.inletTemp); - Real64 const heatRate = inputs.massFlowRate * cp * (inputs.inletTemp - outletTemp); - Real64 const historyHeatRatePerBorehole = (config.numBoreholes > 1) ? (heatRate / numBoreholes) : boreholeHeatRatePerBorehole; + Real64 const heatRateToGround = inputs.massFlowRate * cp * (inputs.inletTemp - outletTemp); + Real64 const historyHeatRatePerBorehole = (config.numBoreholes > 1) ? (heatRateToGround / numBoreholes) : boreholeHeatRateToGroundPerBorehole; energy = historyHeatRatePerBorehole / clampMin(config.boreholeLength, small) * static_cast(inputs.timeStepSeconds); ModelStepOutputs outputs; outputs.outletTemp = outletTemp; - outputs.heatRate = heatRate; - outputs.boreholeHeatRate = boreholeHeatRatePerBorehole * numBoreholes; + outputs.heatRate = -heatRateToGround; + outputs.boreholeHeatRate = -boreholeHeatRateToGroundPerBorehole * numBoreholes; outputs.boreholeWallTemp = bhWallTemp; outputs.avgFluidTemp = 0.5 * (inputs.inletTemp + outletTemp); return outputs; @@ -774,12 +772,12 @@ struct Model::Impl Real64 const doPipe = std::max(config.pipeOuterDiameter, di + small); Real64 const re = 4.0 * massFlowRate / (mu * Constant::Pi * di); Real64 const pr = cp * mu / k; + Real64 constexpr nuLow = 4.01; Real64 nu = 0.0; if (re < 2000.0) { - nu = 4.01; + nu = nuLow; } else if (re < 4000.0) { - Real64 const nuLow = 4.01; Real64 const nuHigh = turbulentNusselt(re, pr); Real64 const sf = 1.0 / (1.0 + std::exp(-(re - 3000.0) / 150.0)); nu = (1.0 - sf) * nuLow + sf * nuHigh; diff --git a/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc index 05f4168b306..fa2cdd941dc 100644 --- a/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc +++ b/tst/EnergyPlus/unit/GroundHeatExchangers/GLHECModel.unit.cc @@ -339,6 +339,52 @@ TEST(GLHECModel, RepeatedCallsAtSameTimeDoNotAccumulate) EXPECT_NEAR(first.boreholeWallTemp, second.boreholeWallTemp, 1.0e-9); } +TEST(GLHECModel, HeatRateMatchesEnergyPlusSignConvention) +{ + auto const fluid = makeFluidFuncs(); + Model model(makeConfig(), fluid); + model.reset(15.0); + + ModelStepInputs inputs; + inputs.timeSeconds = 600; + inputs.timeStepSeconds = 600; + inputs.massFlowRate = 1.0; + inputs.inletTemp = 20.0; + inputs.farFieldGroundTemp = 15.0; + + auto const outputs = model.simulate(inputs); + + EXPECT_NEAR(outputs.heatRate, inputs.massFlowRate * fluid.cp(inputs.inletTemp) * (outputs.outletTemp - inputs.inletTemp), 1.0e-6); +} + +TEST(GLHECModel, MultiBoreholeFieldsUsePerBoreholeFlow) +{ + auto const fluid = makeFluidFuncs(); + + ModelConfig singleCfg = makeConfig(); + singleCfg.numBoreholes = 1; + + ModelConfig fieldCfg = makeConfig(); + fieldCfg.numBoreholes = 4; + + Model singleModel(singleCfg, fluid); + Model fieldModel(fieldCfg, fluid); + singleModel.reset(15.0); + fieldModel.reset(15.0); + + ModelStepInputs inputs; + inputs.timeSeconds = 600; + inputs.timeStepSeconds = 600; + inputs.massFlowRate = 1.0; + inputs.inletTemp = 20.0; + inputs.farFieldGroundTemp = 15.0; + + auto const singleOut = singleModel.simulate(inputs); + auto const fieldOut = fieldModel.simulate(inputs); + + EXPECT_GT(fieldOut.outletTemp, singleOut.outletTemp); +} + TEST(GLHECModel, AggregationTuningChangesLongTermResponse) { ModelConfig coarse = makeConfig(); @@ -386,6 +432,7 @@ TEST(GLHECModel, PrototypeTrajectoryRegressionCheckpoints) }; // Reference checkpoints from prototype GLHEC run (temps2.csv) at 600 s timesteps. + // Prototype heat rates are positive to the ground; EnergyPlus reports the opposite sign convention. std::vector const checkpoints = { {0, 16.1, 19.6839, 2995.95, 444.649, 16.1}, {1, 16.1021, 19.6861, 2983.52, 422.268, 16.5592}, @@ -439,8 +486,8 @@ TEST(GLHECModel, PrototypeTrajectoryRegressionCheckpoints) EXPECT_NEAR(outputs.boreholeWallTemp, cpRef.boreholeWallTemp, 0.5); if (step >= 500) { - EXPECT_NEAR(outputs.heatRate, cpRef.heatRate, 150.0); - EXPECT_NEAR(outputs.boreholeHeatRate, cpRef.boreholeHeatRate, 150.0); + EXPECT_NEAR(outputs.heatRate, -cpRef.heatRate, 150.0); + EXPECT_NEAR(outputs.boreholeHeatRate, -cpRef.boreholeHeatRate, 150.0); } ++nextCheckpoint; }